Which Genome Analysis Software to Use on Your Own WGS Data

If you have your own whole-genome sequencing data and you want to analyze it yourself, the stack is: bwa-mem2 or Dragen for alignment, DeepVariant for small variants, Manta plus CNVnator (or Dragen’s SV/CNV callers) for structural variants, Ensembl VEP for annotation, and bcftools/hail for filtering and querying. Everything except Dragen is free and open source. All of it runs on a laptop with enough disk, or better, on one cloud instance for a few hours. The hard part is not the software. It is knowing which of the ~4-5 million variants in your VCF mean anything, and the answer for most of them is that nobody knows yet.

The files you should have, and what to do first

Ask your sequencing provider for FASTQ (raw reads) and CRAM or BAM (aligned reads), not just a VCF. A VCF alone is lossy: you cannot re-call with a better model, you cannot inspect read pileups at a site you care about, and you cannot rescue regions the original pipeline filtered.

Typical sizes for a 30x human genome: two gzipped FASTQ files at ~40-50 GB total, a CRAM at ~15-20 GB (BAM is ~60-80 GB), and a gVCF at ~1-2 GB. Budget 500 GB of working disk if you plan to realign.

First thing to check is whether the data is what you paid for. Run:

samtools stats sample.cram --reference GRCh38.fa > stats.txt
mosdepth --by 1000 --fast-mode sample sample.cram

Look for mean coverage (should be ≥30x for a “30x” genome, and providers often quote raw yield rather than mapped depth), the fraction of the genome at ≥20x (want >90%), duplicate rate (<15% is normal, >25% suggests a low-complexity library), and insert size distribution. mosdepth writes a .mosdepth.summary.txt and a per-region BED in a couple of minutes.

Also confirm the reference build. GRCh38 and T2T-CHM13 are both in use, and coordinates are not interchangeable. samtools view -H sample.cram | grep '@SQ' | head shows contig names and lengths: chr1 at 248,956,422 bp is GRCh38; 248,387,328 bp is CHM13v2.0.

Alignment and variant calling

If you got FASTQ and want to run the pipeline yourself:

bwa-mem2 mem -t 32 -R '@RG\tID:s1\tSM:sample\tPL:ILLUMINA' \
  GRCh38.fa r1.fq.gz r2.fq.gz \
| samtools sort -@ 8 -m 2G -o sample.bam -
samtools index sample.bam

bwa-mem2 is roughly 2-3x faster than original bwa mem with identical output. On 32 vCPUs a 30x genome aligns in about an hour. Skip MarkDuplicates only if you are using a PCR-free library and a caller that tolerates it; otherwise run it.

For small variants (SNVs and indels), we use DeepVariant rather than GATK HaplotypeCaller for germline single samples. DeepVariant recasts variant calling as image classification over read pileups and consistently posts higher indel F1 on Genome in a Bottle truth sets, especially in homopolymers, which is where most false positives live.

run_deepvariant \
  --model_type=WGS \
  --ref=GRCh38.fa \
  --reads=sample.bam \
  --output_vcf=sample.vcf.gz \
  --output_gvcf=sample.g.vcf.gz \
  --num_shards=32

Expect ~4.5-5.0M variants against GRCh38 for one human, of which ~4-4.5M are SNVs. If you get 6M, something is wrong with your filtering or your reference. GATK is still the right choice for joint calling across a cohort or for somatic work (Mutect2), and the GATK Best Practices documentation remains the clearest published account of what each step does.

Structural variants matter more than most people expect and are systematically underreported. Short reads miss a large fraction of large insertions and repeat-mediated events; long-read platforms recover many of them. Run Manta for breakpoint-based SV, CNVnator or cnvpytor for read-depth CNV, and treat the union as a candidate list rather than a result. If you have the option, ordering long reads (PacBio HiFi or ONT) changes what you can see, since a single genome carries on the order of 20,000 SVs and short-read pipelines detect a minority of them.

Annotation: VEP, and what to do with 4.7 million rows

Ensembl VEP is the annotator we use. It resolves transcript consequences properly (canonical vs MANE Select), takes plugins, and runs offline from a cached database.

