How to Analyze Your Own Whole-Genome Sequencing Data

By the end of this you will have: a CRAM aligned to GRCh38 (about 40–60 GB for a 30x genome), a per-sample gVCF and a filtered VCF with roughly 4.5 million SNVs and 600,000–900,000 indels, structural variant calls from Manta, repeat expansion genotypes from ExpansionHunter, a VEP-annotated TSV you can query, and a QC report with coverage, contamination, Ti/Tv, and duplicate rate. You need the raw FASTQs (2x150 bp, PCR-free library, ~100 GB gzipped for 30x), a Linux box with 16+ cores and 64 GB RAM, about 2 TB of free disk, and Docker or Singularity. Alignment takes 4–8 hours on 32 cores. DeepVariant takes 2–5 hours on CPU, or under an hour on a GPU. If you only have a consumer array or an “interpreted report” PDF, stop here and go get the FASTQs or CRAM. Reports are downstream opinions about a data file you should hold yourself.

1. Verify the raw data before you spend compute on it

Check the md5s your provider gave you, then confirm the read count and length match the claimed coverage.

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

Coverage math: (num_reads × read_length) / 3.1e9. For 30x with 2x150 you need roughly 310 million read pairs. If seqkit reports 180 million pairs, you have a 17x genome, and no amount of downstream cleverness fixes that. Also run fastqc or falco and look at two plots specifically: per-base N content (a spike means a bad cycle) and the duplication estimate (a PCR-free library should read under 10%).

Confirm the library is PCR-free. PCR-amplified WGS has GC-biased coverage that drops out first-exon GC-rich regions and inflates duplicates, and it distorts the allele balance DeepVariant relies on.

2. Pick a reference and stick to it

Use the GRCh38 no-alt analysis set with decoys and HLA contigs (GCA_000001405.15_GRCh38_no_alt_analysis_set.fna). Not hg19, not the full ALT-containing build unless you are running ALT-aware post-processing.

wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
gunzip GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
samtools faidx GCA_000001405.15_GRCh38_no_alt_analysis_set.fna
bwa-mem2 index GCA_000001405.15_GRCh38_no_alt_analysis_set.fna   # ~80 GB RAM, ~1 hr

Why no-alt: ALT contigs are near-duplicates of primary chromosomes, so reads from those regions get MAPQ 0 and disappear from variant calls unless the caller is ALT-aware. Decoys (chrUn_KN707*, EBV) soak up reads that would otherwise mismap into real genes.

T2T-CHM13v2.0 is a better reference in the strict sense: it closes the gaps and improves calling in segmental duplications and acrocentric arms. The problem is annotation. Almost every clinical resource, gnomAD frequency file, and ClinVar record you will want to join against is on GRCh38 coordinates. We align to GRCh38 as the primary. If a specific question involves a hard region (SMN1/SMN2, the HLA locus or KIR), we realign that region to CHM13 separately.

3. Align, mark duplicates, write CRAM

One pipe, no intermediate BAM. The -K 100000000 flag fixes the batch size so the output is deterministic across thread counts, and -Y uses soft clipping for supplementary alignments, which the SV callers want.

REF=GCA_000001405.15_GRCh38_no_alt_analysis_set.fna
RG='@RG\tID:L001\tSM:sample1\tPL:ILLUMINA\tPU:HXXXX.1.ATCACG\tLB:lib1'

bwa-mem2 mem -t 30 -K 100000000 -Y -R "$RG" "$REF" R1.fastq.gz R2.fastq.gz \
  | samtools fixmate -m -u -@ 4 - - \
  | samtools sort -u -@ 8 -m 3G -T /tmp/sort - \
  | samtools markdup -@ 8 --reference "$REF" --output-fmt CRAM \
      --write-index -f markdup.stats - sample1.cram

CRAM at reference-compressed default level is roughly 40–60 GB for a 30x genome versus 90–110 GB for BAM, and it is lossless with respect to the alignment. Keep the reference FASTA forever. A CRAM without its reference is close to unreadable.

If you have multiple lanes, align each lane with its own @RG and merge before markdup, or use samtools merge then markdup. Do not concatenate FASTQs from different flowcells under one read group, you lose the ability to diagnose a bad lane.

4. Run QC and decide whether to proceed

Four numbers matter before you look at a single variant.

mosdepth -n --fast-mode --by 1000 -t 4 --fasta "$REF" qc/sample1 sample1.cram
samtools stats -@ 8 --reference "$REF" sample1.cram > qc/sample1.stats
VerifyBamID --SVDPrefix resource/1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference "$REF" --BamFile sample1.cram --NumThread 8

• Mean coverage: target ≥30x. mosdepth’s *.mosdepth.summary.txt gives per-chromosome means.
• Callable fraction: from *.mosdepth.global.dist.txt, the fraction of the genome at ≥10x should be above 0.95. Below 0.90 and you will be chasing false negatives all day.
• Duplicate rate: from markdup.stats. PCR-free should be under 10%, often 2–5%.
• Contamination: VerifyBamID2 FREEMIX. Under 0.01 is fine. Above 0.03 and you will see inflated heterozygosity and spurious low-VAF calls. Stop and get the sample re-run.

