1. GATK4 Germline WGS/WES Pipeline — Pipelines

Pipeline 1 of 16

DNA-Seq • Germline Authority: Broad Institute / nf-core/sarek

Standard Blueprint v4.5

1. GATK4 Whole-Genome / Exome Germline Variant Discovery

The international clinical and research standard for short-read sequencing data. Aligns paired-end Illumina/MGI FASTQ reads against GRCh38, performs Base Quality Score Recalibration (BQSR), discovers high-confidence germline SNPs and indels via local de novo haplotype re-assembly, and enables scalable joint multi-sample cohort genotyping.

Metric / Requirement	Specification	Recommended Target
Primary Engine	GATK 4.5.0.0 + BWA-MEM2	Nextflow (nf-core/sarek) or Snakemake
Benchmark Runtime	~8 – 12 hours (30× WGS)	32 CPU cores / 64 GB RAM
Standard Container	docker pull broadinstitute/gatk:4.5.0.0	Singularity / Docker / Apptainer
Input → Deliverable	Paired-end FASTQ.gz + Reference FASTA	Recalibrated BAM + Filtered VCF.gz

Workflow Chain: Raw FASTQ → fastp → BWA-MEM2 → MarkDuplicatesSpark → BaseRecalibrator → HaplotypeCaller (-ERC GVCF) → GenomicsDBImport → GenotypeGVCFs → VariantFiltration / VEP

Pipeline Architecture 1.1: Broad Institute GATK4 Germline Best Practices Workflow DAG with QC checkpoints and joint cohort genotyping.

⚡ Jump to Step: 1. QC & Trimming ↓ 2. Alignment ↓ 3. MarkDuplicates ↓ 4. BQSR ↓ 5. HaplotypeCaller ↓ 6. Joint Genotyping ↓ 7. Filtering & Annotation ↓

✂️

Step 1: Raw FASTQ Quality Profiling & Adapter Trimming

fastp v0.23+

Scans raw paired-end sequencing reads, auto-detects adapters, trims low-quality 3′ bases (Phred < 20), removes poly-G artifacts from NextSeq/NovaSeq, and generates interactive HTML quality metrics.

fastp \
  --in1 sample_R1.fastq.gz --in2 sample_R2.fastq.gz \
  --out1 trimmed_R1.fastq.gz --out2 trimmed_R2.fastq.gz \
  --detect_adapter_for_pe \
  --trim_poly_g \
  --cut_front --cut_tail --cut_mean_quality 20 \
  --length_required 50 \
  --thread 16 \
  --html fastp_report.html --json fastp_report.json

Key Flags: --trim_poly_g handles two-color chemistry false G artifacts; --cut_tail clips degrading read ends; --length_required 50 discards overly short inserts.

🧬

Step 2: Reference Genome Alignment & Read Group Assignment

BWA-MEM2

Aligns reads to the GRCh38 human reference using SSE/AVX-512 hardware-accelerated FM-index matching. Injects critical RG (Read Group) headers necessary for downstream GATK sample tracking.

bwa-mem2 mem -t 32 -M \
  -R "@RG\tID:FLOWCELL_L001\tSM:SAMPLE_01\tPL:ILLUMINA\tLB:LIB_WGS_1\tPU:UNIT_01" \
  /references/GRCh38.d1.vd1.fa \
  trimmed_R1.fastq.gz trimmed_R2.fastq.gz | \
  samtools view -bS -@ 8 - | \
  samtools sort -@ 8 -m 4G -o sample.sorted.bam -

samtools index -@ 8 sample.sorted.bam

Key Flags: -M flags shorter split hits as secondary for Picard compatibility; -R embeds sample identity; samtools sort -m 4G manages memory per thread.

📑

Step 3: PCR Duplicate Marking (Spark-Accelerated)

MarkDuplicatesSpark

Identifies independent read pairs that share identical 5′ alignment coordinates to prevent PCR amplification bias from skewing variant allele frequencies.

gatk --java-options "-Xmx32G -XX:+UseParallelGC" MarkDuplicatesSpark \
  -I sample.sorted.bam \
  -O sample.dedup.bam \
  -M sample.dup_metrics.txt \
  --create-output-bam-index true \
  --spark-master local[16]

Key Flags: MarkDuplicatesSpark executes 4× faster than legacy Picard by utilizing memory-mapped parallel threads; --create-output-bam-index true generates the BAI index concurrently.

⚖️

Step 4: Base Quality Score Recalibration (BQSR)

GATK BQSR

Employs machine learning to adjust reported Phred quality scores based on empirical error covariates (read group, cycle number, dinucleotide context) excluding known polymorphism sites (dbSNP & 1000G).

