How to Run a Genome Analysis Pipeline on Your Own WGS Data

By the end of this you will have a CRAM aligned to GRCh38, a small-variant VCF with roughly 4–5 million SNVs and indels, a structural-variant VCF, a quality report you can defend, and an annotated table you can query. You need the raw data (FASTQ or an unaligned CRAM) and a machine with at least 16 cores and 64 GB RAM. You also need about 500 GB of free disk for a 30x genome, and Docker or Singularity. Everything below runs on one box. If you would rather not manage the box, the same steps run on nf-core/sarek with a cloud executor, and we say where.

A pipeline in bioinformatics is just a directed graph of command-line tools with file formats as edges: FASTQ to BAM/CRAM to VCF to annotated TSV. The interesting part is not the graph, it is that different graphs give different answers on the same reads. Two pipelines can disagree on 1–2% of variant calls, which is tens of thousands of positions in a genome. The disagreement concentrates in exactly the regions people care about: indels, homopolymers and segmental duplications. The Functional Equivalence specification was written to make large projects agree on these choices, and it is the reason we pin versions and reference builds instead of using whatever apt gives us.¹

1. Verify what you were given

Before compute, audit the delivery. You want to know read length, insert size, lane structure, and whether the file has been touched.

md5sum -c checksums.md5
zcat sample_R1.fastq.gz | head -8
zcat sample_R1.fastq.gz | awk 'NR%4==2{print length($0)}' | head -1000 | sort | uniq -c

A 30x PE150 human genome is roughly 900 million read pairs, around 100–120 GB of gzipped FASTQ across both mates. If your files total 40 GB, you were sold something closer to 10x, and low coverage silently degrades heterozygote sensitivity long before it breaks anything visibly. Check the header line for instrument and flowcell so you can build read groups:

@A00123:456:HXXXXDSXX:1:1101:1234:5678 1:N:0:ATCACGAT

That gives you ID=HXXXXDSXX.1, PU=HXXXXDSXX.1.ATCACGAT, PL=ILLUMINA, LB=lib1 and SM=yourid. Read groups matter: duplicate marking and base quality recalibration both partition on them, and a missing RG tag will make GATK refuse to run at step 5.

2. Pin the reference and the container set

Use the GRCh38 full analysis set with decoys and without ALT contigs unless you are prepared to run an ALT-aware post-processing step. The hs38DH build (GCA_000001405.15_GRCh38_no_alt_analysis_set.fna plus the HLA/decoy additions) is what the Functional Equivalence spec settled on, and matching it means your VCF can be compared to public callsets without liftover.¹

samtools faidx GRCh38_full_analysis_set_plus_decoy_hla.fa
bwa-mem2 index GRCh38_full_analysis_set_plus_decoy_hla.fa   # ~80 GB peak RAM, ~1 hour
gatk CreateSequenceDictionary -R GRCh38_full_analysis_set_plus_decoy_hla.fa

Pull containers with explicit tags, never :latest. Record the digests in a versions.txt you keep next to the outputs. We treat that file as part of the result. Six months from now, the difference between DeepVariant 1.6.1 and 1.8.0 on your indel calls will matter more than you expect.

3. Align, mark duplicates, write CRAM

bwa-mem2 is the default we reach for. It is a drop-in for BWA-MEM with the same alignment output and roughly 2x the throughput. DRAGMAP (the open reimplementation of the DRAGEN mapper) is the alternative and is faster still with slightly better indel behavior in repeats, at the cost of a larger hash-table index and less community tooling around it.

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:HXXXXDSXX.1\tSM:yourid\tLB:lib1\tPL:ILLUMINA\tPU:HXXXXDSXX.1.ATCACGAT' \
  GRCh38_full_analysis_set_plus_decoy_hla.fa sample_R1.fastq.gz sample_R2.fastq.gz \
| samtools sort -@ 8 -m 4G -o sample.sorted.bam -

-K 100000000 fixes the batch size so output is deterministic regardless of thread count. -Y uses soft clipping for supplementary alignments, which structural-variant callers need. Both flags are part of the FE spec.¹

