How to Analyze Your Own Genome: A Working Guide to Variant Analysis Tools

By the end of this guide you will have the following.

• an aligned CRAM of your own whole genome against GRCh38
• a small-variant VCF with roughly 4 to 5 million single-nucleotide variants and short insertions and deletions
• a measured estimate of how accurate that VCF is in the regions you care about
• an annotated version carrying gene consequences, population frequencies, and clinical assertions
• a reduced candidate list small enough to read

You will also have a separate structural-variant call set, a pharmacogenomic star-allele report, and polygenic scores computed under a defined reference panel. You need the raw reads as paired FASTQ or an unaligned CRAM, roughly 60 to 120 GB for a 30x genome. You also need about 500 GB of scratch disk, 32 GB of RAM and 16 or more cores. A GPU cuts DeepVariant runtime by an order of magnitude but is not required. Finally you need Docker or Apptainer and a conda or pixi environment. Everything below is measurement and interpretation. Nothing here is medical advice, and anything you would act on clinically needs confirmation in an accredited laboratory and a conversation with a clinician or genetic counselor.

1. Understand what you are about to measure

Variant analysis is the process of identifying the positions where your genome differs from a reference sequence, genotyping those differences, and attaching enough biological and population context to each one to decide which are worth a second look. The scale is set by human diversity: any two unrelated people are about 99.9% identical at the nucleotide level, which sounds like near-identity until you multiply the remaining fraction across three billion bases. Population-scale sequencing established the working numbers, with an individual genome carrying on the order of three to four million single-nucleotide variants, several hundred thousand short indels, and a few thousand larger structural changes.¹ Of those, only a few hundred are rare protein-truncating variants and a handful sit in genes with well-established clinical consequence.²

That arithmetic drives every design decision downstream. The pipeline is a funnel, and the interesting failure modes are at both ends: variants you never called because the caller could not see them, and variants you called correctly but drowned in four million neighbors.

2. Verify the inputs and pick a reference

Before compute, check integrity. Verify checksums against whatever your provider published, then confirm the files are not truncated:

samtools quickcheck -v *.cram && echo "cram ok"
seqkit stats -a -j 8 sample_R1.fastq.gz sample_R2.fastq.gz
fastqc -t 8 sample_R1.fastq.gz sample_R2.fastq.gz

Look at mean read length, duplicate estimate, and per-base quality tail. Also check total bases, where 30x of a 3.1 Gb genome is roughly 93 Gbp. If total bases divided by 3.1e9 is below about 25, expect measurably worse indel sensitivity and heterozygote recall.

For the reference, we use the GRCh38 analysis set with ALT contigs excluded from calling but present for decoy purposes, specifically GCA_000001405.15_GRCh38_no_alt_analysis_set.fna. The reason is practical rather than aesthetic: nearly every annotation resource you will want, including gnomAD, ClinVar coordinates, and the Polygenic Score Catalog harmonized files, is published on GRCh38. T2T-CHM13 is the better sequence and resolves regions GRCh38 misrepresents, but the annotation ecosystem around it is thinner, so we treat it as a second alignment for specific questions rather than the primary.

3. Align, mark duplicates, and store as CRAM

bwa-mem2 index ref.fna
bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:L1\tSM:SELF\tPL:ILLUMINA\tLB:lib1' \
  ref.fna sample_R1.fastq.gz sample_R2.fastq.gz \
| samtools fixmate -m -u -@ 4 - - \
| samtools sort -u -@ 8 -m 3G -T /scratch/st - \
| samtools markdup -@ 8 --reference ref.fna -O CRAM -o sample.cram -
samtools index sample.cram
samtools stats -@ 8 --reference ref.fna sample.cram > sample.stats

The -K 100000000 flag fixes the batch size so alignment is deterministic across thread counts, which matters when you want to reproduce a call set a year later. -Y uses soft clipping for supplementary alignments, which most structural-variant callers expect. CRAM with the reference cuts storage to roughly 15 to 20 GB per genome from 60 to 80 GB of BAM.

Check samtools stats for error rate (typically 0.2 to 0.5% for Illumina), insert size mode (350 to 500 bp for standard libraries), and duplicate fraction. Above roughly 15% duplicates, your effective depth is lower than the nominal figure and you should recompute coverage from non-duplicate bases.

4. Call small variants with two independent callers

