Genome Analysis Tools: What to Run on Your Own WGS Data

If you have your own whole-genome sequencing data and want to analyze it yourself, the short answer is this. Align with bwa-mem2 or minimap2. Call small variants with DeepVariant (or GATK HaplotypeCaller if you need gVCF joint-calling semantics). Call structural variants with Manta plus a read-depth caller. Normalize and annotate with bcftools norm and Ensembl VEP. Check quality with mosdepth, somalier, and verifybamid2 before you believe anything. Everything else sits downstream of that spine. That includes the browsers, the annotation databases and the polygenic score tools. This page covers each step, the flags we would use, and the places where a personal genome analysis quietly goes wrong.

Nothing here is diagnostic. A variant that looks interesting in a VCF is a hypothesis about your biology and nothing more. Clinical interpretation of a specific variant requires a clinical-grade laboratory and a genetics professional, and the data you generate at home does not substitute for either.

The file formats you will handle

Six formats carry almost everything:

• FASTQ (.fq.gz): raw reads plus per-base Phred qualities. A 30x human genome is roughly 100 GB of gzipped paired FASTQ.
• BAM/CRAM: aligned reads. CRAM is reference-compressed and typically 30-50% the size of the equivalent BAM, which matters when you are keeping this forever. Always store the exact reference FASTA next to the CRAM, because CRAM cannot be decoded without it.
• VCF/BCF: variant calls. BCF is the binary form and is much faster for bcftools operations.
• gVCF: per-base confidence records including reference blocks, used when calls will later be merged across samples.
• BED: interval lists, used for confident regions, target panels, and masking.
• GFF3/GTF: gene models, used by annotators and browsers.

Sequencing costs fell fast enough that the analysis, not the data generation, became the bottleneck, and the tool ecosystem grew accordingly ¹. A survey of variant analysis tools a decade ago already counted more than two hundred distinct packages across alignment, calling, and annotation ². Most of them you can ignore.

Reference genome: pick one and be rigid about it

Use GRCh38 with the analysis set that includes alt contigs, decoys, and HLA sequences (the GCA_000001405.15_GRCh38_no_alt_analysis_set variant if your caller is not alt-aware). Do not mix references. The single most common failure in personal genome projects is a coordinate mismatch: a VCF called against GRCh37, annotated against a GRCh38 gene model, and interpreted against a ClinVar release keyed to a third build. Chromosome naming (chr1 versus 1) will bite you in the same way, silently producing zero overlaps in a bedtools intersect.

T2T-CHM13 is a better reference in the sense of completeness, and it resolves segmental duplications and false duplications that corrupt calls in GRCh38 regions. We still call on GRCh38 as the primary, because nearly all annotation resources (gnomAD, ClinVar and PGS Catalog scoring files) are keyed to it. We lift to CHM13 only when investigating a region known to be broken in GRCh38.

Alignment

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:S1\tSM:S1\tLB:L1\tPL:ILLUMINA' \
  GRCh38.fa r1.fq.gz r2.fq.gz \
| samtools sort -@ 8 -m 2G -o S1.bam -
samtools index S1.bam

-K 100000000 fixes the chunk size so output is deterministic regardless of thread count. -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 or will mis-handle the sample without it.

For duplicate marking use samtools markdup (after samtools fixmate -m) or gatk MarkDuplicates. On a PCR-free library, duplicate rates above roughly 5-10% suggest a library complexity problem worth raising with whoever prepared it. Skip GATK base quality score recalibration for a modern PCR-free Illumina library called with DeepVariant. DeepVariant’s model was trained on uncalibrated reads and BQSR provides no measurable benefit while costing hours.

Small variant calling

Two reasonable choices.

DeepVariant, run through its Docker or Singularity image:

run_deepvariant \
  --model_type=WGS \
  --ref=GRCh38.fa \
  --reads=S1.bam \
  --output_vcf=S1.vcf.gz \
  --output_gvcf=S1.g.vcf.gz \
  --num_shards=32

