bio-workflows-fastq-to-variants

What this skill does

## Version Compatibility

Reference examples tested with: BWA-MEM2 2.2.1+, Ensembl VEP 111+, GATK 4.5+, bcftools 1.19+, fastp 0.23+, samtools 1.19+

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `<tool> --version` then `<tool> --help` to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

# FASTQ to Variants Workflow

**"Call variants from my whole-genome or exome FASTQ files"** -> Orchestrate fastp QC, BWA-MEM2 alignment, duplicate marking, BQSR, GATK HaplotypeCaller variant calling, and VQSR/hard filtering to produce filtered VCF output.

Complete pipeline from raw DNA sequencing FASTQ files to filtered variant calls.

## Workflow Overview

```
FASTQ files
    |
    v
[1. QC & Trimming] -----> fastp
    |
    v
[2. Alignment] ---------> bwa-mem2
    |
    v
[3. BAM Processing] ----> sort, markdup, index
    |
    v
[4. Variant Calling] ---> bcftools (primary) or GATK
    |
    v
[5. Filtering] ---------> Quality filters
    |
    v
Filtered VCF
```

## Primary Path: BWA + bcftools

### Step 1: Quality Control with fastp

```bash
# Single sample
fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
    -o sample_R1.trimmed.fq.gz -O sample_R2.trimmed.fq.gz \
    --detect_adapter_for_pe \
    --qualified_quality_phred 20 \
    --length_required 50 \
    --html sample_fastp.html

# Batch processing
for sample in sample1 sample2 sample3; do
    fastp -i ${sample}_R1.fastq.gz -I ${sample}_R2.fastq.gz \
        -o trimmed/${sample}_R1.fq.gz -O trimmed/${sample}_R2.fq.gz \
        --detect_adapter_for_pe \
        --html qc/${sample}_fastp.html
done
```

**QC Checkpoint 1:** Check fastp reports
- Q30 bases >85% (DNA typically higher quality than RNA)
- Adapter content <1%
- No unusual GC distribution

### Step 2: BWA-MEM2 Alignment

```bash
# Index reference (once)
bwa-mem2 index reference.fa

# Align with read group info
for sample in sample1 sample2 sample3; do
    bwa-mem2 mem -t 8 \
        -R "@RG\tID:${sample}\tSM:${sample}\tPL:ILLUMINA\tLB:lib1" \
        reference.fa \
        trimmed/${sample}_R1.fq.gz \
        trimmed/${sample}_R2.fq.gz \
    | samtools view -bS - > aligned/${sample}.bam
done
```

**QC Checkpoint 2:** Check alignment stats
```bash
samtools flagstat aligned/${sample}.bam
```
- Mapped reads >95%
- Properly paired >90%

### Step 3: BAM Processing

```bash
for sample in sample1 sample2 sample3; do
    # Sort by coordinate
    samtools sort -@ 8 -o aligned/${sample}.sorted.bam aligned/${sample}.bam

    # Mark duplicates (samtools method)
    samtools fixmate -m aligned/${sample}.sorted.bam - | \
        samtools sort -@ 8 - | \
        samtools markdup -@ 8 - aligned/${sample}.markdup.bam

    # Index
    samtools index aligned/${sample}.markdup.bam

    # Cleanup intermediate
    rm aligned/${sample}.bam aligned/${sample}.sorted.bam
done
```

**QC Checkpoint 3:** Check duplication rate
```bash
samtools flagstat aligned/${sample}.markdup.bam | grep "duplicates"
```
- WGS: <30% duplicates
- Exome/targeted: <50% duplicates

### Step 4: Variant Calling with bcftools

```bash
# Single sample calling
bcftools mpileup -Ou -f reference.fa aligned/sample1.markdup.bam | \
    bcftools call -mv -Oz -o variants/sample1.vcf.gz

# Multi-sample calling (joint calling)
bcftools mpileup -Ou -f reference.fa \
    aligned/sample1.markdup.bam \
    aligned/sample2.markdup.bam \
    aligned/sample3.markdup.bam | \
    bcftools call -mv -Oz -o variants/cohort.vcf.gz

bcftools index variants/cohort.vcf.gz
```

### Step 5: Variant Filtering

```bash
# Basic quality filter
bcftools filter -Oz \
    -e 'QUAL<20 || DP<10 || MQ<30' \
    -o variants/cohort.filtered.vcf.gz \
    variants/cohort.vcf.gz

# More stringent filter
bcftools filter -Oz \
    -e 'QUAL<30 || DP<10 || DP>200 || MQ<40 || MQB<0.1' \
    -s "LowQual" \
    -o variants/cohort.filtered.vcf.gz \
    variants/cohort.vcf.gz

# Stats
bcftools stats variants/cohort.filtered.vcf.gz > variants/vcf_stats.txt
```

