How to Analyze Your Own Genome: From FASTQ to a Short List of Variants

By the end of this you will have a CRAM aligned to GRCh38, a small-variant VCF with roughly 4-5 million variants, and an SV call set. You will have a measured sense of how much you can trust each of them, from running the same pipeline against a truth set. You will also have a filtered table of maybe 20-100 variants worth a second look. You need the items below.

• your raw reads (paired FASTQ or an unaligned/aligned CRAM from the sequencing provider)
• about 500 GB of free disk
• 32 GB RAM
• 16+ cores
• Docker or Singularity
• a package manager (I use micromamba)

A GPU is optional and cuts DeepVariant from hours to minutes. Budget a weekend for the first run and an evening for every run after that.

A variant is a position where your DNA differs from the reference sequence that the field agreed to use as a coordinate system. Most are single-nucleotide substitutions, some are small insertions or deletions, a smaller number are larger rearrangements. The early catalogs established the scale of common variation: 1.42 million SNPs mapped across the genome at roughly one every 1.9 kb.¹ Your personal set is mostly those common variants, with a long tail of rare ones that are the interesting part and the hard part.

1. Check what you were given

Before any compute, confirm the inputs. Sequencing providers deliver inconsistent things. You may get gzipped FASTQ, or an aligned BAM against an old reference. You may get a CRAM without the reference it was compressed against, or sometimes a VCF only.

md5sum -c checksums.md5
seqkit stats -a sample_R1.fastq.gz sample_R2.fastq.gz

Look at read length, total bases and Q30 fraction. A read length of 150 bp paired is standard. For total bases, 30x of a 3.1 Gb genome is ~93 Gbase, so ~310 M read pairs at 2x150. The quality fraction should be >90%. If you got a CRAM, find out which reference build and which specific FASTA:

samtools head input.cram | grep '^@SQ' | head -3
samtools head input.cram | grep '^@PG'

The @SQ lines carry M5 checksums for each contig. Those tell you whether you have GRCh38 with alt contigs, without, with or without HLA decoys. Get this wrong and CRAM decoding silently fails or produces garbage.

Our recommendation on reference: use GRCh38 with the analysis set that includes decoys and alt contigs, and align with alt-awareness. T2T-CHM13 is a better assembly and resolves centromeres, acrocentric arms and segmental duplications that GRCh38 gets wrong. Yet almost every annotation resource, ClinVar record and population frequency database is keyed to GRCh38 coordinates. Run CHM13 as a second pass if you care about the regions GRCh38 cannot represent. Do not make it your primary.

2. Align

If starting from FASTQ:

bwa-mem2 mem -t 16 -K 100000000 -Y \
  -R '@RG\tID:rg1\tSM:sample1\tLB:lib1\tPL:ILLUMINA' \
  GRCh38_full_analysis_set_plus_decoy_hla.fa \
  sample_R1.fastq.gz sample_R2.fastq.gz \
  | samtools sort -@4 -m 4G -o sample.bam -
samtools index sample.bam

-K 100000000 fixes the chunk size so output is deterministic regardless of thread count. -Y uses soft clipping for supplementary alignments, which SV callers want. The read group is not optional: GATK and several downstream tools will refuse to run without SM.

Mark duplicates:

samtools markdup -@8 --reference GRCh38.fa \
  <(samtools collate -@4 -O -u sample.bam | samtools fixmate -@4 -m -u - -) \
  sample.markdup.cram --write-index --output-fmt cram

CRAM against a pinned reference gets you roughly 40-50% of BAM size. Keep the exact FASTA forever.

3. Measure coverage and contamination before you call anything

Most bad variant calls trace back to something you could have caught here.

mosdepth -t 4 -n --fast-mode --by 1000 -Q 20 sample sample.markdup.cram

Read sample.mosdepth.summary.txt. You want mean autosomal depth close to what you paid for and a callable fraction (bases at ≥10x, MQ≥20) above 95%. Check chrX and chrY depth to confirm the sample’s sex chromosomes match expectation. A mismatch usually means a sample swap at the lab.

VerifyBamID --SVDPrefix resource/1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference GRCh38.fa --BamFile sample.markdup.cram

FREEMIX above 0.02 means ~2% of your reads came from someone else. At that level you start getting false heterozygous calls at low allele fraction, and any rare-variant analysis becomes unreliable. Contamination is more common than people assume and it is the single most useful QC number.

4. Call small variants

Use DeepVariant. It is more accurate than GATK HaplotypeCaller on Illumina data out of the box, particularly on indels, and it needs no recalibration step or hand-tuned filters.

docker run -v "$PWD":/data google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/GRCh38.fa \
  --reads=/data/sample.markdup.cram \
  --output_vcf=/data/sample.dv.vcf.gz \
  --output_gvcf=/data/sample.dv.g.vcf.gz \
  --num_shards=16

Emit the gVCF even if you only want the VCF. The gVCF records reference confidence at non-variant positions, which is how you distinguish “reference here” from “no data here”. That distinction matters when a variant you were looking for is absent: without a gVCF you cannot tell whether you are homozygous reference or the region had zero coverage.

