How to Analyze Your Own Whole Genome Data

At the end of this you will have: a CRAM aligned to GRCh38 with per-base coverage statistics, a small-variant VCF with roughly 4–5 million variants called at 99%+ F1 against GIAB truth sets, a structural variant call set, a repeat-expansion genotype table, a pharmacogenomic diplotype report, and a handful of TSV files small enough to hand to an analysis agent without shredding its context window. You need the raw data (FASTQ or CRAM, not a PDF report), about 1 TB of free disk for a 30x genome, 64 GB of RAM, 16–32 cores, and a weekend. Everything below runs on Linux with Docker or Singularity available. If you only have a consumer array (23andMe, AncestryDNA), none of this applies: that is ~700k genotyped sites, not a genome, and the failure modes are completely different.

One boundary up front, and it holds for the whole piece. This is measurement and interpretation of your own data. It is not a diagnosis, and no output here should change a medication, a dose, or a screening schedule on its own. Variants that look consequential need confirmation in a CLIA/CAP-accredited clinical lab and interpretation with a genetic counselor or physician. Research-grade calls and clinical calls are different objects with different error profiles, and the gap between them is exactly where people hurt themselves.

1. Get the actual files, and check what you were given

Ask your provider for FASTQ (raw reads) or CRAM (aligned, reference-compressed) plus the MD5 checksums. A 30x human genome is roughly 100–120 GB as gzipped paired FASTQ, or 15–20 GB as CRAM against GRCh38. If someone offers you only a VCF, you can still do steps 7–10, but you lose the ability to re-call variants, inspect coverage, or look at a read pileup when a result looks wrong. That ability is the main reason to own your own data.

Two practical notes. Clinical WGS ordered for a diagnostic indication is sometimes covered by insurance, particularly for suspected rare disease in children, where it has largely displaced sequential single-gene and panel testing ¹. Elective sequencing in a healthy adult is generally out of pocket. And “30x” is a mean. The figure is not a floor. Verify it yourself in step 4.

md5sum -c checksums.md5
seqkit stats -a *.fastq.gz          # read count, mean length, Q20/Q30 fraction

Expect ~800M–1B read pairs at 2x150 for 30x. Q30 above 90% is normal for a healthy Illumina run. If the Q30 fraction is in the 70s, the run had problems and your downstream variant quality will show it.

2. Skip the trimming, run the QC report

Modern aligners soft-clip adapter sequence, and aggressive trimming costs you more than it buys on WGS. We run fastp with filtering disabled, purely to produce the report.

fastp -i R1.fastq.gz -I R2.fastq.gz --disable_adapter_trimming \
      -Q -L -A --thread 8 --html qc.html --json qc.json

Look at three things in the JSON: duplication rate (above ~20% on a PCR-free library means the library was low-complexity and your effective coverage is lower than advertised), insert size distribution (a mode around 350–450 bp is typical), and per-cycle quality decay.

3. Align to GRCh38, no-alt analysis set

Use GCA_000001405.15_GRCh38_no_alt_analysis_set.fna. The no-alt version keeps the decoy and HLA contigs but drops the ALT haplotypes, which avoids a class of mapping-quality-zero artifacts that silently delete variants in immune and other polymorphic regions. Do not align to GRCh37 unless you have a specific reason: almost every current annotation resource ships GRCh38 first, and liftover between builds is a reliable source of wrong coordinates.

bwa-mem2 index GRCh38_no_alt.fna

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

-K 100000000 fixes the batch size so the run is deterministic. -Y uses soft clipping for supplementary alignments, which structural variant callers want. The read group is not optional: GATK and most downstream tools will refuse to run without SM.

Mark duplicates and convert to CRAM:

samtools collate -@ 8 -O -u S1.sorted.bam tmp_collate \
| samtools fixmate -@ 8 -m -u - - \
| samtools sort -@ 8 -m 2G -u - \
| samtools markdup -@ 8 -f markdup.stats --reference GRCh38_no_alt.fna \
    -O CRAM,version=3.0 - S1.cram
samtools index S1.cram

On 32 cores this is 6–10 hours. CRAM against the same reference is about 40% the size of the equivalent BAM. Keep the reference FASTA forever: CRAM is unreadable without the exact reference it was compressed against, and this is the single most common way people lose access to their own data.

4. QC the alignment before you believe any variant

This step is where most amateur analyses go wrong. Four checks, in order.

Coverage, with mosdepth:

mosdepth -n --fast-mode --by 1000 -t 4 S1 S1.cram
mosdepth --quantize 0:1:5:15:150: -t 4 S1.q S1.cram

