Building a Genomics Pipeline for Your Own Whole Genome

By the end of this you will have five artifacts from one 30x whole genome.

• a CRAM aligned to GRCh38 with duplicates marked (about 15–20 GB)
• a small-variant VCF with roughly 4.5–5.0 million SNVs and indels
• a structural variant VCF with a few thousand calls
• a QC bundle (coverage, contamination, Ti/Tv and sex check)
• a flat annotated table you can query in DuckDB or hand to an agent

You need the raw data (paired FASTQ, typically 100–120 GB gzipped at 30x, or a CRAM from the sequencing provider), a Linux box with at least 32 cores and 64 GB RAM, about 1 TB of scratch disk, and Docker or Singularity. A single genome takes roughly a day of wall clock on that hardware. Nothing here is a diagnosis. A variant in a VCF is a measurement, and any clinical question about one belongs with a genetic counselor or physician who can order a validated confirmatory test.

1. Verify what you were given before you compute on it

Check the checksums the lab gave you, then look at the reads themselves. Most pipeline disasters start with an unnoticed truncated file or a lane that got concatenated twice.

md5sum -c sample.md5
seqkit stats -a -j 8 *.fastq.gz

For 30x on a 3.1 Gb genome with 2x150 reads, expect roughly 620 million read pairs. If seqkit reports an average read length well under 150 the provider trimmed for you, and if the R1 and R2 counts differ at all, the files are corrupt. Run fastqc or fastp --html once to look at per-base quality and adapter content, but do not trim WGS reads for a standard human pipeline. BWA soft-clips adapters and low-quality tails, and aggressive trimming costs you more in alignment sensitivity than it buys.

If you were given a CRAM instead, confirm which reference it points to before anything else:

samtools view -H in.cram | grep -E '^@SQ' | head -3
samtools view -H in.cram | grep -E '^@RG|^@PG'

A CRAM is only decodable with the exact reference FASTA used to compress it. Mismatched references produce silently wrong bases. They do not produce errors.

2. Pick a reference and stop revisiting the decision

Use the GRCh38 full analysis set with decoys and HLA contigs (GCA_000001405.15_GRCh38_no_alt_analysis_set_maskedGRC_exclusions_v2.fasta or the _plus_decoy_hla build). The decoy sequences absorb reads from sequences absent in the primary assembly, which cuts false positives in the regions those reads would otherwise pile into. The ALT-aware question is a real tradeoff: including ALT contigs with bwa-kit’s postalt step improves calls in HLA and a handful of other polymorphic loci, and adds a step whose output most downstream tools handle awkwardly. We run no-ALT plus decoys for everything except targeted HLA work, where we type separately with a dedicated caller.

There is a good argument for T2T-CHM13 as a second alignment, particularly for centromeric and acrocentric regions that simply do not exist in GRCh38. There is also a practical argument against making it primary: nearly every annotation database, from gnomAD to ClinVar, is coordinate-keyed to GRCh38. Align to GRCh38 for the callset you will interpret, and run CHM13 as a supplementary pass if you care about the newly resolved regions.

The broader point: two pipelines that differ in reference, aligner version, or duplicate-marking tool will not give identical VCFs. The functional equivalence specification exists precisely so that large projects can mix pipelines and still merge callsets, and it pins down aligner, reference build, and post-processing to the level of specific flags ¹. Pin your versions in a container digest and record them, so a callset you produce next year is comparable to this one.

3. Align, sort, mark duplicates and compress to CRAM

bwa-mem2 is the default. It is the same algorithm as bwa mem with a rewritten index and SIMD inner loops, roughly 2–3x faster, and the index needs about 25–30 GB resident for human.

bwa-mem2 index -p GRCh38 GRCh38_full_analysis_set_plus_decoy_hla.fa

bwa-mem2 mem -t 30 -K 100000000 -Y \
  -R '@RG\tID:S1.L001\tSM:S1\tLB:lib1\tPL:ILLUMINA\tPU:HXXXXXX.1.S1' \
  GRCh38 S1_R1.fastq.gz S1_R2.fastq.gz \
| samtools sort -@ 6 -m 2G -T /scratch/sort -o S1.sorted.bam -

Two flags matter. -K 100000000 fixes the number of bases processed per batch so output does not depend on thread count, which is what makes runs reproducible ¹. -Y uses soft clipping for supplementary alignments, which structural variant callers expect. Set the read group properly: SM is the sample identity that propagates into every VCF, and getting it wrong is the most common cause of a sample swap that survives all the way to interpretation.

Then mark duplicates. samtools markdup is fast and adequate for a single sample:

samtools collate -@ 8 -O -u S1.sorted.bam /scratch/collate \
| samtools fixmate -@ 8 -m -u - - \
| samtools sort -@ 8 -m 2G -u -T /scratch/sort2 - \
| samtools markdup -@ 8 -d 2500 -f S1.markdup.stats --write-index -O CRAM \
    --reference GRCh38_full_analysis_set_plus_decoy_hla.fa - S1.cram

