How to Run a Whole Genome Sequencing Analysis Pipeline on Your Own Data

By the end of this you will have the following from a pair of FASTQ files.

• a CRAM aligned to GRCh38
• a small-variant VCF with roughly 4.5–5 million SNVs and 800k–1M indels
• a structural variant VCF with a few thousand calls
• a per-base coverage profile
• an annotated TSV you can load into pandas or DuckDB and query

You need: paired-end 30x WGS FASTQs (about 100–120 GB gzipped), a Linux box with 32+ cores and 64 GB RAM (128 GB if you want to build the bwa-mem2 index yourself), 1.5 TB of free disk, Docker or Singularity, and about 12–20 hours of wall clock. Everything below is command-line. Nothing here is a diagnosis, and several steps end with “take this to a clinical geneticist,” which is the correct ending.

1. Verify what you were shipped before you compute on it

Check the MD5s the lab gave you. If they did not give you MD5s, ask.

md5sum -c checksums.md5

Then confirm the FASTQs are intact and paired. A truncated gzip stream will pass a naive read and fail 6 hours into alignment.

for f in sample_R1.fastq.gz sample_R2.fastq.gz; do
  gzip -t "$f" && echo "$f ok"
done
seqkit stats -a -j 8 sample_R*.fastq.gz

seqkit stats -a gives you read count, mean length, and Q20/Q30 rates. For a 2x150 NovaSeq run expect Q30 above 90% and read counts within 1% of each other between R1 and R2. Compute your expected coverage before you trust the lab’s claim: total bases divided by 3.1e9. If you were sold 30x and seqkit says 75 Gbp total, you have 24x, and that matters because indel sensitivity degrades noticeably below 25x.

Also check the read headers for the instrument and lane. You need them for the read group, and multiple lanes mean multiple FASTQ pairs that should be aligned separately and merged. Do not concatenate them blindly.

2. Pick a reference and stop thinking about it

Use the GRCh38 no-ALT analysis set. The plain GRCh38 primary assembly will not do, nor will hg19, nor a UCSC hg38.fa you found in a Dropbox link.

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
gunzip GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
samtools faidx GCA_000001405.15_GRCh38_no_alt_analysis_set.fna

The no-ALT set matters because ALT contigs are highly similar to the primary chromosomes. An aligner that is not ALT-aware will split reads across them, mapping quality collapses to 0, and you silently lose variant calls in regions like the MHC. The analysis set also hard-masks the pseudoautosomal regions on chrY and includes the EBV decoy, which soaks up reads that would otherwise mismap.

T2T-CHM13v2.0 is a better reference by most measures, but the tooling ecosystem is still GRCh38-native. That ecosystem includes gnomAD frequencies, ClinVar coordinates, VEP caches and most polygenic score weight files. Our recommendation: run GRCh38 as the primary, and if you care about specific regions (segmental duplications, centromeric genes, the full KIR or SMN loci) do a second pass against CHM13 rather than switching wholesale.

3. Align with bwa-mem2, mark duplicates, write CRAM

Building the bwa-mem2 index takes about 80 GB of RAM and an hour. If you have less, use bwa mem (roughly 2x slower at runtime but the index needs about 6 GB) or download a prebuilt index.

bwa-mem2 index GCA_000001405.15_GRCh38_no_alt_analysis_set.fna

Align, fixmate, sort, mark duplicates, and emit CRAM in one stream. Writing an intermediate BAM costs you 80 GB of disk and an extra hour of I/O for nothing.

REF=GCA_000001405.15_GRCh38_no_alt_analysis_set.fna
RG='@RG\tID:L001\tSM:sample01\tLB:lib1\tPL:ILLUMINA\tPU:HXXXX.1'

bwa-mem2 mem -t 28 -K 100000000 -Y -R "$RG" "$REF" \
    sample_R1.fastq.gz sample_R2.fastq.gz \
  | samtools fixmate -m -@ 4 -u - - \
  | samtools sort -@ 8 -m 3G -u -T /scratch/sorttmp - \
  | samtools markdup -@ 8 --reference "$REF" --output-fmt CRAM \
      -f markdup_stats.txt - sample01.cram

samtools index -@ 8 sample01.cram

Flag notes. -K 100000000 fixes the input chunk size so the output is deterministic regardless of thread count, which you want if you ever need to reproduce a call. -Y uses soft clipping for supplementary alignments, which Manta and most SV callers expect. samtools fixmate -m adds the mate score tags that markdup needs. The -u flags keep the pipe uncompressed between stages, saving a lot of CPU.

The CRAM will be 15–20 GB at 30x. Keep it, delete the FASTQs only after you have verified you can regenerate reads from the CRAM (samtools fastq) and only if you are comfortable losing original quality-score binning. We keep the FASTQs.

4. QC the alignment, and take the numbers seriously

Four checks. Run them all, they take minutes.

samtools stats -@ 8 --reference "$REF" sample01.cram > sample01.stats
mosdepth -t 8 -n --fast-mode --by 1000 sample01 sample01.cram