vep -i sample.vcf.gz -o sample.vep.vcf --vcf --compress_output bgzip \
  --cache --offline --dir_cache ~/.vep --assembly GRCh38 \
  --everything --mane --pick_allele_gene \
  --plugin CADD,whole_genome_SNVs.tsv.gz \
  --plugin AlphaMissense,AlphaMissense_hg38.tsv.gz \
  --custom clinvar.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT,CLNDN \
  --custom gnomad.genomes.v4.joint.sites.vcf.bgz,gnomAD,vcf,exact,0,AF,AF_popmax

--pick_allele_gene is what keeps the output to one consequence per gene per allele instead of dozens of transcript rows. --everything turns on SIFT, PolyPhen, regulatory features, frequency, and more.

Then filter. A first pass that gets you from millions to a few hundred:

bcftools view -i 'FILTER="PASS"' sample.vep.vcf.gz \
| filter_vep --filter "gnomAD_AF < 0.001 or not gnomAD_AF" \
             --filter "IMPACT is HIGH or CADD_PHRED > 25 or ClinVar_CLNSIG match pathogenic"

Two rules we hold to. First, weight ClinVar by review status: CLNREVSTAT of criteria_provided,_multiple_submitters,_no_conflicts or better (2+ stars) is a different category of evidence from a single unreviewed submission. A large share of “pathogenic” single-submitter entries do not hold up. Second, treat in-silico scores (CADD, AlphaMissense, REVEL) as ranking aids, not verdicts. They are trained on and correlated with each other, and they are not calibrated to your prior.

The intellectual limit is well documented. The bottleneck in human genomics moved from generating sequence to interpreting it, and most variants in a genome have no established functional consequence¹. Large-scale functional assays and population-scale association data are closing that gap, but the gap is still the dominant fact about your VCF.

What the software cannot do for you

Three failure modes account for most bad self-interpretation.

Low-complexity and segmental duplication regions produce false calls. The GIAB “difficult regions” BED tracks and LCR-hs38.bed exist for this reason. Intersect your candidates with them before believing anything. Genes like PMS2, SMN1, CYP2D6, and HLA have paralogs that short-read aligners cannot disambiguate, and generic pipelines get them wrong. Use specialized callers there (for example Paraphase or Cyrius) or accept that you have no call.

Penetrance is not what a database implies. Pathogenic variants ascertained through affected families look far more penetrant than the same variants found incidentally in an unselected population. If you find a “pathogenic” allele in your own data with no family history, the correct posterior is much softer than ClinVar’s label suggests.

Research-grade data is not clinical data. A variant found in your own analysis needs orthogonal confirmation in a CLIA/CAP laboratory before it means anything medically, and interpretation of a clinically significant finding belongs with a genetic counselor or medical geneticist. We do the analysis ourselves and still take findings to a clinician. Sequencing cost dropped by orders of magnitude over twenty years²³, and error and coverage characteristics of short-read platforms have been understood since the earliest population-scale efforts⁴. Neither fact makes a home-run VCF a diagnosis.

A stack we would run

• Reference: GRCh38 with the analysis set (no ALT contigs unless your caller is ALT-aware), plus T2T-CHM13 as a secondary alignment for regions GRCh38 handles badly.
• Align: bwa-mem2, or Dragen if you have access and want 30 minutes instead of 6 hours.
• Small variants: DeepVariant, gVCF output so you can distinguish “reference” from “no coverage.”
• SV/CNV: Manta plus cnvpytor, long reads if available.
• Annotate: VEP offline cache, gnomAD v4 joint frequencies, 2-star-plus ClinVar, CADD and AlphaMissense as ranking.
• Query: bcftools for one-offs, Hail if you are joining across samples or timepoints.
• Visualize: IGV on the CRAM. Always look at the reads before believing a call.

Total compute for one genome from FASTQ to annotated VCF: roughly 6-8 CPU-hours on 32 cores, under $10 on spot instances.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. 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 ↩

2. Simon T. Bennett, C.L. Barnes, Anthony J. Cox, et al. Toward the $1000 Human Genome. Pharmacogenomics, 2005. https://doi.org/10.1517/14622416.6.4.373 ↩

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

4. 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 ↩