Also do a sex check from mosdepth per-chromosome means: chrX at roughly half autosomal coverage with chrY present indicates XY, chrX at autosomal level with chrY near zero indicates XX. A mismatch against expectation usually means a sample swap at the lab. An interesting biological finding is not the usual explanation.

5. Call small variants with DeepVariant

singularity run -B /data:/data docker://google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/"$REF" \
  --reads=/data/sample1.cram \
  --output_vcf=/data/sample1.vcf.gz \
  --output_gvcf=/data/sample1.g.vcf.gz \
  --num_shards=30 \
  --intermediate_results_dir=/data/tmp_dv

We use DeepVariant over GATK HaplotypeCaller for single samples. It needs no BQSR step, no hard-filter tuning, and its indel precision in homopolymers is better out of the box. The tradeoff: DeepVariant’s model is trained on specific chemistries, so if your data is an unusual read length or an older chemistry, validate against a truth set before trusting it. GATK’s advantage shows up in joint calling across a cohort, which is not your situation with one genome.

Emit the gVCF even if you never joint-call. It records reference confidence at every position, which is the only way to distinguish “no variant here” from “no data here” later. That distinction is the single most common interpretive error people make with their own genome.

Validate against GIAB if you can. If your provider sequenced a control, or you can get HG002 data from the same platform, run hap.py against the v4.2.1 truth VCF and confidence BED. Expect SNV F1 above 0.995 and indel F1 above 0.99 on a good 30x run.

6. Sanity-check the variant set

bcftools view -f PASS sample1.vcf.gz -Oz -o sample1.pass.vcf.gz && bcftools index -t sample1.pass.vcf.gz
bcftools stats -F "$REF" sample1.pass.vcf.gz > qc/sample1.bcfstats

What to look for:

• Total PASS variants: 4.3–5.0 million. African-ancestry genomes land near the top, European near 4.4–4.6 million. Under 4 million means coverage or filtering problems. Over 5.5 million means contamination.
• Ti/Tv: 2.0–2.1 genome-wide, 2.9–3.1 restricted to coding exons. A genome-wide Ti/Tv below 1.9 means false positives leaking through.
• Het/hom-alt ratio: roughly 1.3–1.6 for European ancestry, higher for African. A high ratio with a high FREEMIX confirms contamination.
• Novel fraction: the proportion of your variants absent from dbSNP should be around 1–3%. Deep sequencing of a population not well represented in reference panels turns up far more apparently novel rare variants, which is a property of the panels rather than of you ¹.

If you have sequenced relatives, Mendelian inconsistency rate is the strongest single-number quality metric available. The original family-quartet analysis used exactly this to separate true de novo events from sequencing error, and it works the same way on your data ².

7. Structural variants, CNVs, and repeat expansions

Small-variant calling covers maybe 90% of the base pairs that differ between you and the reference, but a much smaller share of the total sequence affected. Skipping SVs leaves real biology on the table ³.

# Manta: deletions, duplications, inversions, translocations
configManta.py --bam sample1.cram --referenceFasta "$REF" --runDir manta_out
manta_out/runWorkflow.py -j 16

# Repeat expansions at known loci (HTT, FMR1, C9orf72, ATXN1-3, RFC1, ...)
ExpansionHunter --reads sample1.cram --reference "$REF" \
  --variant-catalog variant_catalog/grch38/variant_catalog.json \
  --output-prefix eh_sample1

# CYP2D6 copy number and hybrid alleles
Cyrius/star_caller.py --manifest cram_manifest.txt --genome 38 --prefix cyrius --outDir cyrius_out

Short reads cap what you can see. A 30x Illumina genome resolves deletions well, duplications less well, balanced inversions and insertions poorly, and anything inside a segmental duplication barely at all. ExpansionHunter genotypes repeats up to roughly the fragment length reliably and gives a lower bound beyond that, which means a called expansion is a signal to get a clinical PCR or long-read assay, never a conclusion on its own.

Visualize before you believe. samplot plot -n sample1 -b sample1.cram -c chr7 -s 100000 -e 120000 -t DEL -o del.png takes ten seconds and kills most false positives.

8. Annotate

vep -i sample1.pass.vcf.gz -o sample1.vep.vcf.gz --vcf --compress_output bgzip \
  --cache --offline --dir_cache $HOME/.vep --assembly GRCh38 --fork 8 \
  --everything --pick_allele_gene --buffer_size 20000 \
  --custom file=clinvar_20250601.vcf.gz,short_name=ClinVar,format=vcf,type=exact,fields=CLNSIG%CLNREVSTAT%CLNDN \
  --plugin CADD,snv=whole_genome_SNVs.tsv.gz,indels=gnomad.genomes.r4.0.indel.tsv.gz \
  --plugin SpliceAI,snv=spliceai_scores.raw.snv.hg38.vcf.gz,indel=spliceai_scores.raw.indel.hg38.vcf.gz

