Skip to content

Which genome software should you run on your own WGS data?

Oak
A tall laboratory machine drawing a glowing strand through five colored glass stages, where fine threads align against a reference and a few flare red.

If you have your own 30x whole-genome data and want to work with it, the stack is short: align with bwa-mem2 (or use the vendor’s DRAGEN output if you got one), call small variants with DeepVariant, normalize and filter with bcftools, annotate with Ensembl VEP, and look at anything interesting in IGV. Then run a handful of specialized callers for the regions the general pipeline gets wrong. Geneious, GenomeTools, and the UCSC Genome Browser are all real software. The first two are aimed at molecular biology and genome annotation workflows, and UCSC is a browser. It is not a pipeline. None of them takes you from FASTQ to a variant list you’d trust.

The stack we’d run

Start from the rawest thing you have. If your provider gives you FASTQ or unaligned BAM, keep it. If they only give you a VCF, you are stuck with their filtering decisions and you cannot distinguish “reference” from “no coverage here.”

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

-K 100000000 fixes the chunk size so the output is deterministic across thread counts. -Y uses soft clipping for supplementary alignments, which downstream structural-variant callers want.

Then small variants:

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

Ask for the gVCF even if you think you only want the VCF. The gVCF carries per-base confidence, which is the only way to answer “does my data cover this position, or is it absent from the call set because the reads never got there?” That question comes up constantly once you start looking up specific variants you read about.

On a 32-core machine, alignment is roughly two to four hours for a 30x genome, DeepVariant another two to five. Budget disk before compute: 30x paired-end FASTQ is about 80-110 GB gzipped, the sorted BAM is 60-90 GB, and CRAM with the reference on hand is 15-20 GB. Convert and keep CRAM:

samtools view -@ 8 -T GRCh38.fa -C -o sample.cram sample.bam

Then normalize before you compare against anything:

bcftools norm -m -any -f GRCh38.fa -Ou sample.vcf.gz \
| bcftools norm -d exact -Oz -o sample.norm.vcf.gz
bcftools index -t sample.norm.vcf.gz

Unnormalized indels are the single most common reason a variant you know you have fails to match a database record.

Reference genome: GRCh38 or T2T-CHM13

Use GRCh38 with the analysis set that masks false duplications and includes decoy contigs (the no_alt_analysis_set FASTA). Every public annotation resource is keyed to it, including ClinVar, gnomAD and the VEP cache. T2T-CHM13 is a better reference in the sense that it is complete, and alignment against it resolves segmental duplications and centromeric regions that GRCh38 simply lacks, but the annotation ecosystem is still catching up. Our practice: GRCh38 as the primary, CHM13 as a second alignment when a specific locus looks messy.

Do not accept hg19 output and lift it over unless you have no choice. Liftover drops and mis-maps a small percentage of sites, and it does so in exactly the regions that are already hard.

What a general-purpose caller misses

This is the part that separates a usable personal genome from a VCF you glance at once. Short-read alignment against a linear reference fails predictably in paralogous and repetitive sequence, which has been the known hard edge of human genome sequencing since short reads took over.1

  • Structural variants and CNVs: run Manta plus a read-depth caller. DeepVariant will not report a 40 kb deletion.
  • Short tandem repeat expansions: ExpansionHunter with the known-locus catalog. Repeat expansions are invisible to standard calling because the reads do not span them.
  • Pharmacogenes: CYP2D6 sits next to a pseudogene and needs a dedicated caller (Cyrius). PharmCAT will take a VCF plus that call and produce star-allele diplotypes. Interpretation of those results belongs to a clinician. It is not your job or an LLM’s.
  • SMN1/SMN2, PMS2/PMS2CL, HBA1/HBA2 and GBA/GBAP1 are near-identical paralogs. Standard calls in these genes are unreliable, and a clinical lab uses orthogonal assays for them.
  • Mitochondrial DNA: needs its own calling mode for heteroplasmy. Mutect2 in mitochondrial mode works.

Quality control you should run before believing anything