The quantize output gives you a BED of callable territory. On a good 30x genome, roughly 95% of the autosomes should be at ≥15x. If your ≥15x fraction is 85%, your sensitivity for heterozygous variants is meaningfully degraded and you should say so when interpreting any negative result.

Contamination, with VerifyBamID2:

VerifyBamID --SVDPrefix 1000g.phase3.100k.b38.vcf.gz.dat \
            --Reference GRCh38_no_alt.fna --BamFile S1.cram

FREEMIX under 0.01 is clean. Above 0.03, you will see spurious low-allele-fraction heterozygotes and you should stop and ask for a resequence.

Sex and relatedness, with somalier:

somalier extract -d extracted/ --sites sites.hg38.vcf.gz -f GRCh38_no_alt.fna S1.cram
somalier relate extracted/*.somalier

Check that inferred sex chromosome copy number matches expectation. A mismatch usually means a sample swap, and a sample swap invalidates everything downstream.

Finally, samtools stats S1.cram | grep ^SN for error rate and mapped fraction. Mapped should be above 99%, error rate around 0.2–0.5%.

5. Call small variants with DeepVariant

For a single sample on Illumina short reads, DeepVariant is our default over GATK HaplotypeCaller. It needs no base quality score recalibration step, no hard-filter tuning, and it produces better indel accuracy out of the box.

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

Three to five hours on 32 cores. Keep the gVCF: it distinguishes “reference at this position with 40x support” from “no data here,” which is the difference between a real negative and an absence of evidence. Expect 4.4–5.0M variants, of which ~4M are SNVs, with a transition/transversion ratio near 2.0 genome-wide. A Ti/Tv far from 2.0 means something is wrong with the call set.

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

Normalize before doing anything else. Left-align indels and split multi-allelics, or your joins against ClinVar and gnomAD will silently miss.

bcftools norm -f GRCh38_no_alt.fna -m -any -Oz -o S1.norm.vcf.gz S1.vcf.gz
bcftools index -t S1.norm.vcf.gz

6. Structural variants and repeat expansions

Short-read SV calling is the weakest part of this pipeline, and one of the standing arguments for long reads ². Run it anyway, but treat every call as a hypothesis to be confirmed visually in IGV.

configManta.py --bam S1.cram --referenceFasta GRCh38_no_alt.fna --runDir manta_S1
manta_S1/runWorkflow.py -j 16

Manta gives deletions, duplications, insertions, inversions, and translocations from split reads and discordant pairs. Expect a few thousand calls, dominated by common polymorphic deletions. Filter to PASS and intersect against a population SV catalog to drop the ones everyone carries.

Repeat expansions are invisible to standard variant callers and are worth genotyping explicitly:

ExpansionHunter --reads S1.cram --reference GRCh38_no_alt.fna \
  --variant-catalog variant_catalog/hg38/variant_catalog.json \
  --output-prefix S1_eh

Short reads cannot size expansions longer than the fragment length with any precision. ExpansionHunter is a screen, not a measurement, and any allele near or above a pathogenic threshold requires a clinical lab assay (repeat-primed PCR or Southern blot) ordered by a physician.

7. Annotate

Ensembl VEP with cache, offline, plus custom annotations for gnomAD and ClinVar:

vep -i S1.norm.vcf.gz -o S1.vep.vcf.gz --vcf --compress_output bgzip \
  --cache --offline --dir_cache "$HOME/.vep" --assembly GRCh38 \
  --fasta GRCh38_no_alt.fna --everything --pick_allele_gene --fork 8 \
  --custom file=gnomad.genomes.v4.1.sites.vcf.gz,short_name=gnomADg,format=vcf,type=exact,fields=AF%AF_grpmax \
  --custom file=clinvar.vcf.gz,short_name=ClinVar,format=vcf,type=exact,fields=CLNSIG%CLNREVSTAT%CLNDN

This takes an hour or two and produces a large VCF. --pick_allele_gene collapses transcript consequences to one per gene per allele, which is what you want for a first pass, though you lose transcript-specific detail you may need later for a specific gene.

8. Reduce to something a human can read

Five million variants is not an analysis. The reduction that matters is: high-confidence pathogenic assertions in ClinVar, plus rare protein-truncating variants in genes with established disease associations.

bcftools +split-vep S1.vep.vcf.gz -d -A tab \
  -f '%CHROM\t%POS\t%REF\t%ALT\t%SYMBOL\t%Consequence\t%ClinVar_CLNSIG\t%ClinVar_CLNREVSTAT\t%gnomADg_AF\t[%GT\t%GQ\t%DP]\n' \
  -i 'ClinVar_CLNSIG~"athogenic" && ClinVar_CLNREVSTAT~"criteria_provided" && FMT/GQ>20 && FMT/DP>10' \
  > clinvar_hits.tsv

You will get somewhere between 20 and 80 rows. Nearly all will be heterozygous carrier states for autosomal recessive conditions, which carry reproductive relevance and essentially no personal health relevance. Restricting to the ACMG secondary-findings gene list gives a much smaller, more actionable set: about 2% of unselected individuals carry a reportable variant in those genes. Interpretation of any hit belongs with a clinician. Variant classification shifts over time, laboratories disagree, and penetrance in an unselected adult is usually far lower than the penetrance reported in ascertained families ³. A “pathogenic” label derived from a research pipeline is a reason to consult. The label is not a conclusion.

Where the pipeline itself is the source of error, the fix is validation against known-truth samples. Public health and clinical labs formalize this by running reference material through the identical pipeline and measuring concordance before trusting any patient result ⁴. Download the GIAB HG002 data, run your exact commands, and compare with hap.py against the v4.2.1 truth set. If you have not done that, you do not know your own false positive rate.

9. Pharmacogenomics and polygenic scores

PharmCAT takes your VCF and produces star-allele diplotypes and guideline-linked annotations. It covers these genes and others.

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

pharmcat_pipeline S1.norm.vcf.gz -o pharmcat_out

CYP2D6 needs a dedicated caller (Cyrius or StellarPGx) because of the pseudogene and hybrid alleles. Do not change any medication or dose based on this output. Bring the report to your prescriber or pharmacist.

Polygenic scores with PLINK2:

plink2 --vcf S1.norm.vcf.gz --score weights.txt 1 2 4 header cols=+scoresums --out prs_out

A raw score for one person is a number with no scale. To place yourself, score the 1000 Genomes samples with the same weight file and compute your percentile within the ancestry group that matches the GWAS training population. Scores trained in European-ancestry cohorts lose a large fraction of their predictive accuracy when applied to individuals of African ancestry. That attenuation is a property of the score. The attenuation is not a property of the person ².

10. Make it queryable

Keep the VCF and CRAM as source of truth in cold storage. For daily work, flatten the annotated variants into Parquet and query with DuckDB.

bcftools +split-vep S1.vep.vcf.gz -d -A tab -f '...' | \
  duckdb -c "CREATE TABLE v AS SELECT * FROM read_csv_auto('/dev/stdin'); COPY v TO 'variants.parquet';"

If you plan to hand this to an AI agent, give it the Parquet file and a schema. Do not give it the VCF. Pasting raw variant lines into a chat window burns context and produces confident nonsense about positions the model never read. Give the agent a query tool and a written description of your coverage gaps.

Common problems

Coverage reported as a mean hides the gaps. GC-rich first exons, segmental duplications, and the MHC routinely drop below 10x even in a 30x genome. Always report a negative result alongside the callable fraction of the gene in question.

Build mixing. A GRCh37 coordinate pasted into a GRCh38 analysis lands in the wrong gene roughly as often as it lands in the right one. Put the build in every filename.

Stale ClinVar. The database changes weekly. A variant classified as pathogenic in 2019 may be a variant of uncertain significance now, or reclassified as benign. Re-annotate before acting on anything, and record the ClinVar release date with your results.

Multi-allelic sites. If you skipped bcftools norm -m -any, your ClinVar join misses every site where you carry one of two alternate alleles. This failure is silent.

Treating research calls as clinical results. Diagnostic-grade WGS involves validated pipelines, defined confirmation policies, and interpretation by trained staff, and the interpretive bottleneck is larger than the sequencing one ¹. Your pipeline reproduces the sequencing part and none of the rest.

Mitochondrial and Y calls. Heteroplasmy needs a dedicated caller (Mutect2 in mitochondria mode). Diploid callers will mis-handle chrM entirely.

Expecting an answer. The most likely outcome of all of this for a healthy adult is a few carrier states, a pharmacogenomic profile worth showing a pharmacist, and several hundred variants of uncertain significance. The value is a permanent, queryable baseline you can re-annotate as knowledge improves. The baseline is not a verdict.

Footnotes

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

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

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

4. Kelly F. Oakeson, Jennifer M. Wagner, Michelle Mendenhall, et al. Bioinformatic Analyses of Whole-Genome Sequence Data in a Public Health Laboratory. Emerging infectious diseases, 2017. https://doi.org/10.3201/eid2309.170416 ↩