From FASTQ to an Annotated Variant Table: Analyzing Your Own Genome

By the end of this you will have a CRAM aligned to GRCh38 and a normalized VCF of roughly 4–5 million small variants. You will have a structural variant and repeat-expansion call set. You will also have a flat annotated table (Parquet or TSV) keyed by variant. The table carries gnomAD frequencies, ClinVar assertions, predicted consequences, and per-site read depth. That table is the thing you query, and the thing an LLM agent can reason over without hallucinating coordinates. You need: your raw data (paired FASTQ, or a CRAM/BAM from the lab), a machine with 16+ cores and 64 GB RAM, about 500 GB of scratch disk, Docker, and a few days of wall clock if you are doing it end to end on one box. Everything below is measurement and interpretation. Nothing here is a diagnosis, and any variant you plan to act on needs confirmation in a CLIA/CAP laboratory and a conversation with a clinical geneticist or genetic counselor.

1. Inventory what you were given

Labs hand back very different things under the label “your genome.” Check before you build anything.

md5sum -c checksums.md5
seqkit stats -a sample_R1.fastq.gz sample_R2.fastq.gz
samtools quickcheck -v sample.cram && echo OK
samtools view -H sample.cram | grep -E '^@SQ' | head
samtools view -H sample.cram | grep -E '^@RG|^@PG'

Things to establish in the first ten minutes:

• Read length and layout. 2×150 bp PCR-free is the standard for a 30x human genome. 2×100 or PCR-positive libraries change your duplicate rate and indel error profile.
• Reference build. @SQ lines with SN:chr1 LN:248956422 mean GRCh38. LN:249250621 means GRCh37/hg19. If the CRAM has no embedded reference, you need the exact FASTA it was aligned against or samtools cannot decode it.
• Coverage. A 30x PCR-free genome is roughly 100–120 GB of gzipped FASTQ, 60–90 GB as BAM, 15–20 GB as CRAM. If your file is 3 GB, you have an exome or a genotyping array export. You do not have a genome.
• Whether you got a gVCF. If yes, keep it, but call your own variants anyway. Caller version and filtering thresholds matter more than most people expect, and you want to know which ones were used.

If what you have is a consumer array export (a ~25 MB text file of ~600k genotypes), most of this guide does not apply. Arrays assay preselected common sites and cannot see the rare coding variants that carry most of the interpretable signal.

2. Pick a reference and align

Use the GRCh38 analysis set with alt contigs and decoys (GCA_000001405.15_GRCh38_no_alt_analysis_set if your downstream tools are alt-unaware, the full analysis set if you will run an alt-aware post-processing step). Do not build a custom reference. Every annotation resource you want to join against is keyed to a public build.

bwa-mem2 index GRCh38_full_analysis_set_plus_decoy_hla.fa

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:L1\tSM:SAMPLE\tLB:lib1\tPL:ILLUMINA' \
  GRCh38_full_analysis_set_plus_decoy_hla.fa \
  sample_R1.fastq.gz sample_R2.fastq.gz \
| samtools fixmate -m -u - - \
| samtools sort -u -@ 8 -T /scratch/tmp -o - \
| samtools markdup -@ 8 --reference GRCh38_full_analysis_set_plus_decoy_hla.fa \
    - sample.cram
samtools index sample.cram

-K 100000000 fixes the batch size so the output is deterministic across thread counts. -Y uses soft clipping for supplementary alignments, which SV callers like Manta want. Skip BQSR unless you are running GATK HaplotypeCaller. DeepVariant does not benefit from it.

Then measure what you got, because coverage is the single number that determines what the rest of the pipeline can and cannot see:

mosdepth --by 1000 --fast-mode -t 4 sample sample.cram
samtools stats sample.cram > sample.stats

Look at mean depth, the fraction of the callable genome at ≥10x and ≥20x, duplicate rate (under 10% for PCR-free), and insert size distribution. At 30x you should be above 95% of the genome at ≥10x. Low-coverage sequencing (1–5x) is a real and useful design for population-scale allele frequency work, but it does not give you reliable individual genotypes at rare sites, which is exactly what personal interpretation depends on ¹.

3. Call small variants

We use DeepVariant for germline SNVs and indels on Illumina data. It is a single container, it needs no BQSR, and its indel precision in homopolymers is better than a default HaplotypeCaller run.

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

Expect on the order of 4–5 million small variants against GRCh38, of which roughly 20% are indels. Runtime is a few hours on 32 cores. Keep the gVCF: it records confident-reference blocks, so later you can distinguish “reference at this position” from “no coverage here,” which matters when a report says you do not carry a variant.

Mitochondrial DNA needs its own path (Mutect2 in mitochondrial mode or a dedicated heteroplasmy caller) because copy number is high and heteroplasmy is a continuous quantity. It is not a genotype.

4. Normalize before you annotate

This step is boring and it is where most homemade pipelines silently break. An indel can be written several ways at the same locus, and joins against ClinVar or gnomAD fail if representations differ.