mosdepth -n --fast-mode --by 1000 -t 4 sample sample.cram
verifybamid2 --SVDPrefix 1000g.phase3.100k.b38.vcf.gz.dat --Reference GRCh38.fa --BamFile sample.cram
bcftools stats sample.norm.vcf.gz > stats.txt

Sanity checks on a good 30x human genome: mean coverage 28-35x with 90%+ of the genome at 20x or better, contamination estimate under 1-2%, transition/transversion ratio around 2.0-2.1 genome-wide, and roughly 4-5 million variant sites against GRCh38. If your SNV count comes in at 3 million or 7 million, something is wrong with the sample, the alignment, or the filtering. Also confirm the reported sex chromosomes match expectation with somalier, which catches sample swaps fast.

Annotation without drowning

vep --cache --offline --assembly GRCh38 --fork 8 \
  --vcf --everything --canonical --pick_allele_gene \
  --plugin CADD,whole_genome_SNVs.tsv.gz \
  -i sample.norm.vcf.gz -o sample.vep.vcf.gz --compress_output bgzip

Four to five million variants is too many to read. The useful filter chain is: restrict to coding and splice consequences, then to gnomAD allele frequency below some threshold, then intersect with ClinVar for anything with a submitted interpretation. That typically leaves a few hundred rows, most of them benign or of uncertain significance. Uncertain significance means uncertain. A variant labeled pathogenic in ClinVar by one submitter with no assertion criteria is weak evidence, and acting on any of this requires a clinician and usually a confirmatory clinical-grade test, because research-grade WGS is not a diagnostic.

Where the tools that rank for “genome software” fit

GenomeTools is a solid, fast toolkit for GFF3 manipulation, annotation sorting, and sequence indexing. You will use it if you are building genome annotations. You will not use it if you are analyzing your own variants. Geneious is commercial desktop software for cloning, primer design, Sanger assembly, and phylogenetics, with a free viewer tier rather than a free full version. The UCSC Genome Browser is where you go after you have a coordinate and want to see conservation, regulatory tracks, and gene models around it. All three are worth knowing. None of them is your pipeline.

The reason this is a personal-computer problem at all is cost. The $1000 genome was a stated engineering target twenty years ago, and by the time short-read platforms hit production scale the bottleneck moved from generating sequence to storing and interpreting it.23 Sequencing is cheap. Making sense of the output is the work.

Questions people also ask

What is the best software for DNA sequencing? For short-read human WGS analysis, the choices are bwa-mem2 or DRAGEN for alignment and DeepVariant or GATK HaplotypeCaller for small variants. Manta and ExpansionHunter cover what those miss, and VEP handles annotation. For long reads, swap in minimap2 and Sniffles2. “Best” depends entirely on read type and question.

Are humans 99 percent the same? Yes, roughly. Any two people match at about 99.9% of single-base positions, which still leaves each genome differing from the reference at 4-5 million sites. Structural variation affects far more total bases than SNVs do, so the “99.9%” figure understates how different two genomes are in practice.

Is there a free version of Geneious? Geneious offers a free viewer with limited functionality and time-limited trials of the full product. For most personal-genome work the open-source stack above covers the same ground at no cost.

What is genome technology? The current production standard is sequencing-by-synthesis short reads, which brought per-genome cost down by orders of magnitude and moved whole-genome sequencing from consortium projects to routine use.3 Long-read platforms now complement it by resolving structural variants and repetitive regions that short reads cannot span.

Footnotes

  1. Daniel C. Koboldt, Li Ding, Elaine R. Mardis, et al. Challenges of sequencing human genomes. Briefings in Bioinformatics, 2010. https://doi.org/10.1093/bib/bbq016 ↩

  2. Simon T. Bennett, C.L. Barnes, Anthony J. Cox, et al. Toward the $1000 Human Genome. Pharmacogenomics, 2005. https://doi.org/10.1517/14622416.6.4.373 ↩

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

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, blood RNA sequencing, and plasma proteomics, delivered as a detailed report plus every data file. Build your profile.