Then mark duplicates and compress:

gatk MarkDuplicatesSpark -I sample.sorted.bam -O sample.md.bam \
  --optical-duplicate-pixel-distance 2500 --conf 'spark.executor.cores=16'

samtools view -C -T GRCh38_full_analysis_set_plus_decoy_hla.fa \
  -@ 16 -o sample.cram sample.md.bam
samtools index sample.cram

Expect 4–8 hours of wall clock on 32 cores for alignment, and a final CRAM of 15–20 GB at 30x. Delete the intermediate BAMs once the CRAM’s read count matches (samtools flagstat on both). Keep the FASTQs or the unaligned CRAM forever. Reference builds change, callers improve, and reads are the only thing you cannot regenerate.

4. Quality-control the alignment before you trust any variant

Four numbers decide whether the rest of the run is meaningful.

mosdepth --by 1000 --fast-mode -t 8 sample sample.cram
samtools stats -@ 8 sample.cram > sample.stats
VerifyBamID --SVDPrefix resource/1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference GRCh38_full_analysis_set_plus_decoy_hla.fa --BamFile sample.cram

• Mean autosomal coverage: want ≥30x for confident heterozygote calls. Below 20x, het sensitivity falls off sharply because you need multiple reads supporting the alt allele.
• Callable fraction: from sample.mosdepth.global.dist.txt, the fraction of the genome at ≥20x should be above 0.90.
• Contamination: FREEMIX from VerifyBamID2 should be under 0.02. Above that, you will see an excess of low-allele-fraction “het” calls that are someone else’s DNA.
• Duplicate rate: under 15% for PCR-free libraries. A 35% duplicate rate means your effective coverage is much lower than the raw number suggests.

Also run a sex check from chrX/chrY coverage ratio and compare it to what you expect. It catches sample swaps, which are the single most common failure in delivered sequencing data and are invisible downstream.

5. Call small variants

Two reasonable choices. We use DeepVariant for a single germline sample and reach for GATK only when we need joint genotyping across a cohort.

run_deepvariant \
  --model_type=WGS \
  --ref=GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --reads=sample.cram \
  --output_vcf=sample.dv.vcf.gz \
  --output_gvcf=sample.dv.g.vcf.gz \
  --num_shards=32

DeepVariant skips base quality score recalibration and local realignment entirely. The CNN operates on pileup images. That removes two steps and two sources of divergence between pipelines. On GIAB benchmarks it reaches SNV F1 above 99.5% and beats GATK on indels, particularly in homopolymer runs, which is where most false positives live.²

The GATK route, if you want it:

gatk HaplotypeCaller -R $REF -I sample.cram -O sample.g.vcf.gz -ERC GVCF \
  --native-pair-hmm-threads 8
gatk GenotypeGVCFs -R $REF -V sample.g.vcf.gz -O sample.raw.vcf.gz

Then filter. For a single sample, VQSR is unreliable because it needs enough variants to fit the Gaussian mixture, so use hard filters or CNNScoreVariants instead. Emit a gVCF either way. The gVCF records reference confidence at non-variant positions, which is the only way to distinguish “reference here” from “no data here” when you later want to check whether a specific position was covered.

Sanity numbers for a single 30x human genome against GRCh38: 4.0–5.0 million total variants, Ti/Tv ratio around 2.0–2.1 genome-wide (2.9–3.0 restricted to coding), het/hom ratio near 1.5–2.0 depending on ancestry. If Ti/Tv drops to 1.5, you have a false-positive problem.

6. Benchmark the pipeline without stopping at the sample

You cannot validate calls on your own genome because there is no truth set for it. You validate the pipeline by running the identical commands on GIAB HG002, for which a truth VCF and high-confidence BED exist.

hap.py HG002_GRCh38_truth.vcf.gz HG002.dv.vcf.gz \
  -f HG002_GRCh38_highconf.bed \
  -r GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --engine=vcfeval -o hg002_bench --threads 16