bcftools norm -m -any -f GRCh38.fa -Ou sample.dv.vcf.gz \
| bcftools norm -d exact -Oz -o sample.norm.vcf.gz
bcftools index -t sample.norm.vcf.gz
bcftools stats sample.norm.vcf.gz | grep -E '^SN|^TSTV'

Check the Ts/Tv ratio: about 2.0–2.1 genome-wide, around 3.0 in coding regions. A genome-wide Ts/Tv of 1.6 means you have a pile of false positive transversions and should look at your filtering before going further.

5. Annotate with VEP

Ensembl VEP with a local cache is the workhorse. Run it offline, because the REST API will not survive 4 million variants.

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=clinvar_20250xxx.vcf.gz,short_name=ClinVar,format=vcf,type=exact,fields=CLNSIG%CLNREVSTAT%CLNDN \
  --custom file=gnomad.genomes.v4.1.sites.vcf.gz,short_name=gnomADg,format=vcf,type=exact,fields=AF%AF_popmax \
  --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 \
  --plugin LOFTEE,loftee_path:$HOME/loftee,human_ancestor_fa:human_ancestor.fa.gz

--pick_allele_gene collapses VEP’s many transcript consequences to one per gene per allele, which keeps the output table tractable. Keep --mane so you know which transcript a clinical report would have used. gnomAD v4 gives you allele frequencies from 730,947 exomes and 76,215 genomes, which is the main tool for dismissing the thousands of rare-looking variants that are simply common in some population.

Most of your variants are non-coding, and their interpretation is far weaker than coding interpretation. Regulatory annotations from ENCODE give you cell-type-specific candidate enhancers and promoters to intersect against. An overlap with a DNase-hypersensitive site is a hypothesis and is not evidence of effect ². Cross-species conservation, the kind of signal first mapped at scale by the mouse genome comparison, remains one of the more useful filters for non-coding sites ³.

6. Reduce to a table you can read

bcftools +split-vep sample.vep.vcf.gz -d -f \
 '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t[%GT\t%DP\t%AD\t%GQ\t]%SYMBOL\t%Consequence\t%IMPACT\t%HGVSc\t%HGVSp\t%CANONICAL\t%MANE_SELECT\t%gnomADg_AF\t%CADD_PHRED\t%SpliceAI_pred_DS_AG\t%LoF\t%ClinVar_CLNSIG\t%ClinVar_CLNREVSTAT\n' \
 -A tab | bgzip > sample.table.tsv.gz

Load into DuckDB and work from SQL rather than shell pipelines:

CREATE TABLE v AS SELECT * FROM read_csv_auto('sample.table.tsv.gz', delim='\t');

SELECT SYMBOL, HGVSp, GT, DP, gnomADg_AF, ClinVar_CLNSIG, ClinVar_CLNREVSTAT
FROM v
WHERE ClinVar_CLNSIG LIKE '%Pathogenic%'
  AND ClinVar_CLNREVSTAT LIKE '%multiple_submitters%'
  AND (gnomADg_AF IS NULL OR gnomADg_AF < 0.01)
  AND DP >= 15 AND GQ >= 30;

That query typically returns a handful of rows in a genome, most of which are carrier-state recessive variants or low-penetrance risk alleles. The one-star and no-assertion ClinVar entries are where nearly all the false alarms live, which is why the review-status filter is not optional. Expect roughly 100–200 predicted loss-of-function variants per genome, the great majority of which are in genes with no known disease association or are rescued by the second allele.

The ACMG secondary findings list (v3.2, 81 genes) is a sensible starting scope for “things a clinical lab would report back regardless of why you sequenced.” Restrict to that gene set first, then widen. Diagnostic-grade WGS interpretation is a structured process run by trained analysts, and it finds real answers in cases where targeted testing did not ⁴⁵. Your version of it at home is a screen that generates questions for a clinician. It is not a report.

7. The regions short reads handle badly

Small-variant calling covers maybe 90% of what is interpretable and misses categories that matter.

Structural variants:

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

Repeat expansions (Huntington’s CAG, FMR1 CGG, C9orf72 GGGGCC, and about 60 others):

ExpansionHunter --reads sample.cram --reference GRCh38.fa \
  --variant-catalog variant_catalog.json --output-prefix sample.eh

Short reads cannot size an expansion longer than the fragment length with any confidence. ExpansionHunter tells you “large” and that is the limit. Sizing needs PCR or long reads in a clinical lab.

Pharmacogenes: run PharmCAT over your VCF, and Cyrius for CYP2D6 specifically, because that locus has structural variation and a pseudogene that generic callers get wrong. PharmCAT emits diplotypes and phenotype calls. Read them as measurements. Any change to a medication belongs to your prescriber or a pharmacist, and the genotype is one input among many.

HLA: HLA-LA gives you G-group typing at class I and II from the alignment.

Polygenic scores: use the PGS Catalog with pgsc_calc, which handles strand harmonization and ancestry-aware standardization.

nextflow run pgscatalog/pgsc_calc -profile docker \
  --input samplesheet.csv --target_build GRCh38 --pgs_id PGS000018 --run_ancestry ref.tar.zst