From samtools stats, pull:

• reads mapped / raw total sequences should be above 0.99 against the no-ALT set with decoy.
• reads properly paired above 0.95.
• insert size average typically 350–450 for a PCR-free library, with standard deviation under 120. A bimodal or very wide insert distribution means library problems and will produce false SV calls.
• error rate under 0.005.
• Duplicate fraction from markdup_stats.txt: under 5% for PCR-free, 10–20% for PCR-positive libraries. Above 25% means you paid for coverage you did not get.

From mosdepth, check sample01.mosdepth.summary.txt for genome-wide mean and the *.thresholds or quantized output for the fraction of the callable genome at ≥10x and ≥20x. Target: ≥95% of the autosomes at 10x or more. Coverage below that concentrates in GC-extreme regions, and the first exon of many genes is GC-rich.

Contamination, which is the check people skip:

VerifyBamID --SVDPrefix resource/1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference "$REF" --BamFile sample01.cram --NumThread 8

FREEMIX above 0.02 means roughly 2% of your reads come from another person. At that level you start seeing spurious heterozygous calls at low allele fraction. Above 0.03, stop and get the sample re-sequenced.

Finally, confirm chromosomal sex from coverage: the chrX/autosome mean depth ratio from mosdepth should be near 1.0 or near 0.5, and chrY near 0 or near 0.5. A ratio that sits in between is a sample identity or aneuploidy signal and is worth investigating before anything downstream.

5. Call small variants with DeepVariant

DeepVariant outperforms hand-tuned GATK hard filters on a single germline sample, requires no known-sites recalibration step, and is one command. That is the whole argument.

docker run --rm -v "$PWD":/data -v /scratch:/scratch \
  google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
    --model_type=WGS \
    --ref=/data/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna \
    --reads=/data/sample01.cram \
    --output_vcf=/data/sample01.dv.vcf.gz \
    --output_gvcf=/data/sample01.dv.g.vcf.gz \
    --num_shards=28 \
    --intermediate_results_dir=/scratch/dv_tmp

On 28 cores this runs 3–6 hours for 30x WGS. Emit the gVCF even if you have no plans for it. It is the only way to distinguish “reference at this position” from “no data at this position,” and you will want that distinction the first time you check whether you carry a specific variant and the VCF is silent.

Then check the output before you look at a single gene:

bcftools stats -F "$REF" sample01.dv.vcf.gz > sample01.vcfstats
grep -E "^SN|^TSTV" sample01.vcfstats

Expected for a single human genome against GRCh38: 4.4–5.0M SNVs and 800k–1M indels. Genome-wide Ti/Tv should be 2.0–2.1, with a het/hom ratio of 1.4–1.7. SNV count varies with ancestry, with African-ancestry genomes at the high end. Ti/Tv below 1.9 means false positives are leaking in. A het/hom ratio far below 1.4 in a non-consanguineous sample usually means allelic dropout from low coverage.

If you want GATK instead, run HaplotypeCaller -ERC GVCF per interval. Follow it with GenotypeGVCFs, then VQSR with the GATK resource bundle. It is more steps, needs the known-sites VCFs, and VQSR is poorly calibrated on a single sample. We do not recommend it for n=1.

6. Call structural variants and repeat expansions separately

Small-variant callers do not find deletions longer than about 50 bp. Structural variants are a meaningful share of the functionally interesting genome, and the standard single-sample combination is Manta plus a read-depth caller.

configManta.py --bam sample01.cram --referenceFasta "$REF" --runDir manta_run
manta_run/runWorkflow.py -j 16

Output lands in manta_run/results/variants/diploidSV.vcf.gz. Expect 5,000–10,000 calls, most of them small deletions and insertions in repetitive regions. Manta is precise on balanced events and breakpoint-resolved deletions and weak on large copy-number changes, so add a depth-based caller if CNVs matter to you.

Short tandem repeats need their own tool because the pathogenic ones (the CAG, CGG, and GAA expansions) are longer than the read and invisible to everything above:

ExpansionHunter --reads sample01.cram --reference "$REF" \
  --variant-catalog variant_catalog/grch38/variant_catalog.json \
  --output-prefix sample01.eh

Short-read ExpansionHunter genotypes are reliable for repeats shorter than the fragment length and give a lower bound beyond that. Any expansion near or above a clinical threshold needs orthogonal confirmation by a laboratory method. This is a place where a short-read WGS result is a screen. The result is not an answer.

7. Annotate with VEP against a local cache

Run VEP offline. The REST API is rate-limited and will take days for 5M variants.

vep --offline --cache --dir_cache $HOME/.vep --assembly GRCh38 \
  --fasta "$REF" --fork 16 --vcf --compress_output bgzip \
  --input_file sample01.dv.vcf.gz --output_file sample01.vep.vcf.gz \
  --everything --check_existing --pick_allele_gene \
  --custom $HOME/annot/clinvar_20250601.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT,CLNDN \
  --plugin dbNSFP,$HOME/annot/dbNSFP4.7a_grch38.gz,REVEL_score,CADD_phred,AlphaMissense_score \
  --plugin SpliceAI,snv=$HOME/annot/spliceai_scores.raw.snv.hg38.vcf.gz,indel=$HOME/annot/spliceai_scores.raw.indel.hg38.vcf.gz