GATK HaplotypeCaller, if you want the classic workflow:

gatk HaplotypeCaller -R GRCh38.fa -I S1.bam \
  -O S1.g.vcf.gz -ERC GVCF --native-pair-hmm-threads 8

GATK’s local de novo reassembly of active regions is what made it accurate on indels, and the MapReduce-style sharding is what made it tractable on whole genomes ³. For a single sample, we prefer DeepVariant. It is more accurate on indels and in homopolymers, and it needs no cohort for filtering (VQSR degrades badly on n=1). The hard-filter fallback for a single-sample GATK VCF is crude. Expect roughly 4-5 million variants against GRCh38, of which about 3.5 million are SNVs and several hundred thousand are indels. If your count is far outside that range, suspect contamination, low coverage, or a reference mismatch before you suspect biology.

Structural variants, repeats, and the things short reads get wrong

Small variants are the easy part. A large fraction of functional variation is structural, and short-read WGS recovers only some of it ⁴.

• manta for split-read and discordant-pair SVs (deletions, duplications, inversions and some insertions).
• smoove (Lumpy plus filtering) or delly as a second caller. Take the intersection when you care about precision, the union when you care about sensitivity.
• cnvpytor or GATK gCNV for read-depth copy number, which catches large events the paired-end callers miss.
• ExpansionHunter for known repeat expansions, run with the Illumina variant catalog. Short reads cannot size expansions longer than the fragment length reliably, so treat large calls as flags for orthogonal testing. Do not treat them as answers.

Long reads (PacBio HiFi or ONT) change this picture entirely and resolve SVs, phasing, and segmental duplications that short reads cannot ⁵. If SVs are your interest and budget allows, sequence long.

Normalization before annotation

This step is skipped constantly and it corrupts every lookup that follows. The same indel can be written several ways, and a database join on an unnormalized representation silently misses.

bcftools norm -m -any -f GRCh38.fa -c x S1.vcf.gz \
| bcftools norm -d exact -Ob -o S1.norm.bcf
bcftools index S1.norm.bcf

The -m -any flag splits multiallelic records into biallelic rows. The -f flag left-aligns and trims against the reference. The -c x flag errors out on REF mismatches, which is how you discover you used the wrong build.

Annotation

Ensembl VEP is the one we run, with a local cache so it is fast and reproducible:

vep -i S1.norm.bcf -o S1.vep.tsv --tab --cache --offline \
  --assembly GRCh38 --everything --pick_allele_gene \
  --plugin CADD,whole_genome_SNVs.tsv.gz \
  --custom gnomad.genomes.v4.sites.vcf.gz,gnomADg,vcf,exact,0,AF

--pick_allele_gene gives one consequence per allele per gene, which keeps the table joinable. snpEff is a lighter alternative with a simpler install. ANNOVAR is still widely used but its licensing and database staleness make it a poor default.

Then filter. A useful first pass on a personal genome: keep variants with gnomAD allele frequency below 0.001, a high- or moderate-impact consequence, and genotype quality above 20. That usually leaves a few hundred to a few thousand rows. Almost all of them mean nothing. Interpretation frameworks exist for exactly this reason, and a rare protein-truncating variant in a gene with no established disease association is not evidence of anything ⁴. If a variant lands in ClinVar as pathogenic in a gene with real clinical consequence, stop analyzing and take it to a genetic counselor for confirmatory testing in a certified lab.

QC you should run before any of the above matters

• mosdepth --fast-mode --by 1000 S1 S1.bam: coverage distribution. For a 30x nominal genome, check that the mean callable coverage is near 30 and that the fraction of the genome above 10x is above roughly 95%.
• verifybamid2: contamination estimate against 1000 Genomes panel allele frequencies. Anything above about 1-2% contamination distorts heterozygous calls and inflates variant counts.
• somalier extract then somalier relate: confirms sex, ancestry projection, and that the sample is who you think it is. Sample swaps happen.
• bcftools stats plus plot-vcfstats: check the Ti/Tv ratio (expect around 2.0-2.1 genome-wide) and the heterozygous/homozygous ratio, which is ancestry-dependent.
• multiqc . to collect the above into one report.