-d 2500 sets the optical duplicate pixel distance for patterned flowcells (NovaSeq, NextSeq). Expect a duplicate rate of 5–15% for a normal PCR-free library. Above 25% means the library was over-amplified and your effective coverage is lower than the nominal 30x.

CRAM at reference-based compression runs about 15–20 GB for a 30x genome versus 60–80 GB for BAM, and it is lossless for the fields you care about if you keep all quality scores. Do not use quality binning to save space until you have finished calling.

4. Call small variants

Two credible paths.

DeepVariant is what we run for single-sample germline WGS. It is a CNN over pileup images, it does not need base quality score recalibration, and it produces fewer indel false positives than the alternatives on Illumina data out of the box.

docker run --rm -v "$PWD":/data google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --reads=/data/S1.cram \
  --output_vcf=/data/S1.dv.vcf.gz \
  --output_gvcf=/data/S1.dv.g.vcf.gz \
  --num_shards=30 \
  --intermediate_results_dir=/scratch/dv

On 30 CPU cores this is on the order of a day, with call_variants dominating. A GPU cuts that step to well under an hour. Always write the gVCF even for one sample. It is what lets you distinguish “reference at this position” from “no data at this position,” which is the difference between a meaningful negative and an artifact.

The GATK path is HaplotypeCaller -ERC GVCF, then GenotypeGVCFs, then VQSR or hard filters. The GATK framework’s design assumptions are worth understanding even if you use DeepVariant: local reassembly around candidate sites, recalibration of the base qualities that the instrument reports, and a variant-level quality model trained on known-site overlap rather than per-site thresholds ². That framework also makes joint calling across many samples materially more sensitive at low coverage, which is why cohort projects use it ². For one genome at 30x, joint calling buys you nothing and DeepVariant is simpler.

Whichever you use, normalize before you compare anything:

bcftools norm -m -any -f GRCh38.fa -Oz -o S1.norm.vcf.gz S1.dv.vcf.gz
bcftools index -t S1.norm.vcf.gz

Left-aligned, biallelic-split records are the only form in which two VCFs can be meaningfully intersected. Skip this and your concordance numbers will be wrong by several percent, entirely in the indels.

5. Call structural variants and repeats separately

Small-variant callers ignore anything longer than about 50 bp. Structural variants are a small count with outsized effect, and they need their own tools.

configManta.py --bam S1.cram --referenceFasta GRCh38.fa --runDir manta_S1
manta_S1/runWorkflow.py -j 24

Manta takes an hour or two and emits deletions, duplications, insertions, inversions, and breakends. Expect a few thousand PASS calls per genome, dominated by deletions in the 100 bp to 10 kb range. Precision for short-read SV calling is meaningfully worse than for SNVs, so treat any single call as a hypothesis requiring orthogonal confirmation.

For tandem repeat expansions, ExpansionHunter with the standard variant catalog genotypes the known pathogenic loci. It is genotyping specific pre-specified repeats, not discovering new ones, and repeat lengths near or beyond the read length are estimates with wide confidence intervals.

6. QC the alignment and the callset, with actual thresholds

Coverage:

mosdepth --fast-mode --no-per-base --by 1000 -t 8 S1 S1.cram

Check mean autosomal depth (aim for 30x or better), and the fraction of the callable genome at 10x or more (should exceed 0.95 for a good 30x PCR-free library). Coverage varies across the genome and is not uniform. GC-rich first exons and promoters drop, and high-identity segmental duplications collect reads with MAPQ 0 that most callers discard. Hard regions like these are a known and persistent limitation of short-read genome sequencing rather than a defect in your run ³.

Contamination:

verifybamid2 --SVDPrefix resource/1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference GRCh38.fa --BamFile S1.cram --NumThread 8

FREEMIX above 0.02 means roughly 2% of reads come from another individual, which inflates heterozygous calls and low-allele-fraction artifacts. Above 0.03, re-sequence rather than interpret.

Callset shape:

bcftools stats -F GRCh38.fa -s - S1.norm.vcf.gz > S1.stats
plot-vcfstats -p plots/ S1.stats

Reference numbers for a 30x human genome start with 4.5–5.0 million SNVs plus indels in total. The transition/transversion ratio is near 2.0–2.1 genome-wide, and the heterozygous/homozygous-alternate ratio is around 1.4–1.7 depending on ancestry. Over 95% of SNVs are present in dbSNP. A Ti/Tv of 1.7 means false positives. A variant count of 3.5 million means a European-ancestry sample undercalled, or a coverage problem. A count of 6 million means contamination or a broken filter.

Sex check: compute mean depth on chrX non-PAR and chrY from the mosdepth output, and the heterozygosity rate on chrX non-PAR. XY samples show chrX at roughly half autosomal depth with near-zero chrX heterozygosity. Mismatches against the reported sex are the fastest way to catch a sample swap.