--pick_allele_gene collapses to one consequence per allele per gene, which keeps the output at a size you can work with. --everything turns on the frequency, regulatory, and canonical-transcript flags in one go.

Allele frequency is your primary filter, and the reference panel you filter against has to reflect your ancestry. Deep WGS of 1,070 Japanese individuals identified millions of variants absent from the global databases of the time, a large fraction of them population-specific ¹. If you filter with a panel that does not include people like you, common benign variants in your population will show up as ultra-rare and land at the top of your candidate list.

8. Reduce to something you can read

bcftools +split-vep sample01.vep.vcf.gz -f \
  '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t%SYMBOL\t%Consequence\t%IMPACT\t%gnomADg_AF\t%ClinVar_CLNSIG\t%REVEL_score\t%CADD_phred\t%SpliceAI_pred[\t%GT\t%DP\t%VAF]\n' \
  -d -A tab | gzip > sample01.annot.tsv.gz

Load it and apply three passes.

Pass one, ClinVar: keep variants with CLNSIG containing Pathogenic or Likely_pathogenic and CLNREVSTAT of at least two submitters with no conflicts. In a healthy adult genome this yields a handful of carrier-state findings and occasionally something in an ACMG secondary-findings gene. Cross-reference the ACMG SF v3.2 list (81 genes, mostly cancer predisposition, cardiomyopathy, arrhythmia, and familial hypercholesterolemia) because those are the ones with established management pathways ².

Pass two, loss of function: IMPACT == "HIGH" and gnomAD popmax AF < 0.001, filtered to genes where loss of function is a known disease mechanism. Most HIGH-impact calls in any genome are annotation artifacts on the last exon or on transcripts nobody expresses. Check the variant in IGV before you believe it.

Pass three, pharmacogenomics: run PharmCAT or PyPGx on the VCF rather than reading star alleles out by hand. Diplotype calling for CYP2D6 in particular depends on structural variation and hybrid alleles that a position-level VCF cannot represent. PharmCAT output is a report, and any medication decision based on it belongs with your prescriber. That decision does not belong with you and a terminal.

Everything else is a variant of uncertain significance. There will be thousands. The VUS burden is the main practical obstacle to clinical WGS interpretation, and reclassification over time is normal and expected ³. Store your annotated table with the ClinVar release date and re-annotate every 6–12 months. Variants move categories.

9. Validate the pipeline itself

If you have the option, run one GIAB sample (HG002 is the usual choice) through exactly the commands above and compare to the truth set:

hap.py HG002_GRCh38_truth.vcf.gz HG002.dv.vcf.gz \
  -f HG002_GRCh38_confident.bed -r "$REF" -o happy_out --engine=vcfeval

On 30x Illumina you should see SNV F1 above 0.995 and indel F1 above 0.99 inside the confident regions. If your numbers are much lower, your genome is not the reason. The problem is in your pipeline. This is the single highest-value step people skip.

Common problems

Coverage looks fine but calls are sparse on chrX. Check whether your sample is XY and whether DeepVariant was run in a haploid-aware mode. Male chrX outside the PAR is haploid, and diploid calling there produces low-quality het calls.

MAPQ 0 pileups in the MHC, SMN1/SMN2, PMS2, and the HBA cluster. These are segmental duplications. Short reads cannot resolve them and you will get both false positives and false negatives. Treat any finding in those loci as unresolved until confirmed by a targeted method.

A gene shows zero coverage. Before concluding deletion, check the mosdepth per-base output and whether the gene is on an ALT contig or in a region masked in the analysis set. Real homozygous deletions have clean breakpoints and flanking discordant read pairs, both of which Manta will have called.

Ti/Tv of 1.7 and 6M SNVs. You aligned to the wrong reference build or mixed builds between alignment and calling. Coordinates will look plausible and results will be garbage.

The VCF has no VAF field. DeepVariant writes VAF in recent versions but older ones do not. Derive it from AD instead: bcftools query -f '[%AD]' and divide.

You found something that looks serious. Confirm it visually in IGV, check it against a second caller, then take it to a clinical geneticist or genetic counselor with the raw data. Research-grade WGS is not a diagnostic test, and confirmatory testing in an accredited laboratory is the standard path before anything acts on a result ³. Diagnostic interpretation is a clinical process that depends on phenotype, family history, and segregation, none of which are in your VCF ².

Footnotes

1. Masao Nagasaki, Jun Yasuda, Fumiki Katsuoka, et al. Rare variant discovery by deep whole-genome sequencing of 1,070 Japanese individuals. Nature Communications, 2015. https://doi.org/10.1038/ncomms9018 ↩

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

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