This is the step where tool choice has the largest effect on your results, and the reason is documented. Independent pipelines applied to the same reads produce call sets that overlap substantially but not completely. A meaningful fraction of calls is unique to a single pipeline, concentrated in indels and low-complexity sequence.³ Later comparisons of GATK, FreeBayes, and SAMtools on whole genomes found the same pattern, with caller-specific sensitivity and precision tradeoffs rather than one uniformly correct answer.⁴ Recent exome benchmarking against gold-standard datasets reaches a similar conclusion and shows that the ranking depends on variant class and region.⁵

We run DeepVariant as the primary caller and GATK HaplotypeCaller as the secondary, then compare.

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

# Secondary: GATK4, single sample, hard filters
gatk --java-options "-Xmx16g" HaplotypeCaller \
  -R ref.fna -I sample.cram -O sample.hc.g.vcf.gz \
  -ERC GVCF --native-pair-hmm-threads 16
gatk GenotypeGVCFs -R ref.fna -V sample.hc.g.vcf.gz -O sample.hc.vcf.gz
gatk VariantFiltration -V sample.hc.vcf.gz -O sample.hc.filt.vcf.gz \
  --filter-name "QD2" --filter-expression "QD < 2.0" \
  --filter-name "FS60" --filter-expression "FS > 60.0" \
  --filter-name "MQ40" --filter-expression "MQ < 40.0" \
  --filter-name "SOR3" --filter-expression "SOR > 3.0"

Two notes on GATK. Variant Quality Score Recalibration needs a cohort, so on a single genome you use hard filters or a model-based single-sample filter, and both are weaker than cohort VQSR. Base quality score recalibration adds hours and changes very little on modern Illumina chemistry with 30x depth, so we skip it for single-sample work. If you want a reproducible, resource-conscious GATK4 workflow rather than assembling one yourself, OVarFlow packages the standard steps in a Snakemake pipeline with modest memory requirements.⁶

Then look at the disagreement explicitly:

bcftools isec -p isec_dv_hc -Oz sample.dv.vcf.gz sample.hc.filt.vcf.gz
# 0000.vcf.gz = DV only, 0001.vcf.gz = HC only, 0002/0003 = shared

The DV-only and HC-only sets are your uncertainty budget. If a variant you care about appears in only one, treat it as unconfirmed until you have looked at the reads in IGV or confirmed it by an independent assay.

5. Benchmark the pipeline before you trust it

Running a benchmark on a known sample is the single highest-value hour in this whole process, because it converts vague confidence into per-class numbers you can quote. Download the Genome in a Bottle HG002 reads and run them through exactly your pipeline. Then evaluate against the GIAB truth VCF and high-confidence BED with hap.py using the vcfeval engine, which handles representational differences in indels correctly:

docker run --rm -v "$PWD":/data pkrusche/hap.py /opt/hap.py/bin/hap.py \
  /data/HG002_GRCh38_1_22_v4.2.1_benchmark.vcf.gz \
  /data/hg002.dv.vcf.gz \
  -f /data/HG002_GRCh38_1_22_v4.2.1_benchmark_noinconsistent.bed \
  -r /data/ref.fna --engine=vcfeval --threads 16 \
  --stratification /data/stratification.tsv \
  -o /data/happy_hg002

Read the output by class and by stratum. Do not read it as a single headline number. Modern callers exceed 99.5% F1 for SNVs inside high-confidence regions, and indels are several tenths of a percent worse. Both degrade sharply in segmental duplications, homopolymers, and tandem repeats. Framing discovery as an optimization problem with explicit sensitivity and precision targets tuned against known truth is the approach we endorse, and it has been formalized for personal genomes.⁷ Save the stratified table next to your VCF. When a result later looks surprising, the first question is whether the variant sits in a region where your pipeline was measured to be unreliable.

6. Normalize and merge before annotating

Annotation tools assume left-aligned, biallelic records. Skipping normalization is the most common source of silently wrong annotations, because a variant written differently from the way ClinVar writes it will simply fail to match.

bcftools norm -f ref.fna -m -any -c w -Ou sample.dv.vcf.gz \
| bcftools norm -d exact -Ou \
| bcftools view -f PASS -i 'FMT/GQ>=20 && FMT/DP>=10' -Oz -o sample.norm.vcf.gz
bcftools index -t sample.norm.vcf.gz
bcftools stats sample.norm.vcf.gz | head -40

Sanity-check the totals. A 30x genome should give roughly 4.0 to 5.0 million variants after this filter, with an SNV-to-indel ratio near 8:1 and a transition-to-transversion ratio near 2.0 to 2.1 genome-wide. A Ti/Tv well below 2.0 usually means false positives are leaking through.