If you have a genotyping array on the same person, intersect it. Concordance at shared sites should exceed 99%. That single number catches almost every class of pipeline error at once.

7. Annotate

VEP with a local cache, offline, is the workhorse.

vep --offline --cache --dir_cache /data/vep --assembly GRCh38 \
  --fasta GRCh38.fa --species homo_sapiens \
  --input_file S1.norm.vcf.gz --format vcf \
  --vcf --compress_output bgzip --output_file S1.vep.vcf.gz \
  --fork 16 --buffer_size 50000 \
  --everything --pick_allele_gene --canonical --mane \
  --custom file=/data/gnomad.genomes.v4.1.sites.vcf.gz,short_name=gnomADg,format=vcf,type=exact,fields=AF%AF_grpmax%nhomalt \
  --plugin CADD,snv=/data/whole_genome_SNVs.tsv.gz

Then layer ClinVar as a separate annotation pass so you can update it independently. ClinVar changes weekly and you will want to re-run just that step:

bcftools annotate -a clinvar_20260901.vcf.gz \
  -c INFO/CLNSIG,INFO/CLNREVSTAT,INFO/CLNDN \
  -Oz -o S1.vep.clinvar.vcf.gz S1.vep.vcf.gz

Pay attention to CLNREVSTAT. A “pathogenic” call with no assertion criteria and one submitter is a different object from one reviewed by an expert panel, and the review status is the field that tells you which you have.

8. Flatten it into something queryable

A VCF with a 40-field CSQ string is not a good substrate for exploration or for an agent. Flatten it once:

bcftools +split-vep S1.vep.clinvar.vcf.gz \
  -d -s worst \
  -f '%CHROM\t%POS\t%REF\t%ALT\t%FILTER\t%SYMBOL\t%Gene\t%Consequence\t%IMPACT\t%MANE_SELECT\t%HGVSc\t%HGVSp\t%gnomADg_AF\t%CADD_PHRED\t%CLNSIG\t%CLNREVSTAT[\t%GT\t%DP\t%GQ]\n' \
  -i 'FILTER="PASS"' \
| gzip > S1.flat.tsv.gz

Load that into DuckDB and write it back out as Parquet. A 5-million-row table with 19 columns is around 150–250 MB as Parquet, small enough to keep in memory and fast enough that interactive queries return in under a second.

CREATE TABLE v AS SELECT * FROM read_csv('S1.flat.tsv.gz', delim='\t', header=false);
COPY v TO 'S1.parquet' (FORMAT PARQUET, COMPRESSION ZSTD);

From there the useful queries are simple: high-impact consequences in MANE Select transcripts with gnomAD allele frequency below 0.001, ClinVar entries with two or more stars, or pharmacogenomic star alleles from a dedicated caller like PharmCAT. Handing the Parquet file plus the QC summary to an agent works well because the agent can compute over the whole table instead of guessing from a text excerpt. The interpretation limits do not change: a rare, high-impact-annotated variant in your file is a hypothesis about a measurement, and acting on it requires a clinician and an orthogonal confirmatory assay from an accredited lab.

Common problems

Reference and contig-naming mismatch. A chr1 VCF intersected with a 1 annotation file returns zero overlaps and no error. Check with bcftools view -h | grep contig on both files before you merge anything.

Read group missing or wrong. Callers will run, and the sample name in the output will be SAMPLE or the file basename. Set SM at alignment time. Do not leave it until afterward.

Un-normalized VCFs. The same indel written two valid ways will not match between your callset and ClinVar. Always run bcftools norm -m -any -f ref before comparison or annotation.

Low-complexity and segmental-duplication regions. Genes with high-identity paralogs (SMN1/SMN2, PMS2/PMS2CL, and the like) produce systematically unreliable short-read calls, and this is a property of the data rather than the pipeline ³. Flag those regions and treat calls inside them as uninformative without a targeted assay.

Silent disk exhaustion. DeepVariant’s intermediate make_examples output for a 30x genome runs to hundreds of gigabytes. Point --intermediate_results_dir at real scratch, not /tmp.

Comparing across pipeline versions without controlling for them. Aligner version, reference build, and duplicate-marking tool all shift the callset by a percent or more. If you re-sequence in two years, either re-run the old data through the new pipeline or use an explicitly equivalent configuration ¹.

Assuming more coverage fixes everything. Sequencing throughput and cost per genome have improved by orders of magnitude, and the remaining hard problems are mapping ambiguity, repeat structure, and interpretation rather than raw read count ⁴. Going from 30x to 60x improves indel sensitivity modestly and does nothing for a region where every read maps at MAPQ 0.

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. Mark A. DePristo, Eric Banks, Ryan Poplin, et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics, 2011. https://doi.org/10.1038/ng.806 ↩ ↩²

3. 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 ↩ ↩²

4. Sang Tae Park, Jayoung Kim. Trends in Next-Generation Sequencing and a New Era for Whole Genome Sequencing. International Neurourology Journal, 2016. https://doi.org/10.5213/inj.1632742.371 ↩