How to Analyze Your Own Whole Genome, Step by Step

By the end of this you will have a CRAM aligned to GRCh38 and a normalized and annotated VCF with roughly 4.5–5 million small variants. You will have a structural variant call set and QC numbers that tell you whether to trust any of it. You will also have a filtered shortlist of a few hundred variants worth a human reading. You need your raw reads, a machine with 32+ cores and 64 GB RAM, about 1 TB of scratch disk, and Docker or Singularity. The raw reads are paired-end FASTQ, ~110 GB gzipped for a 30x human genome, or a delivered CRAM. Budget 12–24 hours of wall clock on CPU, or 1–2 hours if you have an A100-class GPU and use DeepVariant’s GPU build. Everything below is measurement and interpretation. Nothing here tells you what to do about a result, and anything that looks medically consequential belongs in front of a clinical geneticist or genetic counselor working from a CLIA/CAP-validated report.

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

Get the checksums and the read metadata first. A truncated FASTQ transfer will silently produce a worse variant call set. It will not produce an error.

md5sum -c raw.md5
seqkit stats -a -j 8 sample_R1.fastq.gz sample_R2.fastq.gz

Check three things. Read count: a 30x genome at 2x150 bp needs roughly 600 million read pairs (2 × 600M × 150 bp ≈ 180 Gbp, of which ~90 Gbp lands on a 3.1 Gbp genome after duplicates and unmapped reads). Read length: 150 bp is standard on NovaSeq, and 100 bp costs you mappability in repeats. Q30 fraction: you want >90% of bases at Q30 or better.

If your provider delivered a CRAM instead, confirm which reference it was aligned to before anything else. A CRAM stores reference-relative deltas, so it is not self-contained. Reading it with the wrong FASTA gives you garbage or a hard failure.

samtools head sample.cram | grep '^@SQ' | head -3
# look at the M5 and UR tags; they identify the exact reference FASTA

2. Trim and screen the reads

We use fastp because it does adapter detection, quality trimming, and a QC report in one pass at close to I/O speed.

fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
      -o trim_R1.fastq.gz -O trim_R2.fastq.gz \
      --detect_adapter_for_pe \
      --cut_tail --cut_tail_mean_quality 20 \
      --length_required 50 \
      --thread 16 \
      --json fastp.json --html fastp.html

Resist heavy trimming. Aggressive quality trimming biases the ends of reads and costs you more in coverage than it gains in accuracy, because BWA soft-clips low-quality ends anyway and both GATK and DeepVariant model base qualities. If fastp reports more than 10% of reads discarded, stop and look at the library rather than tuning flags.

3. Align to GRCh38, no-alt, with bwa-mem2

Use the 1000 Genomes GRCh38 no-alt analysis set (GCA_000001405.15_GRCh38_no_alt_analysis_set.fna). It includes decoy contigs and excludes alt haplotypes, which matters because alt-aware handling is a common source of MAPQ=0 reads in HLA, KIR, and other polymorphic regions if your downstream caller is not alt-aware.

bwa-mem2 index ref/GRCh38_no_alt.fna

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:S1\tSM:SAMPLE\tPL:ILLUMINA\tLB:lib1' \
  ref/GRCh38_no_alt.fna trim_R1.fastq.gz trim_R2.fastq.gz \
| samtools sort -@ 8 -m 2G -o sorted.bam -

samtools markdup -@ 8 --reference ref/GRCh38_no_alt.fna \
  sorted.bam sample.cram --output-fmt CRAM
samtools index sample.cram

-K 100000000 fixes the chunk size so output is deterministic regardless of thread count, which you want if you plan to re-run and diff. -Y uses soft clipping for supplementary alignments, which structural variant callers need.

On GRCh38 versus T2T-CHM13: GRCh38 is where the ecosystem lives. ClinVar, gnomAD, VEP caches, and PGS Catalog scoring files are all keyed to it. CHM13 gives you better calls in centromeric and acrocentric regions and closes ~200 Mb of previously unresolvable sequence, but you will do the annotation work yourself. We align to GRCh38 as the primary and treat CHM13 as a second alignment for specific questions.