7. Annotate with VEP, gnomAD, ClinVar, and one or two predictors

We prefer Ensembl VEP run offline against a local cache. It is scriptable and versioned, and it gives you per-transcript consequences rather than a single collapsed verdict.

vep -i sample.norm.vcf.gz -o sample.vep.vcf.gz --vcf --compress_output bgzip \
  --cache --offline --dir_cache "$VEP_CACHE" --assembly GRCh38 \
  --fasta ref.fna --species homo_sapiens --fork 16 \
  --everything --pick_allele_gene --transcript_version \
  --custom file=gnomad.genomes.v4.1.sites.vcf.gz,short_name=gnomADg,format=vcf,type=exact,fields=AF%AF_popmax%nhomalt \
  --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.indel.tsv.gz \
  --plugin AlphaMissense,file=AlphaMissense_hg38.tsv.gz \
  --plugin SpliceAI,snv=spliceai_scores.raw.snv.hg38.vcf.gz,indel=spliceai_scores.raw.indel.hg38.vcf.gz

Be deliberate about predictors. Computational pathogenicity scores are correlated with each other, trained on overlapping data, and calibrated for missense variants in reasonably well-studied genes. Stacking six of them does not give you six independent opinions. We carry CADD as a general deleteriousness prior, AlphaMissense for missense, and SpliceAI for anything near a splice junction. All three are ranking signals rather than verdicts. For clinical-grade interpretation the field uses structured evidence frameworks in which computational prediction is one weighted line of evidence among many, and knowledge-driven platforms exist precisely because the annotation step alone does not resolve a case.⁸

8. Reduce four million variants to a readable list

Filtering is where the funnel closes. slivar is the tool we use, because its expressions are readable JavaScript over the annotated fields and it runs over a whole genome in minutes.

slivar expr --vcf sample.vep.vcf.gz --pass-only \
  -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' \
  --sample-expr 'hom:sample.hom_alt && sample.GQ > 20 && sample.DP > 12' \
  -o sample.rare.vcf.gz

Then split into two tracks. The first is a frequency-and-consequence track: rare (population maximum allele frequency under 0.1%), predicted loss-of-function or high-impact missense, restricted to a gene list you chose in advance. The second is an assertion track: every position with a ClinVar classification of pathogenic or likely pathogenic and a review status of at least two submitters, regardless of predicted consequence. The second list is short, usually a few dozen entries after you discard low-review-status assertions, and it is where the findings with real interpretive weight live.

Export to a table you can read or hand to an agent:

bcftools +split-vep sample.rare.vcf.gz -f '%CHROM\t%POS\t%REF\t%ALT\t%SYMBOL\t%Consequence\t%HGVSc\t%HGVSp\t%gnomADg_AF\t%ClinVar_CLNSIG\t%CADD_PHRED\t%am_pathogenicity\t[%GT\t%DP\t%GQ]\n' \
  -d -A tab > sample.candidates.tsv

Anything on the assertion track that appears clinically consequential should be confirmed on a clinically validated platform and reviewed with a clinician. Research-grade sequencing has a nonzero false-positive rate even at high quality, and penetrance for most variants in an unselected individual is lower than the case-ascertained literature implies.

9. Add the variant classes short reads miss

Small-variant calling ignores most structural variation. Short reads are blind or nearly blind to large insertions, many repeat expansions, and complex rearrangements. Run at least one structural-variant caller and one repeat caller:

manta configManta.py --bam sample.cram --referenceFasta ref.fna --runDir manta_out
manta_out/runWorkflow.py -j 16
ExpansionHunter --reads sample.cram --reference ref.fna \
  --variant-catalog variant_catalog_hg38.json --output-prefix sample.eh

Treat short-read structural calls as hypotheses. Precision for deletions is decent, for insertions and inversions much worse. If structural variation or repeat content matters for your question, long reads are the correct instrument, and the tooling around long-read genomes (alignment with minimap2, small variants with DeepVariant or Clair3, structural variants with Sniffles or cuteSV, plus direct methylation calling from the raw signal) is now mature enough to run end to end.⁹ Long reads also resolve phasing across most genes, which short reads do not.

10. Layer on pharmacogenomics and polygenic scores

Two analyses use different tooling than everything above and are worth running because they are well-defined and reproducible. Pharmacogenomic star-allele calling needs haplotype-aware logic over specific loci rather than per-variant annotation, which is why PharmCAT and pypgx exist:

pharmcat_pipeline sample.norm.vcf.gz -o pharmcat_out
pypgx run-ngs-pipeline CYP2D6 cyp2d6-pipeline \
  --variants sample.norm.vcf.gz --depth-of-coverage doc.zip --control-statistics cs.zip

Read the PharmCAT report as a genotype-to-phenotype summary with the evidence level attached to each guideline. CYP2D6 in particular involves copy number and hybrid alleles that short reads call poorly, so an indeterminate call there is common and appropriate. Any prescribing decision belongs with a clinician and pharmacist. The report does not make that decision.

For polygenic scores, use the Polygenic Score Catalog’s pipeline so the scoring, allele harmonization, and ancestry adjustment are all explicit:

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

Interpret the output as a percentile within a defined reference population and nothing more. Scores derived predominantly from European-ancestry cohorts transfer poorly to other ancestries, and a percentile is a distributional statement about a population rather than a prediction about you.

Common problems

Annotation mismatches from unnormalized input are the most frequent silent failure. If a variant you expect to carry a ClinVar assertion shows none, re-check that you ran bcftools norm with the same reference FASTA used for calling and that multiallelics were split.

Reference build confusion is the second. gnomAD v4 on GRCh38, a CADD file built for GRCh37, and coordinates copied from a paper without stating the build will produce annotations that look plausible and are wrong. Record the build and the version of every resource in the VCF header.

Low-complexity regions generate most of your false positives. Homopolymers, short tandem repeats, and segmental duplications account for a disproportionate share of caller disagreement, which is exactly what pipeline comparison studies find.³ Keeping a low-complexity BED alongside your call set and flagging overlaps prevents hours of chasing artifacts.

Coverage dropouts are easy to miss because a VCF shows no variant where there were no reads. Generate a callable-regions BED with mosdepth --quantize 0:1:10:30: and check that the gene you are interested in is fully covered before concluding it carries no variants.

Somatic thinking applied to germline data is a conceptual error worth naming. Somatic variant calling assumes low variant allele fractions, tumor-normal pairing, and different quality management. It also uses a different tool family with different filters.¹⁰ A germline pipeline run at 30x will not find low-frequency clonal variation, and it is not designed to.

Finally, resist over-interpreting single variants. The reason population sequencing reshaped human genetics is that effect sizes and penetrance only become estimable at scale, and most of what an individual genome contains is unremarkable by design.²¹ The value of your own call set is in the small number of well-supported findings and in the longitudinal comparisons you can make against it later.

Footnotes

1. The 1000 Genomes Project Consortium. A map of human genome variation from population-scale sequencing. Nature, 2010. https://doi.org/10.1038/nature09534 ↩ ↩²

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. Jason O’Rawe, Tao Jiang, Guangqing Sun, et al. Low concordance of multiple variant-calling pipelines: practical implications for exome and genome sequencing. Genome Medicine, 2013. https://doi.org/10.1186/gm432 ↩ ↩²

4. Anna Supernat, Oskar Valdimar Vidarsson, Vidar M. Steen, et al. Comparison of three variant callers for human whole genome sequencing. Scientific Reports, 2018. https://doi.org/10.1038/s41598-018-36177-7 ↩

5. Matthew Wong, Bryan Liew, Melissa Hum, et al. Benchmarking of variant calling software for whole-exome sequencing using gold standard datasets. Scientific Reports, 2025. https://doi.org/10.1038/s41598-025-97047-7 ↩

6. Jochen Bathke, Gesine Lühken. OVarFlow: a resource optimized GATK 4 based Open source Variant calling workFlow. BMC Bioinformatics, 2021. https://doi.org/10.1186/s12859-021-04317-y ↩

7. Gareth Highnam, Jason J. Wang, Dean Kusler, et al. An analytical framework for optimizing variant discovery from personal genomes. Nature Communications, 2015. https://doi.org/10.1038/ncomms7275 ↩

8. Dvir Dahary, Yaron Golan, Yaron Mazor, et al. Genome analysis and knowledge-driven variant interpretation with TGex. BMC Medical Genomics, 2019. https://doi.org/10.1186/s12920-019-0647-8 ↩

9. Bioinformatics for human long-read whole genome sequencing. Nature Reviews Methods Primers, 2026. https://doi.org/10.1038/s43586-026-00528-w ↩

10. Xiaoyu He, Shanyu Chen, Ruilin Li, et al. Comprehensive fundamental somatic variant calling and quality management strategies for human cancer genomes. Briefings in Bioinformatics, 2020. https://doi.org/10.1093/bib/bbaa083 ↩