A raw PRS number is meaningless without the reference distribution and without matching ancestry. Scores derived in European-ancestry cohorts lose substantial predictive accuracy in other populations, and the percentile is the only interpretable output.

8. Freeze the artifacts and plan reanalysis

Keep the CRAM, the gVCF, the normalized VCF, the annotation database versions, and a manifest of tool versions and hashes. Then rerun annotation, not alignment, on a schedule. ClinVar reclassifies variants continuously, gnomAD grows, and gene-disease associations get added. Periodic reanalysis of unchanged sequence data produces new findings, which is a well-documented pattern in diagnostic programs ⁵. A yearly vep rerun against a fresh cache and a diff of the ClinVar significance column costs an afternoon.

One note on framing. The sequencing itself has been cheap for years, and the “thousand-dollar genome” number refers to reagents and instrument time, not to interpretation, storage, clinical follow-up, or the counseling required when something ambiguous turns up ⁶⁷. Most of the real cost of a genome is downstream of the machine ⁸. The same reasoning applies to what you find, and interpretation is the expensive part. The ethical weight of incidental findings about disease risk, carrier status, and relatives is real. It should be decided before you look and not afterward ⁹.

Common problems

Build and contig mismatches. chr1 versus 1, GRCh37 versus GRCh38, different contig ordering in the header. Symptoms: VEP annotating nothing, bcftools erroring on the reference, or annotations landing on the wrong gene. Fix with bcftools annotate --rename-chrs and check with bcftools view -h | grep contig. Liftover from GRCh37 fails on a few percent of sites and inverts strand on some, so prefer realigning from FASTQ.

Skipping normalization. Multi-allelic records and left-unaligned indels mean your ClinVar join misses real hits. Always bcftools norm -m -any -f ref first.

Reading absence as reassurance. A variant not in your VCF may be a variant in a region with 3x coverage. Before concluding you do not carry something, check the gVCF and mosdepth output for that interval. Genes with recent segmental duplications routinely have unreliable short-read calls. Examples are PMS2, SMN1/SMN2, CYP2D6 and GBA.

Trusting single-submitter ClinVar entries. Filter on CLNREVSTAT for two-star or better before you spend an evening on a gene.

Over-reading VUS. A variant of uncertain significance is a statement about missing evidence. CADD and SpliceAI scores are prioritization tools. They are not classifications. Somatic-cancer intuitions do not transfer either: driver-versus-passenger logic in tumor genomes rests on recurrence across large cohorts and matched normal tissue ¹⁰, and none of that framework applies to a germline variant in one person.

Sample identity errors. Confirm the X/Y coverage ratio matches your expectation and, if you have array data or a relative’s genome, check concordance with bcftools gtcheck. Sample swaps happen, and they are easier to catch at the start than after a week of analysis.

Disk and memory. VEP with --everything --fork 16 and CADD will use 30+ GB of RAM and produce an annotated VCF several times the size of the input. Write intermediates to fast local scratch. Do not write them to network storage.

Footnotes

1. Runyang Nicolas Lou, Arne Jacobs, Aryn Wilder, et al. A beginner’s guide to low-coverage whole genome sequencing for population genomics. 2020. https://doi.org/10.22541/au.160689616.68843086/v1 ↩

2. The ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human genome. Nature, 2012. https://doi.org/10.1038/nature11247 ↩

3. Mouse Genome Sequencing Consortium. Initial sequencing and comparative analysis of the mouse genome. Nature, 2002. https://doi.org/10.1038/nature01262 ↩

4. Sarra L. Ryan, John F. Peden, Zoya Kingsbury, et al. Whole genome sequencing provides comprehensive genetic testing in childhood B-cell acute lymphoblastic leukaemia. Leukemia, 2023. https://doi.org/10.1038/s41375-022-01806-8 ↩

5. David R. Murdock, Jill A. Rosenfeld, Brendan Lee. What Has the Undiagnosed Diseases Network Taught Us About the Clinical Applications of Genomic Testing?. Annual Review of Medicine, 2022. https://doi.org/10.1146/annurev-med-042120-014904 ↩ ↩²

6. Timothy Caulfield, Jim Evans, Amy McGuire, et al. Reflections on the Cost of “Low-Cost” Whole Genome Sequencing: Framing the Health Policy Debate. PLoS Biology, 2013. https://doi.org/10.1371/journal.pbio.1001699 ↩

7. Laura Bonetta. Whole-Genome Sequencing Breaks the Cost Barrier. Cell, 2010. https://doi.org/10.1016/j.cell.2010.05.034 ↩

8. Unknown. The ‘Thousand-Dollar Genome’: An Ethical Exploration. Jahrbuch für Wissenschaft und Ethik, 2012. https://doi.org/10.1515/jfwe.2012.273 ↩

9. Wim Pinxten, Heidi Carmen Howard. Ethical issues raised by whole genome sequencing. Best Practice & Research Clinical Gastroenterology, 2014. https://doi.org/10.1016/j.bpg.2014.02.004 ↩

10. The Cancer Genome Atlas Research Network. Comprehensive genomic characterization defines human glioblastoma genes and core pathways. Nature, 2008. https://doi.org/10.1038/nature07385 ↩