Skip to main content

WES/WGS Single-Sample Pipeline

Diagram: Single-Sample WES/WGS Workflow​

WES/WGS single-sample workflow

Purpose​

This pipeline processes one sample at a time and produces a filtered VCF and gVCF organized for downstream tools and project QC. It automatically adapts to:

  • WES: restricted to an exome interval list
  • WGS: whole genome (no interval restriction)
WES interval source

The interval resource depends on the software stack. Legacy gatk-3.5 Bash WES uses Agilent SureSelect hg19 BED files, while current gatk-4.6 native WES uses the GATK bundle / Broad b37 exome interval list. See the FAQ for details.

What the Pipeline Does​

1. Alignment & Read Groups​

  • Align paired-end FASTQ files using BWA-MEM.
  • Add read groups (sample, library, lane, platform) required by GATK.
  • Output: lane-level BAMs with correct RG tags.

In the bundled BWA 0.7.18, the default batch target is 10,000,000 input bases per thread. Changing -t/--threads therefore changes the batch boundaries, which can alter paired-end insert-size estimation and, in a small number of cases, read placement. CBIcall fixes the batch size at -K 40000000: 10,000,000 bases multiplied by the recommended four-thread setting. This preserves the validated four-thread behavior while allowing other thread counts to produce the same alignments; -t/--threads still controls parallelism. The behavior is defined in the BWA 0.7.18 source (default batch target and effective batch size); the newer bwa-mem3 guide describes the same batching behavior.

2. Lane Merging​

  • Merge all lane BAMs for the same sample into a single BAM.
  • Ensures duplicate marking and BQSR operate on the full dataset.

3. Duplicate Marking​

  • Use GATK MarkDuplicates on the merged BAM.
  • Flags PCR/optical duplicates to prevent them from inflating support for artefactual variants.

4. Base Quality Score Recalibration (BQSR)​

  • Two-step process: BaseRecalibrator then ApplyBQSR.
  • Uses known variant databases (dbSNP, Mills, 1000G indels) to model and correct systematic base-quality errors.
  • Output: recalibrated BAM used for variant calling.

Set export_mtdna_bam: true to also write exports/mtdna/<id>-DNA_MIT.bam and its index. This small handoff file is the required input for a later mtDNA run and is preserved when cleanup_bam: true.

5. Variant Calling (HaplotypeCaller, gVCF)​

  • Run GATK HaplotypeCaller in GVCF mode (-ERC GVCF).
  • WES: uses exome intervals; WGS: full genome.
  • Output: <id>.hc.g.vcf.gz (per-sample gVCF).

6. GenotypeGVCFs (Raw VCF)​

  • Run GATK GenotypeGVCFs on the sample gVCF.
  • Output: <id>.hc.raw.vcf.gz (raw VCF with SNPs and indels).

7. Variant Quality Score Recalibration (VQSR)​

  • If there are enough variants (SNPs and indels), build VQSR models:
    • VariantRecalibrator for SNPs and indels separately.
    • Uses multiple annotations (QD, MQ, FS, MQRankSum, ReadPosRankSum).
  • Output: recalibration VCFs and tranche files.

8. Apply VQSR or Hard Filters​

  • If models exist:
    • Apply SNP VQSR.
    • Then apply INDEL VQSR.
    • Output: <id>.hc.vqsr.vcf.gz.
  • If not:
    • Skip directly to hard filters on the raw VCF or post-SNP VQSR VCF.

9. Generate Final QC VCF​

  • Run VariantFiltration with the GATK 4.6 hard filters below.
  • Output: <id>.hc.QC.vcf.gz (final QC VCF).
Filter nameExpression
LowQUALQUAL < 30.0
QD2QD < 2.0, when QD is present
FS60FS > 60.0
MQ40MQ < 40.0
MQRS-12.5MQRankSum < -12.5, when MQRankSum is present
RPRS-8ReadPosRankSum < -8.0, when ReadPosRankSum is present
QD2_indelQD < 2.0, when QD is present
FS200FS > 200.0
RPRS-20ReadPosRankSum < -20.0, when ReadPosRankSum is present

10. Coverage & Sex Determination​

  • Extract reads from the selected QC coverage contig from raw and recalibrated BAMs.
  • Compute coverage statistics in a region-first tabular file. The region defaults to chr1 and can be changed with qc_coverage_region.
  • Infer sample sex from VCF-derived depth proxies using a dedicated script. The sex file includes a # METHOD comment plus X_AUTOSOME_RATIO, X_MINUS_Y_DEPTH, THRESHOLD, and DECISION fields so the final call is auditable.
  • Outputs:
    • 03_stats/<id>.coverage.txt
    • 03_stats/<id>.sex.txt

Output Files​

FileMeaning
02_varcall/<id>.hc.g.vcf.gzPer-sample gVCF (HaplotypeCaller)
02_varcall/<id>.hc.raw.vcf.gzRaw VCF after GenotypeGVCFs
02_varcall/<id>.hc.vqsr.vcf.gzVCF after VQSR (if VQSR was applied)
02_varcall/<id>.hc.QC.vcf.gzFinal QC-filtered VCF (recommended)
03_stats/<id>.coverage.txtCoverage metrics (region, sample, mode, read counts, depth, duplicate and target percentages)
03_stats/<id>.sex.txtSex determination result
exports/mtdna/<id>-DNA_MIT.bamOptional mtDNA-only BAM for a later MToolBox run
exports/mtdna/<id>-DNA_MIT.bam.baiIndex for the optional mtDNA-only BAM
logs/<id>.logMain pipeline log

When to Use This Pipeline​

  • Standard research WES or WGS processing.
  • Generating gVCFs for cohort joint genotyping.
  • Producing filtered single-sample VCFs for downstream review or interpretation.
Workflow source files

The WES and WGS WDL entrypoints share the same implementation logic, with pipeline-specific behavior selected by the validated pipeline value.