If your data is Oxford Nanopore or PacBio HiFi, change --model_type and consider Clair3 as an alternative. If you have long reads, add --phase_reads for HiFi or run whatshap phase afterwards. Phasing tells you whether two heterozygous variants in the same gene sit on the same chromosome copy or opposite copies, which changes their interpretation completely.

5. Benchmark your own pipeline

This is the step people skip and it is the one that converts your VCF from a list of assertions into a measurement with error bars. Download the GIAB HG002 reads, run them through the exact pipeline above, and compare to the truth set. The current Q100 diploid assembly-based benchmark covers essentially the whole genome including previously excluded segmental duplications and the MHC, so you can measure performance in the hard regions rather than only the easy ones.²

hap.py HG002_benchmark.vcf.gz HG002.dv.vcf.gz \
  -f HG002_benchmark.bed \
  -r GRCh38.fa \
  -o hg002_vs_dv \
  --engine=vcfeval --threads 16

Expect SNV F1 around 0.995-0.999 and indel F1 several tenths of a percent lower, with the gap widening in homopolymers and tandem repeats. Now stratify, because the genome-wide number hides everything that matters:

hap.py ... --stratification GRCh38-stratifications.tsv

The stratification BED files split results by homopolymer length, segmental duplication, low mappability and GC extremes. You will find that your indel precision in homopolymers ≥10 bp is far below the headline figure. When a variant later shows up in one of those regions, you will know how much to discount it.

6. Normalize

Two VCFs can encode the same variant differently. Fix the representation before joining against any database.

bcftools norm -m -any -f GRCh38.fa -c s --atomize sample.dv.vcf.gz \
  | bcftools norm -d exact \
  | bcftools view -f PASS -Oz -o sample.norm.vcf.gz
bcftools index -t sample.norm.vcf.gz

-m -any splits multiallelic records into one line each. -f left-aligns indels. --atomize breaks MNPs into their component SNVs so they match database entries. -c s sets the REF allele from the FASTA and warns on mismatch. If you see many REF mismatches, you have the wrong reference and everything upstream needs redoing.

7. Annotate

Ensembl VEP, cached and offline. Online mode is too slow for a whole genome.

vep -i sample.norm.vcf.gz -o sample.vep.vcf.gz \
  --cache --offline --dir_cache $HOME/.vep --assembly GRCh38 \
  --vcf --compress_output bgzip --fork 16 \
  --everything --pick_allele_gene --canonical --mane \
  --custom file=gnomad.genomes.v4.1.sites.vcf.gz,short_name=gnomADg,format=vcf,type=exact,fields=AF%AF_afr%AF_amr%AF_eas%AF_nfe%AF_sas \
  --custom file=clinvar_20250601.vcf.gz,short_name=ClinVar,format=vcf,type=exact,fields=CLNSIG%CLNREVSTAT%CLNDN \
  --plugin dbNSFP,dbNSFP4.7a_grch38.gz,REVEL_score,AlphaMissense_score,CADD_phred \
  --plugin SpliceAI,snv=spliceai_scores.snv.hg38.vcf.gz,indel=spliceai_scores.indel.hg38.vcf.gz

Notes on the flags. --mane marks the MANE Select transcript, which is the one clinical labs report against, so prefer it when a variant has different consequences on different transcripts. --pick_allele_gene keeps one consequence per allele per gene instead of the dozens VEP produces by default. The --custom annotations should point at the actual gnomAD and ClinVar files rather than the versions bundled in the cache, because both move fast and the cache lags.

On the in-silico predictors: REVEL, AlphaMissense and CADD are trained on overlapping data and correlate heavily, so agreement between them is weaker evidence than it looks. SpliceAI is the more useful one because it predicts a specific mechanism (splice gain or loss at a given distance) rather than a generic pathogenicity score. A SpliceAI delta score above 0.5 near a canonical splice site is worth following up. Treat all of these as hypothesis generators, never as answers.

8. Triage down to a readable list

Filtering order matters. Start with frequency, because it removes 99% of rows.

bcftools +split-vep sample.vep.vcf.gz -f \
  '%CHROM\t%POS\t%REF\t%ALT\t%SYMBOL\t%Consequence\t%gnomADg_AF\t%ClinVar_CLNSIG\t%ClinVar_CLNREVSTAT\t%REVEL_score\t%SpliceAI_pred\t[%GT\t%GQ\t%DP]\n' \
  -d -A tab \
  -i 'FORMAT/GQ>=20 && FORMAT/DP>=10' \
  > sample.annotated.tsv

Then in pandas or duckdb, three passes:

1. Known and reviewed. Keep anything with CLNSIG containing Pathogenic or Likely_pathogenic and CLNREVSTAT of at least two stars (criteria_provided,_multiple_submitters,_no_conflicts or better). One-star and no-assertion ClinVar records include a large amount of stale and wrong material. This pass usually returns a handful of rows.
2. Rare and predicted damaging. Keep gnomADg_AF < 0.001 in every population, protein-truncating or high-SpliceAI consequence, in a gene with a known loss-of-function mechanism. Check the population-specific AFs, not just the global one: a variant at 0.0001 globally can be at 0.02 in one ancestry group and therefore benign. Reference panels are ancestry-dependent, and deep sequencing of single populations turns up large numbers of variants absent from broader catalogs.³ If your ancestry is underrepresented in gnomAD, expect this pass to return more rows and more of them to be noise.
3. Carrier and pharmacogenetic status. Run PharmCAT on the normalized VCF for star-allele calls across CYP2C19, CYP2D6, DPYD and TPMT, plus SLCO1B1 and others. CYP2D6 needs a structural-variant-aware caller (Cyrius or StellarPGx) because copy number and hybrid alleles are common there and a plain VCF cannot represent them.

Nothing in this output is a diagnosis. A pathogenic-looking variant in your VCF needs orthogonal confirmation in a CLIA/CAP laboratory and interpretation by a clinical geneticist or genetic counselor, who will weigh your family history and phenotype. Pharmacogenetic calls are for discussion with a physician or pharmacist, never for changing a medication on your own. The false-positive rate on a single unconfirmed research-grade call is high enough that acting on one is a bad idea.

9. Structural variants and copy number

Small-variant callers miss everything above ~50 bp. Run at least one SV caller and one CNV caller.

configManta.py --bam sample.markdup.cram --referenceFasta GRCh38.fa --runDir manta_out
manta_out/runWorkflow.py -j 16

cnvpytor -root sample.pytor -rd sample.markdup.cram
cnvpytor -root sample.pytor -his 1000 10000
cnvpytor -root sample.pytor -partition 1000 10000
cnvpytor -root sample.pytor -call 1000 > sample.cnv.tsv

Short-read SV calls have precision in the 60-90% range depending on type and size, far worse than SNVs. Deletions call well, insertions and inversions call badly, and anything in a segmental duplication should be assumed wrong until proven otherwise. Cross-check candidates against gnomAD-SV for frequency and look at the region in IGV before you believe anything.

10. Keep the run reproducible

Wrap the whole thing in Nextflow or Snakemake and pin every container by digest. Do not pin by tag. Record the reference FASTA md5, the VEP cache version, the gnomAD and ClinVar release dates, and the tool versions in a JSON next to the output. Re-annotate with fresh ClinVar every six months. Reinterpretation changes classifications: variants of uncertain significance get resolved in both directions over time, and the sequencing data does not change but the meaning of it does. If you did nothing else after your first run, re-running step 7 annually against current databases is the highest-value repeat work.

Common problems

Reference mismatch. bcftools norm -c s reports thousands of REF mismatches. You aligned to one FASTA and are normalizing against another, or you mixed GRCh38 with and without alt contigs. Fix the FASTA, do not force through with -c w.

chr prefix chaos. Your VCF says chr1, your annotation database says 1, and the join returns zero rows. Standardize with bcftools annotate --rename-chrs and check the contig lines in the header of every file you touch.

Depth looks fine but the callable fraction is low. Usually high duplicate rate (a low-complexity library over-amplified) or a large fraction of MQ0 reads. Check samtools flagstat for duplicate percentage. Above 15% means you paid for coverage you did not get.

Too many candidate variants after filtering. Almost always a frequency filter applied only to global AF, or a missing GQ/DP filter, or an ancestry mismatch with the reference panel.³ Add population-specific AF thresholds and require GQ≥20.

A ClinVar pathogenic hit in a gene you have never heard of. Check the review status and the submitter count first. Then check whether the associated condition is dominant or recessive and whether your genotype matches. A single heterozygous variant in a recessive gene makes you a carrier and does not make you affected. Then check whether it falls in a region your HG002 benchmark says is unreliable. Most of these evaporate. The ones that do not belong in front of a clinician.

Indels in homopolymers. If a candidate is a 1 bp indel inside a run of 8+ identical bases, your stratified hap.py numbers already told you the precision there is poor. Require long-read or capillary confirmation.

Missing variant where you expected one. Check the gVCF. Zero-depth and homozygous-reference look identical in a plain VCF and mean opposite things.

Footnotes

1. The International SNP Map Working Group, Cold Spring Harbor Laboratories:, Ravi Sachidanandam, et al. A map of human genome sequence variation containing 1.42 million single nucleotide polymorphisms. Nature, 2001. https://doi.org/10.1038/35057149 ↩

2. Nancy F. Hansen, Nathan Dwarshuis, Hyun Joo Ji, et al. A complete diploid human genome benchmark for personalized genomics. Cell, 2026. https://doi.org/10.1016/j.cell.2026.06.016 ↩

3. Masao Nagasaki, Jun Yasuda, Fumiki Katsuoka, et al. Rare variant discovery by deep whole-genome sequencing of 1,070 Japanese individuals. Nature Communications, 2015. https://doi.org/10.1038/ncomms9018 ↩ ↩²