How to Build a Genome Analyzer You Trust

A “genome analyzer” is a pipeline: reads in FASTQ, alignment in CRAM, variants in VCF, annotation in a table, and a filtering strategy that takes roughly 4.5 million variants down to a few dozen worth a second look. The pipeline is not one program. If you have your own 30x whole-genome data, you can run the whole thing on a 32-core cloud instance for under $20 of compute in about six hours. The hard part is not running the tools. The hard part is knowing your error floor well enough to say which of those few dozen calls are real.

Here is the stack we’d use and why.

What you should have in hand

Before anything, check what your provider gave you. There is a ladder of usefulness:

• FASTQ (.fq.gz): raw reads plus base qualities. Best. You can realign to any reference, including T2T-CHM13, and re-call with a newer caller in five years. Roughly 40–60 GB per 30x genome, gzipped.
• CRAM/BAM: aligned reads. Still good. CRAM is reference-compressed and runs about 15–20 GB at 30x, versus 60–90 GB for BAM. You can re-call variants but you’re stuck with someone else’s aligner and reference unless you reconvert to FASTQ.
• VCF only: the caller’s conclusions, no evidence. You can annotate but you cannot check a single call. Treat a VCF-only delivery as a lossy summary.
• A consumer genotyping array text file (600k–1M sites): not a genome. About 0.02% of your positions, chosen for common variants, with essentially no power for rare coding variants.

Ask for FASTQ plus CRAM. Keep the FASTQ in cold object storage and work from CRAM.

The pipeline we’d run

Reference first, because it determines everything downstream. Use the GRCh38 no-ALT analysis set. Not because it’s the best assembly (CHM13 is) but because ClinVar, gnomAD, VEP caches, and every PRS weight file are keyed to GRCh38 coordinates. Running a second alignment against CHM13 is useful for segmental duplications and the acrocentric arms, but keep GRCh38 as the coordinate system of record.

Alignment:

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:L1\tSM:me\tPL:ILLUMINA\tLB:lib1' \
  GRCh38_no_alt.fa r1.fq.gz r2.fq.gz \
| samtools sort -@ 8 -m 4G -o me.bam -
samtools markdup -@ 8 me.bam me.markdup.bam
samtools view -@ 8 -T GRCh38_no_alt.fa -C -o me.cram me.markdup.bam

-K 100000000 fixes the batch size so output is deterministic regardless of thread count, which matters when you want a reproducible hash. -Y uses soft clipping for supplementary alignments, which keeps structural variant callers happy.

Small variants: run DeepVariant, not GATK HaplotypeCaller, if you’re calling a single Illumina genome.

run_deepvariant --model_type=WGS --num_shards=32 \
  --ref=GRCh38_no_alt.fa --reads=me.cram \
  --output_vcf=me.vcf.gz --output_gvcf=me.g.vcf.gz

DeepVariant’s advantage is that it was trained on the same kind of base-quality and context artifacts that reversible-terminator chemistry produces, so it handles homopolymer indels and low-complexity regions better than a hand-tuned hard filter ¹. GATK is the right answer when you have a cohort and can use joint genotyping across samples. For n=1, it isn’t.

Coverage QC, before you look at a single variant:

mosdepth -n --fast-mode --by 1000 me me.cram

Look at the fraction of the callable genome at ≥20x. Below about 90%, your “no variant here” calls are not evidence of anything. Also run verifyBamID2 for contamination and check the estimated alpha is under 0.01, and run somalier or a quick X/Y coverage ratio to confirm the sample’s reported sex matches the data. Sample swaps happen.

Then the things the small-variant caller misses:

• Structural variants: manta plus smoove/lumpy, intersected. Single-sample SV calling has a high false discovery rate. Anything you care about needs visual confirmation in IGV.
• Copy number: CNVpytor or GATK gCNV on read depth.
• Short tandem repeats: ExpansionHunter with the standard catalog. Note that FMR1 and other GC-rich repeats are poorly covered by PCR-free short reads and expansions are routinely missed.
• Pharmacogenes: PharmCAT for star-allele calling. CYP2D6 needs pangu or Cyrius because of the CYP2D7 hybrid structure. These calls describe metabolizer phenotype, nothing more. Any use in prescribing belongs with a physician and a pharmacist.

Measure your error floor

This is the step almost everyone skips. Your pipeline has a false positive rate and you should know it numerically before you interpret anything.

Download the GIAB HG002 benchmark (v4.2.1 small-variant VCF plus its high-confidence BED), run the same pipeline end to end on the HG002 reads, and score with hap.py using rtg vcfeval as the comparison engine:

hap.py HG002_benchmark.vcf.gz me_hg002.vcf.gz \
  -f HG002_highconf.bed -r GRCh38_no_alt.fa \
  --engine=vcfeval -o hg002_eval