Read the stratified output. Do not rely on the headline F1. Performance inside the high-confidence regions tells you little about the segmental duplications and tandem repeats that sit outside them, and a meaningful fraction of medically relevant genes overlap those hard regions.² Do this once, keep hg002_bench.summary.csv, and rerun it whenever you bump a tool version. It is the only empirical claim you can make about your own numbers.

7. Structural variants, repeats, and mitochondria

Short-read SV calling is the weakest part of any short-read pipeline, but ignoring it means missing a whole class of variation. Run at least two callers and take the intersection for anything you act on.

configManta.py --bam sample.cram --referenceFasta $REF --runDir manta_out
manta_out/runWorkflow.py -j 16

smoove call --name sample --fasta $REF -p 16 --genotype sample.cram

ExpansionHunter --reads sample.cram --reference $REF \
  --variant-catalog variant_catalog.json --output-prefix sample.eh

Manta is fast and good at deletions and duplications above ~50 bp. smoove (Lumpy plus svtyper) catches a partly different set. ExpansionHunter genotypes known repeat expansions from a catalog, which short-read callers otherwise miss completely. Expect a few thousand SV calls and a high false-positive rate in the small-deletion range. Mitochondrial calls need their own path (GATK’s Mutect2-based mitochondrial mode) because the copy number is high and heteroplasmy means allele fractions are continuous, not 0/0.5/1.

8. Annotate

Annotation is where a VCF becomes something you can query. Use VEP with a local cache so it is reproducible and offline.

vep -i sample.dv.vcf.gz -o sample.vep.vcf.gz --vcf --compress_output bgzip \
  --cache --offline --dir_cache /refs/vep --assembly GRCh38 \
  --everything --pick_allele_gene --fork 16 \
  --custom /refs/gnomad.genomes.v4.1.sites.vcf.gz,gnomADg,vcf,exact,0,AF,AF_popmax \
  --plugin CADD,/refs/whole_genome_SNVs.tsv.gz \
  --plugin SpliceAI,snv=/refs/spliceai_scores.raw.snv.hg38.vcf.gz,indel=/refs/spliceai_scores.raw.indel.hg38.vcf.gz \
  --plugin dbNSFP,/refs/dbNSFP4.7a_grch38.gz,REVEL_score,AlphaMissense_score

Then flatten to a table you can load into DuckDB or pandas:

bcftools +split-vep sample.vep.vcf.gz -d -f \
  '%CHROM\t%POS\t%REF\t%ALT\t%Consequence\t%SYMBOL\t%gnomADg_AF\t%CADD_PHRED\t%SpliceAI_pred_DS_AG[\t%GT\t%DP]\n' \
  -A tab > sample.annot.tsv

The frequency filter does most of the work. Of 4.5 million variants, roughly 20,000–30,000 are protein-altering, and only a few hundred are both rare (gnomAD popmax AF < 0.001) and predicted damaging. Population frequency is the strongest single filter you have, and the large reference panels are what made individual genome interpretation tractable at all.³

Layer ClinVar on top with bcftools annotate, and treat the review status field as load-bearing. A “Pathogenic” call with one submitter and no assertion criteria is a different object from a three-star expert-panel classification.

Anything you find here is a hypothesis. It is not yet a result. Research-grade sequencing is not diagnostic. A variant that looks pathogenic needs orthogonal confirmation in a CLIA/CAP laboratory and interpretation with a genetics clinician or certified genetic counselor, who will weigh your family history and phenotype in ways a VCF cannot. The same applies to pharmacogenomic star alleles from PharmCAT: they are informative, and they are not a basis for changing anything on your own.

9. Add the RNA layer

DNA tells you what is possible. RNA tells you what is happening in the tissue sampled. For bulk RNA-seq from whole blood or PBMCs, we run two quantifications in parallel because they answer different questions.

# Transcript-level, fast, for expression tables
salmon quant -i gencode_v45_index -l A -1 rna_R1.fq.gz -2 rna_R2.fq.gz \
  -p 16 --validateMappings --gcBias --seqBias -o salmon_out