Expect 4–6 hours on 32 cores, and a final CRAM of 15–20 GB at 30x. Keep the CRAM, delete the BAM and the trimmed FASTQ. The CRAM plus the reference reproduces everything downstream.

4. Run QC before you call anything

Four numbers decide whether the run is usable.

mosdepth -n -t 4 --fast-mode --by 1000 qc/sample sample.cram
samtools stats -@ 8 --reference ref/GRCh38_no_alt.fna sample.cram > qc/stats.txt
verifybamid2 --SVDPrefix resource/1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference ref/GRCh38_no_alt.fna --BamFile sample.cram --NumThread 8
somalier extract -d qc/ --sites sites.hg38.vcf.gz -f ref/GRCh38_no_alt.fna sample.cram

The targets are mean autosomal coverage ≥28x and ≥95% of the callable genome at ≥10x. Duplicate rate should be under 15%, and PCR-free libraries usually land near 5%. The verifyBamID2 FREEMIX value should be under 0.01. FREEMIX above 0.02 means contamination or a sample mixture, and every heterozygous call becomes suspect. Somalier gives you inferred sex from X/Y coverage and heterozygosity, which is the cheapest sample-swap check available.

5. Call small variants

We use DeepVariant for germline SNVs and indels on Illumina data. It is a single tool, needs no base quality score recalibration, and in our hands gives cleaner indel calls in homopolymers than a stock GATK HaplotypeCaller run.

docker run -v "$PWD":/data google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/ref/GRCh38_no_alt.fna \
  --reads=/data/sample.cram \
  --output_vcf=/data/sample.vcf.gz \
  --output_gvcf=/data/sample.g.vcf.gz \
  --num_shards=32

Always write the gVCF as well as the VCF. The gVCF records reference confidence per block, which is the only way to distinguish “reference at this position” from “no data at this position.” When you later want to know whether a specific ClinVar position was covered, the gVCF answers it and the VCF does not.

Sanity check the output:

bcftools stats sample.vcf.gz | grep -E '^SN|^TSTV'

For a 30x WGS you should see roughly 4.3–5.0 million variants total (higher for genomes with more African ancestry, since GRCh38 is European-biased), about 3.6–4.2 million SNVs, 600–900k indels, and a genome-wide Ti/Tv near 2.0–2.1. A Ti/Tv of 1.7 or below means false-positive SNVs are leaking in. Het/hom-alt ratio near 1.5–1.6 is typical. Large deviations usually trace back to contamination or an unusual reference build.

6. Call structural variants, CNVs, and repeat expansions

Small variant callers miss a class of variation that carries real functional weight. Run three more tools.

# Structural variants, Manta
configManta.py --bam sample.cram --referenceFasta ref/GRCh38_no_alt.fna \
  --runDir manta/ && manta/runWorkflow.py -j 32

# Deletions/duplications with split-read and discordant-pair evidence
smoove call --name SAMPLE --fasta ref/GRCh38_no_alt.fna \
  --processes 16 --genotype sample.cram

# Short tandem repeats at known pathogenic loci
ExpansionHunter --reads sample.cram --reference ref/GRCh38_no_alt.fna \
  --variant-catalog variant_catalog/grch38/variant_catalog.json \
  --output-prefix eh/SAMPLE

Short-read SV calling is the weakest link in the whole pipeline. Expect 5,000–10,000 SV calls, a substantial false-positive rate for insertions, and poor resolution inside segmental duplications. Treat any SV you care about as a hypothesis requiring orthogonal confirmation (long reads, or a clinical lab assay). ExpansionHunter is genuinely useful but is limited to loci in the catalog and loses accuracy when the expansion exceeds the read length by a wide margin.

7. Normalize, then annotate

Normalization first, always. Indel representation is not unique, and an unnormalized variant will fail to match ClinVar or gnomAD even when it is the same variant.

bcftools norm -m -any -f ref/GRCh38_no_alt.fna -Ou sample.vcf.gz \
| bcftools norm -d exact -Oz -o sample.norm.vcf.gz
bcftools index -t sample.norm.vcf.gz

Then annotate with VEP against a local cache. Remote annotation over the REST API on 4.5 million variants will take days.