A healthy 30x DeepVariant run lands above 99.5% F1 for SNVs and somewhat lower for indels, with indel errors concentrated in homopolymers and short tandem repeats. Two consequences follow. First, even at 99.5%, several thousand of your ~4.5 million calls are wrong. Any singleton variant you find interesting is therefore more likely a pipeline artifact than a discovery until you look at the reads. Second, error rates are not uniform: a call in a segmental duplication carries a completely different prior than a call in a unique region.

If you have parents or a sibling sequenced, call them too. Mendelian inheritance is the strongest error filter available to a single family, and it was exactly this logic that let the first family-quartet analysis separate genuine de novo mutations from sequencing error and estimate the human per-generation mutation rate ². Trio data also gives you long-range phasing for free, which tells you whether two variants in the same gene are on the same chromosome or opposite ones.

Annotation: 4.5 million down to something readable

vep --offline --cache --assembly GRCh38 --vcf --everything \
  --plugin CADD,whole_genome_SNVs.tsv.gz \
  --plugin SpliceAI,snv=spliceai_snv.hg38.vcf.gz \
  -i me.vcf.gz -o me.annot.vcf.gz --compress_output bgzip

Then filter in order: population frequency (gnomAD v4 AF, using the popmax field, not global), consequence class, and clinical annotation. A useful starting filter is AF_popmax < 0.001 and consequence in the loss-of-function set, intersected with ClinVar pathogenic/likely-pathogenic with at least two-star review status. That typically leaves a handful of heterozygous carrier findings and, in most people, nothing more.

Two failure modes worth naming. ClinVar contains many stale submissions, so a single-submitter “pathogenic” with no review is weak evidence. And a variant’s penetrance in a clinically ascertained family is not its penetrance in an unselected person who sequenced themselves for curiosity. Population-scale reanalysis has repeatedly pushed reported penetrance downward. This is the main reason clinical screening programmes still rely on targeted, validated assays with known analytic and clinical performance rather than treating a genome as a general-purpose screen ³.

If a filter surfaces something with real consequences (a BRCA-class variant, a cardiomyopathy gene or a pharmacogene with dosing implications), the next step is a CLIA-certified confirmatory test ordered through a clinician and a session with a genetic counselor. Research-grade calls are not a clinical result, and we would not act on one.

Questions people also ask

How much does it cost to analyze a genome? Sequencing and compute are now separate line items. The reagent cost of a 30x genome has fallen by roughly four orders of magnitude since the “$1000 genome” was a stated target rather than a shipping product ⁴. Compute is the cheap part: alignment plus DeepVariant on a 32-core instance runs about 5–8 hours, well under $20 spot. Storage is the recurring cost, a few dollars a month per genome in cold object storage if you keep FASTQ. Interpretation done properly is the expensive part, because it is human time.

What does a genome test tell you? Your ~4–5 million variant positions relative to the reference. From these you can derive carrier status for recessive conditions, pharmacogenomic star alleles, and a small set of well-characterized high-penetrance risk variants. You can also derive ancestry and polygenic scores. It does not tell you your gene expression, your protein levels, or what your body is doing today. The genome is the fixed layer. RNA, proteins, and metabolites are the state.

Can genome sequencing detect ADHD? No. ADHD is highly heritable at the population level but has no single causal variant with meaningful predictive power in an individual. Polygenic scores for psychiatric traits explain a small fraction of variance and transfer poorly across ancestries. Diagnosis is clinical, made by a clinician from history and behavioral assessment.

Can AI decode the human genome? AI is already load-bearing in the parts that are well-posed. DeepVariant turns variant calling into image classification, AlphaFold predicts structure from sequence, and SpliceAI predicts splice disruption better than any hand-built rule set. What no model does today is take a novel missense variant in a poorly studied gene and tell you what happens to the person carrying it. That gap is about training data, not model capacity: sequencing outpaced phenotyping, and it has since the beginning ⁵.

Footnotes

1. David Bentley, Shankar Balasubramanian, Harold Swerdlow, et al. Accurate whole human genome sequencing using reversible terminator chemistry. Nature, 2008. https://doi.org/10.1038/nature07517 ↩

2. Jared C. Roach, Gustavo Glusman, Arian F. A. Smit, et al. Analysis of Genetic Inheritance in a Family Quartet by Whole-Genome Sequencing. Science, 2010. https://doi.org/10.1126/science.1186802 ↩

3. Heidi Howard, Bartha Maria Knoppers, Martina C. Cornel, et al. Whole-genome sequencing in newborn screening? A statement on the continued importance of targeted approaches in newborn screening programmes. European Journal of Human Genetics, 2015. https://doi.org/10.1038/ejhg.2014.289 ↩

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

5. L. C. Smith, Leroy Hood. Mapping and Sequencing the Human Genome: How to Proceed. Nature Biotechnology, 1987. https://doi.org/10.1038/nbt0987-933 ↩