Then filter with slivar, which is much faster than parsing VEP output in Python:

slivar expr --vcf sample1.vep.vcf.gz --ped sample1.ped \
  -g gnomad.hg38.v4.zip \
  --info 'INFO.gnomad_popmax_af < 0.001 && variant.FILTER == "PASS"' \
  --sample-expr 'het:sample.het && sample.GQ > 20 && sample.DP > 12' \
  -o rare.vcf.gz

That drops you from 4.5 million variants to a few thousand rare, well-supported ones. The population frequency filter does most of the work: a variant at 5% frequency in gnomAD is not the explanation for anything rare.

9. Look at the defensible subsets first

Three categories in your genome have interpretation good enough to be worth reading, and everything else is weaker than it looks.

Known pathogenic variants in clinically actionable genes. Intersect your rare set with the ACMG secondary-findings gene list (v3.2, roughly 80 genes) and keep only ClinVar Pathogenic or Likely Pathogenic entries with two-star review status or better. This is the category where clinical WGS pipelines find real, actionable results, and it is also where a false positive does the most damage ⁴. Any variant in this category is a finding to take to a clinical geneticist for confirmation on a CLIA-validated assay. Such a finding is not a conclusion. Research-grade WGS is not a diagnosis. Diagnostic use of WGS is governed by validated pipelines, curated variant interpretation, and clinical context, none of which a laptop pipeline reproduces ⁵.

Carrier status for recessive conditions. Same filter, but you are looking at heterozygous P/LP variants in autosomal recessive genes. Interpretable, well-catalogued, relevant to reproductive planning. Discuss with a genetic counselor rather than acting on the VCF.

Pharmacogenomic star alleles. Run PharmCAT on your VCF, plus Cyrius for CYP2D6, which PharmCAT will not resolve from short reads alone. You get diplotypes mapped to CPIC phenotype categories for these genes and others.

• CYP2C19
• CYP2C9
• DPYD
• TPMT
• NUDT15
• SLCO1B1

These are laboratory measurements of your genotype. Any change to a medication is a conversation with the prescribing clinician, who will want a clinical-grade confirmation.

Polygenic scores sit outside those three. You can compute them (plink2 --score, PGS Catalog weights), and the arithmetic is trivial, but the ceiling on interpretation is low: scores transfer poorly across ancestries, and the significance thresholds for whole-genome association work are stricter than the familiar 5×10⁻⁸ because sequencing tests far more variants than an array does ⁶. Treat a PRS as a population-level statistic you happen to have computed on one person.

For the fun side of the data, psmc reconstructs effective population size history over hundreds of thousands of years from a single diploid genome, using only the density of heterozygous sites along each chromosome ⁷. It takes a few hours and tells you nothing about your health, and the method is a nice demonstration of how much history is encoded in one person’s heterozygosity.

Common problems

Chromosome naming mismatch. Your CRAM says chr1, your annotation VCF says 1, and bcftools annotate silently returns zero matches. Check with bcftools view -h file.vcf.gz | grep contig | head -3 on both files. Fix with bcftools annotate --rename-chrs.

Coverage looks fine, a specific gene does not. The following regions all have paralogs or high homology, so short reads there get low MAPQ and calls are unreliable or absent.

• PMS2
• SMN1/SMN2
• HBA1/HBA2
• CYP2D6
• STRC
• NEB
• the HLA region

Check the mosdepth per-base output over the gene before concluding anything about it. Use the paralog-aware caller (Cyrius, SMNCopyNumberCaller) or accept that short-read WGS does not answer that question.

Indels in homopolymers. A run of 8+ identical bases produces recurrent false indels. If a candidate variant sits in a homopolymer or a short tandem repeat, look at the reads in IGV before you spend an hour on it.

Duplicate rate above 20%. Usually a PCR-amplified library or over-sequencing of a low-complexity library. Marking duplicates removes them, so your effective coverage is lower than your raw read count suggested. Recompute coverage after markdup.

A VUS. Most rare variants you find will be Variants of Uncertain Significance, and most of them stay that way. A VUS is a statement about the state of the evidence. The label is not a statement about you. Do not convert one into a conclusion because it sits in a gene with an alarming name.

Storage. FASTQ + CRAM + VCFs + intermediate DeepVariant files run past 1 TB. Keep the CRAM, the reference FASTA, the gVCF, and the filtered VCF. Delete the FASTQs only after you have verified the CRAM reads cleanly with samtools quickcheck and you can regenerate FASTQ from it with samtools fastq.

Footnotes

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

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

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

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. Changjiang Xu, Ioanna Tachmazidou, Klaudia Walter, et al. Estimating Genome‐Wide Significance for Whole‐Genome Sequencing Studies. Genetic Epidemiology, 2014. https://doi.org/10.1002/gepi.21797 ↩

7. Heng Li, Richard Durbin. Inference of human population history from individual whole-genome sequences. Nature, 2011. https://doi.org/10.1038/nature10231 ↩