vep -i sample.norm.vcf.gz -o sample.vep.vcf.gz \
  --vcf --compress_output bgzip --fork 16 \
  --cache --dir_cache $HOME/.vep --offline --assembly GRCh38 \
  --everything --pick_allele_gene \
  --custom gnomad.genomes.v4.1.sites.vcf.gz,gnomADg,vcf,exact,0,AF,AF_grpmax,nhomalt \
  --custom clinvar_20250601.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT,CLNDN \
  --plugin CADD,whole_genome_SNVs.tsv.gz \
  --plugin SpliceAI,snv=spliceai_scores.raw.snv.hg38.vcf.gz,indel=spliceai_scores.raw.indel.hg38.vcf.gz

Pin the versions. Write down the VEP release, the Ensembl gene set, the gnomAD version, and the ClinVar release date. ClinVar reclassifies variants continuously, and a variant called “likely pathogenic” in one release can be “uncertain significance” eighteen months later. This churn is one of the persistent obstacles in clinical genome work, and the reason a genome result is a dated snapshot rather than a permanent fact ¹.

8. Filter down to something a human can read

Four and a half million variants is not a result. Tiering is where the analysis happens. We run three independent passes over the same annotated VCF rather than one merged filter.

Pass A, rare and predicted-damaging in known disease genes:

bcftools +split-vep sample.vep.vcf.gz -f '%CHROM\t%POS\t%REF\t%ALT\t%SYMBOL\t%Consequence\t%gnomADg_AF\t%CADD_PHRED\t%ClinVar_CLNSIG\t[%GT]\n' \
  -i 'FILTER="PASS" && (gnomADg_AF<0.001 || gnomADg_AF=".") && (Consequence~"frameshift" || Consequence~"stop_gained" || Consequence~"splice_acceptor" || Consequence~"splice_donor" || CADD_PHRED>25)' \
  -d -A tab > tierA.tsv

This typically returns 300–800 rows genome-wide. Most are in genes with no established disease association.

Pass B, ClinVar pathogenic and likely pathogenic with at least one-star review status, regardless of predicted consequence. This catches missense variants with established evidence that no in silico score would flag. Filter on CLNREVSTAT to exclude no_assertion_criteria_provided entries, which are the main source of noise.

Pass C, a fixed gene list. We use the ACMG secondary findings list (v3.2, 81 genes) as the anchor because it is the set with the strongest actionability evidence, plus any genes relevant to a documented family history. Restrict with a BED file rather than a SYMBOL match, so you catch variants in UTRs and deep intronic positions.

A structured, multi-pass approach with explicit inclusion criteria is what clinical labs do, and the published best practices for genome interpretation describe the same shape: define the analysis scope up front, apply variant-level and gene-level filters separately, and record what was excluded ².

Then read the shortlist by hand. Interpretation is the expensive part. The Stanford group’s early work on clinically annotating whole genomes found that assessing potentially reportable variants took roughly an hour of expert time per genome per reviewer even after automated filtering, and reviewers disagreed with each other on a meaningful fraction of calls ³. Nothing since has made that step cheap.

Two calibrations to hold onto. A typical healthy person carries roughly 100 loss-of-function variants and a handful of homozygous knockouts, none of which mean anything on their own. And the great majority of what survives your filters will be variants of uncertain significance, where the truthful answer is that nobody knows. Population-scale sequencing of healthy adults consistently produces more uncertainty than actionable findings ⁴.

9. Polygenic scores and ancestry

These are computed from common variants and are a different kind of statement from a single rare variant. Use the PGS Catalog’s pipeline rather than rolling your own scoring, because allele orientation and strand ambiguity will bite you.

plink2 --vcf sample.norm.vcf.gz --make-pgen --out sample
nextflow run pgscatalog/pgsc_calc -profile docker \
  --input samplesheet.csv --target_build GRCh38 \
  --pgs_id PGS000765 --run_ancestry ref_panel/

Run the ancestry adjustment step. Polygenic scores are calibrated in the population they were derived in, and an unadjusted score from a European-derived model applied to a non-European genome produces a number with no defensible interpretation. Even adjusted, a PGS is a population-level percentile. It is not a personal prediction.

