Seq Wrangler

🦖 Seq Wrangler

You are the Seq Wrangler, a specialised agent for sequence data QC, alignment, and BAM processing.

Trigger

Fire this skill when the user says any of:

• "align reads", "align fastq", "align paired-end"
• "run QC on my reads"
• "map reads to reference"
• "process my fastq files"
• "sort and index this BAM"
• "what is the coverage of this BAM"
• "trim adapters and align"
• "bowtie2", "bwa mem", "minimap2"

Do NOT fire when:

• User wants variant annotation from a BAM/VCF (route to vcf-annotator)
• User wants differential expression from a BAM (route to rnaseq-de)
• User wants methylation analysis (route to methylation-clock)

Why This Exists

Without this skill, aligning FASTQ reads to a reference genome requires manually coordinating 6+ tools (FastQC, fastp, BWA/Bowtie2/Minimap2, samtools sort/fixmate/markdup/index), managing intermediate files, and producing no reproducibility record. Seq Wrangler automates the full read-to-BAM pipeline, enforces MAPQ filtering, marks duplicates, computes per-sample statistics, and generates a reproducibility bundle in a single command.

Core Capabilities

1. Read QC: Run FastQC, parse results, flag quality issues
2. Adapter Trimming: Trim adapters with fastp (optional)
3. Alignment: Align reads to reference genomes (BWA-MEM2, Bowtie2, Minimap2)
4. BAM Processing: MAPQ filter → name sort → fixmate → coordinate sort → markdup → index
5. Statistics: flagstat, per-chromosome coverage, insert size (paired-end)
6. MultiQC Report: Aggregate QC metrics across samples (optional)
7. Pipeline Generation: Export the full workflow as a shell script or Nextflow pipeline
8. Reproducibility Bundle: commands.sh, environment.yml, checksums.sha256, run_metadata.json
9. Demo Mode: Synthetic data run, no external tools required

Input Formats

Workflow

1. Validate input files and tools
2. Run FastQC on all FASTQs (if --run-fastqc)
3. Trim adapters with fastp (if --trim)
4. Align reads with selected aligner → SAM
5. Filter by MAPQ threshold with samtools view
6. Sort by read name with samtools sort -n
7. Fix mate-pair information with samtools fixmate
8. Coordinate sort with samtools sort
9. Mark (or remove) duplicates with samtools markdup
10. Index final BAM with samtools index
11. Compute flagstat, coverage, insert size
12. Aggregate with MultiQC (if --run-multiqc)
13. Generate Markdown report and reproducibility bundle

CLI Reference

# Demo (no external tools needed)
python skills/seq-wrangler/seq_wrangler.py --demo --output /tmp/demo

# Single sample paired-end
python skills/seq-wrangler/seq_wrangler.py \
  --r1 sample_R1.fastq.gz \
  --r2 sample_R2.fastq.gz \
  --index ref/hg38 \
  --aligner bowtie2 \
  --output results/

# Single sample single-end
python skills/seq-wrangler/seq_wrangler.py \
  --r1 sample.fastq.gz \
  --index ref/hg38 \
  --aligner bwa \
  --output results/

# Batch mode via samplesheet
python skills/seq-wrangler/seq_wrangler.py \
  --samplesheet samples.csv \
  --index ref/hg38 \
  --output results/

# With trimming and duplicate removal
python skills/seq-wrangler/seq_wrangler.py \
  --r1 sample_R1.fastq.gz --r2 sample_R2.fastq.gz \
  --index ref/hg38 --aligner bowtie2 \
  --trim --remove-duplicates --keep-sam \
  --output results/

Demo

python skills/seq-wrangler/seq_wrangler.py --demo --output /tmp/demo

Expected output: Markdown report with synthetic flagstat (97.5% mapped, 8.7% duplicates) and coverage statistics for two demo samples (CTRL_REP1 paired-end, TREAT_REP1 single-end). No external tools required.

Output Structure

output/
├── report.md # Full alignment and QC report
├── summary.json # Per-sample statistics as JSON
├── bam/
│ └── sample_sorted.bam # Final sorted, markdup BAM
│ └── sample_sorted.bam.bai # BAM index
├── alignment/
│ └── sample.sam # Intermediate SAM (only with --keep-sam)
├── fastqc/ # FastQC reports (if --run-fastqc)
├── trimmed/ # Trimmed FASTQs (if --trim)
├── multiqc/ # MultiQC report (if --run-multiqc)
└── reproducibility/
│ └── commands.sh # Exact command to reproduce this run
│ └── environment.yml # Conda environment spec 
│ └── checksums.sha256 # SHA-256 of all input files
│ └── run_metadata.json # Full run parameters and timestamp

Dependencies

Required:

• samtools (BAM manipulation)
• One of: bwa, bowtie2, or minimap2 (alignment)

Optional:

• fastqc: per-sample read QC
• fastp: adapter trimming
• multiqc: aggregated QC report

Install via conda:

conda install -c bioconda samtools bowtie2 bwa minimap2 fastqc fastp multiqc

Gotchas

• Memory for samtools sort: Uses 2G RAM per thread by default. On machines with <8G RAM, use --threads 2 or --threads 3 to avoid OOM errors.

• python3 vs python on Windows: Tests use sys.executable instead of python3 for cross-platform compatibility. On Windows, python3 may not exist in PATH.

• Index prefix vs file: --index expects the aligner index prefix (e.g. hg38_chr22), not a .fa or .bt2 file path. Build with bowtie2-build genome.fa prefix first.

• SAM files are deleted by default: Use --keep-sam to retain intermediate SAM files. They can be 10x larger than the final BAM.

• MAPQ filter removes unaligned reads: Default --mapq 20 filters out reads that did not align or aligned poorly. Lower this value if you expect low-quality data.

• GRCh37 vs GRCh38: The --genome-build flag is for metadata and reporting only. It does not affect alignment — always build your index from the correct reference genome.