**QC Checkpoint 4:** Check variant stats
- Ti/Tv ratio ~2.1 for whole genome
- Ti/Tv ratio ~2.8-3.0 for exome
- >95% overlap with dbSNP for known sites

## Alternative Path: BWA + GATK HaplotypeCaller

### Step 4 Alternative: GATK Variant Calling (DRAGEN Mode -- Recommended)

```bash
gatk CreateSequenceDictionary -R reference.fa
samtools faidx reference.fa

# DRAGEN-GATK mode: no BQSR needed, improved STR handling and QUAL calibration
gatk HaplotypeCaller \
    -R reference.fa \
    -I aligned/sample1.markdup.bam \
    -O variants/sample1.g.vcf.gz \
    -ERC GVCF \
    --dragen-mode
```

### Step 4 Alternative: GATK Standard Mode (with BQSR)

```bash
gatk BaseRecalibrator \
    -R reference.fa \
    -I aligned/sample1.markdup.bam \
    --known-sites dbsnp.vcf.gz \
    -O recal_data.table

gatk ApplyBQSR \
    -R reference.fa \
    -I aligned/sample1.markdup.bam \
    --bqsr-recal-file recal_data.table \
    -O aligned/sample1.recal.bam

gatk HaplotypeCaller \
    -R reference.fa \
    -I aligned/sample1.recal.bam \
    -O variants/sample1.g.vcf.gz \
    -ERC GVCF

# Joint genotyping (for multiple samples)
gatk GenomicsDBImport \
    -V variants/sample1.g.vcf.gz \
    -V variants/sample2.g.vcf.gz \
    -V variants/sample3.g.vcf.gz \
    --genomicsdb-workspace-path genomicsdb \
    -L intervals.bed

gatk GenotypeGVCFs \
    -R reference.fa \
    -V gendb://genomicsdb \
    -O variants/cohort.vcf.gz
```

### Step 5 Alternative: GATK Variant Filtering

```bash
# Hard filtering (for small cohorts)
gatk VariantFiltration \
    -R reference.fa \
    -V variants/cohort.vcf.gz \
    --filter-expression "QD < 2.0" --filter-name "LowQD" \
    --filter-expression "FS > 60.0" --filter-name "HighFS" \
    --filter-expression "MQ < 40.0" --filter-name "LowMQ" \
    --filter-expression "MQRankSum < -12.5" --filter-name "LowMQRS" \
    --filter-expression "ReadPosRankSum < -8.0" --filter-name "LowRPRS" \
    -O variants/cohort.filtered.vcf.gz

# VQSR (for large cohorts >30 samples)
gatk VariantRecalibrator \
    -R reference.fa \
    -V variants/cohort.vcf.gz \
    --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap.vcf.gz \
    --resource:omni,known=false,training=true,truth=false,prior=12.0 omni.vcf.gz \
    --resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G.vcf.gz \
    --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp.vcf.gz \
    -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
    -mode SNP \
    -O cohort.snp.recal \
    --tranches-file cohort.snp.tranches

gatk ApplyVQSR \
    -R reference.fa \
    -V variants/cohort.vcf.gz \
    -O variants/cohort.vqsr.vcf.gz \
    --recal-file cohort.snp.recal \
    --tranches-file cohort.snp.tranches \
    -mode SNP \
    --truth-sensitivity-filter-level 99.5
```

## Parameter Recommendations

| Step | Parameter | WGS | Exome/Targeted |
|------|-----------|-----|----------------|
| bwa-mem2 | -t | 8-16 | 8 |
| samtools markdup | - | Required | Required |
| bcftools mpileup | -d | 250 (default) | 1000 |
| bcftools mpileup | -q | 20 | 20 |
| bcftools filter | QUAL | >20 | >30 |
| bcftools filter | DP | >10, <2x mean | >20 |
| GATK | intervals | - | Target BED |

## Choosing a Variant Caller

| Criterion | bcftools | GATK HaplotypeCaller | DeepVariant |
|-----------|----------|---------------------|-------------|
| Speed | Fastest | Moderate | Slowest (without GPU) |
| SNP accuracy | Good | Good | Highest |
| Indel accuracy | Moderate (weak in homopolymers) | Good (local assembly) | Highest |
| Joint calling | Basic multi-sample | GVCF workflow (scales to 100k+) | GLnexus (scales to 250k+) |
| Best for | Quick analysis, non-model organisms | Clinical pipelines, GATK ecosystem | Highest accuracy, GPU available |
| VQSR | Not supported | For large cohorts (>30 samples) | Not needed (use GLnexus) |
| DRAGEN mode | N/A | Recommend