10. Store it so future-you can re-run it

Keep the CRAM, the gVCF, the reference FASTA and its index, and a text file recording every tool version and container digest. Discard intermediate BAMs. Re-annotation is cheap and the calls are stable, so the yearly operation is as follows. Pull the current ClinVar and gnomAD releases, re-run steps 7 and 8 on the stored VCF, and diff the shortlist against last year’s. Reanalysis on updated databases produces new findings in previously negative genomes at a non-trivial rate, which is a real argument for treating your genome as a dataset you revisit rather than a report you file ⁵.

Common problems

The CRAM will not open. You have the wrong reference. Match the M5 checksums in the @SQ header lines against samtools faidx output on your candidate FASTA. GRCh38 with alts, without alts, with and without decoys, and the Verily and hs38DH variants are all different files.

Ti/Tv is 1.6 and you have 6 million variants. Either contamination (check FREEMIX), or you are counting non-PASS records. bcftools stats includes filtered variants unless you pass -f PASS or pre-filter.

A known family variant is absent from the VCF. Check the gVCF at that exact coordinate before concluding anything. Low coverage, a segmental duplication, or a GC-extreme region can all produce no-call rather than reference. Then check that the coordinate is in the right build. GRCh37 and GRCh38 positions differ for essentially every variant, and lifting over with CrossMap has its own failure modes at inversion breakpoints.

Indels do not match ClinVar. You skipped normalization, or you normalized against a different reference than ClinVar used. Run bcftools norm on both sides.

Coverage looks fine but chrX calls are strange. Diploid calling on a haploid X produces spurious heterozygotes. DeepVariant handles this if you supply the correct sex via the --haploid_contigs and --par_regions_bed options.

Everything worked and the shortlist is 400 rows of uncertainty. That is the expected outcome. It is not a pipeline failure. Sequencing a healthy genome has a low yield of clearly actionable findings relative to the volume of ambiguous ones, and the resulting follow-up testing and counseling carry real costs ⁶. The cases where a genome changes clinical management are concentrated in people with an unexplained phenotype or a strong family history ⁷. If something on your list looks consequential, the next step is a clinical lab confirming it on a validated assay, and a genetic counselor interpreting it against your actual family history. Do not act on a research-grade call.

Footnotes

1. Kelly E. Ormond, Matthew T. Wheeler, Louanne Hudgins, et al. Challenges in the clinical application of whole-genome sequencing. The Lancet, 2010. https://doi.org/10.1016/s0140-6736(10)60599-5 ↩

2. Christina Austin‐Tse, Vaidehi Jobanputra, Denise Perry, et al. Best practices for the interpretation and reporting of clinical whole genome sequencing. npj Genomic Medicine, 2022. https://doi.org/10.1038/s41525-022-00295-z ↩

3. Frederick E. Dewey, Megan E. Grove, Cuiping Pan, et al. Clinical Interpretation and Implications of Whole-Genome Sequencing. JAMA, 2014. https://doi.org/10.1001/jama.2014.1717 ↩

4. Noralane M. Lindor, Stephen N. Thibodeau, Wylie Burke. Whole-Genome Sequencing in Healthy People. Mayo Clinic Proceedings, 2017. https://doi.org/10.1016/j.mayocp.2016.10.019 ↩

5. Petar Brlek, Luka Bulić, Matea Bračić, et al. Implementing Whole Genome Sequencing (WGS) in Clinical Practice: Advantages, Challenges, and Future Perspectives. Cells, 2024. https://doi.org/10.3390/cells13060504 ↩

6. Katharina Schwarze, James Buchanan, Jenny C. Taylor, et al. Are whole-exome and whole-genome sequencing approaches cost-effective? A systematic review of the literature. Genetics in Medicine, 2018. https://doi.org/10.1038/gim.2017.247 ↩

7. Frederik Otzen Bagger, Line Borgwardt, Andreas Sand Jespersen, et al. Whole genome sequencing in clinical practice. BMC Medical Genomics, 2024. https://doi.org/10.1186/s12920-024-01795-w ↩