Whole-genome analysis is a sequence of standard, composable steps, and the value comes from running them carefully rather than from any single clever tool ⁶.

Questions people also ask

What tools are used to analyze DNA? For short-read human WGS, the stack is as follows.

• bwa-mem2 or minimap2 for alignment
• samtools/bcftools for file manipulation and normalization
• DeepVariant or GATK for small variants
• Manta/smoove/cnvpytor for structural variants
• VEP or snpEff for annotation
• IGV for looking at reads when a call seems implausible
• snakemake or nextflow (with the nf-core/sarek pipeline) to tie it together reproducibly

What are three tools used in genetics research? GATK for variant discovery ³, Ensembl VEP for functional annotation, and PLINK 2 or Hail for population-scale association and polygenic score work. Gene-finding tools built on sequence features have a long history, from early systems like GRAIL onward ⁷, and the field’s tooling has been periodically surveyed as it grew ⁸.

What are genomic tools? Software and assays that measure or interpret genome-scale data: sequencers and arrays on the measurement side, and aligners, variant callers, annotators, and browsers on the interpretation side. Functional genomics adds RNA-seq, chromatin, and methylation assays with their own tool stacks ⁹.

What are the three main types of genetic testing? Targeted single-gene or panel testing, exome sequencing, and whole-genome sequencing. They differ in breadth and in what they can detect: panels are deepest and cheapest, exomes cover about 1-2% of the genome (the coding portion), and WGS covers coding plus non-coding, structural variants, and mitochondrial DNA. Clinical versions of all three run in accredited labs with confirmatory methods.

Can I use my consumer genotyping file instead? A 23andMe or AncestryDNA raw file is array genotypes at roughly 600,000-900,000 positions. The file holds no sequence. It is fine for ancestry and for many polygenic scores after imputation, and it cannot find rare or novel variants, indels, or structural variation, because the array only reports sites it was designed to probe.

Footnotes

1. Rachel L Goldfeder, Dennis P Wall, Muin J Khoury, et al. Human Genome Sequencing at the Population Scale: A Primer on High-Throughput DNA Sequencing and Analysis. American Journal of Epidemiology, 2017. https://doi.org/10.1093/aje/kww224 ↩

2. S. Pabinger, A. Dander, M. Fischer, et al. A survey of tools for variant analysis of next-generation genome sequencing data. Briefings in Bioinformatics, 2013. https://doi.org/10.1093/bib/bbs086 ↩

3. Aaron McKenna, Matthew Hanna, Eric Banks, et al. The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Research, 2010. https://doi.org/10.1101/gr.107524.110 ↩ ↩²

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

5. Sang Tae Park, Jayoung Kim. Trends in Next-Generation Sequencing and a New Era for Whole Genome Sequencing. International Neurourology Journal, 2016. https://doi.org/10.5213/inj.1632742.371 ↩

6. Alexis N. Burian, Wufan Zhao, Te‐Wen Lo, et al. Genome sequencing guide: An introductory toolbox to whole‐genome analysis methods. Biochemistry and Molecular Biology Education, 2021. https://doi.org/10.1002/bmb.21561 ↩

7. Edward C. Uberbacher, Ying Xu, Richard J. Mural. [16] Discovering and understanding genes in human DNA sequence using GRAIL. Methods in Enzymology, 1996. https://doi.org/10.1016/s0076-6879(96)66018-2 ↩

8. Andrew Fortna, Kathleen Gardiner. Genomic sequence analysis tools: a user’s guide. Trends in Genetics, 2001. https://doi.org/10.1016/s0168-9525(01)02229-6 ↩

9. Evelina Gasperskaja, Vaidutis Kučinskas. The most common technologies and tools for functional genome analysis. Acta medica Lituanica, 2017. https://doi.org/10.6001/actamedica.v24i1.3457 ↩