# Genome alignment, for splicing and allele-specific expression
STAR --genomeDir star_idx_gencode_v45 --readFilesIn rna_R1.fq.gz rna_R2.fq.gz \
  --readFilesCommand zcat --runThreadN 16 --outSAMtype BAM SortedByCoordinate \
  --quantMode GeneCounts --waspOutputMode vW --varVCFfile sample.het.vcf

Salmon gives TPMs in minutes. STAR gives you a BAM you can inspect at a locus, junction counts in SJ.out.tab, and WASP-filtered reads for allele-specific expression once you supply your heterozygous sites. ASE is the most interesting thing a single person can do with paired DNA and RNA: it tests whether one haplotype is expressed less than the other at a het site, which is direct evidence that a regulatory variant does something. Linking non-coding variants to expression effects is the central problem in interpreting the 98% of the genome that is not protein-coding, and expression data is the most direct handle an individual has on it.³

For a single sample there is no differential expression. DESeq2 needs replicates. What you can do with one timepoint is check expression rank against a reference panel like GTEx whole blood, and look for genes with zero counts that should be expressed. You can also estimate cell composition by deconvolution (CIBERSORTx or a simple marker-based approach), since most apparent expression change in blood is cell-proportion change. With repeated timepoints of your own, you become your own control, and that is when the data starts being worth the storage.

10. Wrap it so you can run it again

Once the steps work, move them into Nextflow or WDL. nf-core/sarek covers alignment through variant calling with sensible defaults, container pinning, and resumability:

nextflow run nf-core/sarek -r 3.4.4 -profile docker \
  --input samplesheet.csv --genome GATK.GRCh38 \
  --tools deepvariant,manta,vep --outdir results

The reason to bother is -resume and provenance. A 30x genome is 12–20 hours of compute. Having a crash at hour 14 restart from the last completed process rather than from FASTQ pays for the migration on the first failure.

Common problems

bwa-mem2 index gets OOM-killed. The index build needs roughly 80 GB of RAM for GRCh38. Build it once on a large instance, or download a prebuilt index, or fall back to classic bwa index which needs about 6 GB.

Chromosome naming mismatch. chr1 versus 1 between your reference, your dbSNP file, and your annotation cache. Symptom: zero annotated variants and no error message. Fix with bcftools annotate --rename-chrs chr_map.txt and check with bcftools view -h | grep contig | head.

Duplicate rate above 30%. Usually a PCR-amplified library with low input DNA. Effective coverage is the raw coverage times (1 − dup rate), so a 30x run at 35% duplicates is really 20x. You cannot fix this computationally.

Variant count far from 4–5 million. Under 3.5 million usually means over-filtering or low coverage. Over 6 million usually means contamination or a reference mismatch. Check FREEMIX first.

Ti/Tv near 1.0 in a subset of calls. Classic signature of false positives clustered in low-complexity regions. Intersect with the GIAB low-complexity BED and see whether the bad calls disappear.

Two pipelines, two answers. Expected. Aligner choice, duplicate-marking implementation, and caller account for most of it, which is why the Functional Equivalence work specified exact tools and flags rather than just output formats.¹ Pin your versions, record them, and compare like to like.

Assuming the annotated table is a finding. Predicted-damaging is a prior. It is not a conclusion. Most rare, damaging-looking variants in a healthy person’s genome do nothing detectable. Take anything that looks serious to a clinician with the raw data in hand.

Footnotes

1. Allison Regier, Yossi Farjoun, David E. Larson, et al. Functional equivalence of genome sequencing analysis pipelines enables harmonized variant calling across human genetics projects. Nature Communications, 2018. https://doi.org/10.1038/s41467-018-06159-4 ↩ ↩² ↩³ ↩⁴

2. Daniel C. Koboldt, Li Ding, Elaine R. Mardis, et al. Challenges of sequencing human genomes. Briefings in Bioinformatics, 2010. https://doi.org/10.1093/bib/bbq016 ↩ ↩²

3. Tuuli Lappalainen, Alexandra J. Scott, Margot Brandt, et al. Genomic Analysis in the Age of Human Genome Sequencing. Cell, 2019. https://doi.org/10.1016/j.cell.2019.02.032 ↩ ↩²