# Phase A: Build empirical covariate recalibration table
gatk --java-options "-Xmx16G" BaseRecalibrator \
  -R /references/GRCh38.d1.vd1.fa \
  -I sample.dedup.bam \
  --known-sites /references/Homo_sapiens_assembly38.dbsnp138.vcf.gz \
  --known-sites /references/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
  -O recal_data.table

# Phase B: Apply recalibration model to BAM
gatk --java-options "-Xmx16G" ApplyBQSR \
  -R /references/GRCh38.d1.vd1.fa \
  -I sample.dedup.bam \
  --bqsr-recal-file recal_data.table \
  -O sample.recal.bam \
  --create-output-bam-index true

Key Flags: --known-sites masks true biological variation so real variants are not penalized; ApplyBQSR outputs calibrated base qualities in the QUAL string.

🔍

Step 5: De Novo Haplotype Re-Assembly & Per-Sample gVCF Calling

HaplotypeCaller

Reconstructs active genomic regions via local De Bruijn graphs, evaluates haplotype likelihoods using PairHMM, and emits a genomic VCF (gVCF) with reference confidence records.

gatk --java-options "-Xmx24G -XX:GCTimeLimit=50 -XX:GCThresholdGCTimeLimit=80" HaplotypeCaller \
  -R /references/GRCh38.d1.vd1.fa \
  -I sample.recal.bam \
  -O sample.g.vcf.gz \
  -ERC GVCF \
  --native-pair-hmm-threads 16 \
  -L chr1,chr2,chr3,chr4,chr5,chr6,chr7,chr8,chr9,chr10,chr11,chr12,chr13,chr14,chr15,chr16,chr17,chr18,chr19,chr20,chr21,chr22,chrX,chrY,chrM

Key Flags: -ERC GVCF produces a block-compressed gVCF essential for cohort genotyping without the N+1 re-calling bottleneck; --native-pair-hmm-threads leverages vectorized AVX assembly.

🌐

Step 6: Cohort GenomicsDB Import & Joint Genotyping

GenotypeGVCFs

Combines hundreds or thousands of single-sample gVCFs into high-performance TileDB column storage and executes joint cohort genotyping across all loci simultaneously.

# Phase A: Import cohort samples into TileDB
gatk --java-options "-Xmx32G" GenomicsDBImport \
  --genomicsdb-workspace-path cohort_genomicsdb \
  -L /references/wgs_calling_regions.hg38.interval_list \
  --sample-name-map cohort_samples.tsv \
  --reader-threads 8

# Phase B: Perform joint multi-sample genotyping
gatk --java-options "-Xmx32G" GenotypeGVCFs \
  -R /references/GRCh38.d1.vd1.fa \
  -V gendb://cohort_genomicsdb \
  -O cohort_joint_raw.vcf.gz

Key Flags: gendb:// accesses the optimized multi-sample array database; eliminates the classic N+1 bottleneck when adding future clinical samples to the cohort.

🏷️

Step 7: Variant Filtration & Clinical Annotation (VEP / ClinVar)

VariantFiltration • VEP

Filters spurious calls using Broad Institute standard hard-filtering thresholds (QD < 2.0, FS > 60.0, MQ < 40.0) and annotates functional consequences, HGVS nomenclature, gnomAD allele frequencies, and ClinVar significance.

# Phase A: GATK Hard Filtering for SNPs
gatk VariantFiltration \
  -R /references/GRCh38.d1.vd1.fa \
  -V cohort_joint_raw.vcf.gz \
  --filter-expression "QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
  --filter-name "GATK_SNP_FAILED" \
  -O cohort_joint_filtered.vcf.gz

# Phase B: Ensembl Variant Effect Predictor (VEP) Clinical Annotation
vep --input_file cohort_joint_filtered.vcf.gz \
    --output_file cohort_joint_annotated.vcf.gz \
    --format vcf --vcf --compress_output bgzip \
    --species homo_sapiens --assembly GRCh38 \
    --everything --fork 16 \
    --custom /references/clinvar.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT,CLNDN

Key Flags: QD < 2.0 flags variants with low normalized confidence; --custom ClinVar annotates ACMG clinical interpretations.

🚀 Automated Nextflow (nf-core/sarek) Production Command

Execute this complete workflow across any SLURM, AWS Batch, or Google Cloud HPC cluster using the community gold-standard nf-core/sarek container pipeline:

nextflow run nf-core/sarek -r 3.4.4 \
  --input samplesheet.csv \
  --outdir ./results_gatk4_wgs/ \
  --genome GATK.GRCh38 \
  --tools haplotypecaller \
  -profile docker

🗺️ View All Pipelines Next: 2. Somatic Cancer →

← Return to Home Hub

📜 Essential Shell Scripts 📦 Practice Datasets

1. GATK4 Whole-Genome & Exome Germline Discovery