How to Use WGS Extract, and When to Use samtools Instead

By the end of this you will have: a verified alignment file with a known reference build, a coverage report you trust, a small-variant VCF you called yourself, mtDNA and Y haplogroup calls, and a microarray-format text file for third-party upload sites. You need the alignment file your sequencing provider gave you (BAM or CRAM, 15 to 80 GB depending on format), roughly 400 GB of free disk, 16 GB of RAM (32 GB is more comfortable), and a machine you can leave running overnight. Linux or macOS. On Windows, use WSL2 with the data on the Linux filesystem, not /mnt/c, because the 9p translation layer will cost you 3 to 5x on I/O-bound steps. Tools: samtools 1.19+, bcftools, mosdepth, Docker (for DeepVariant), and optionally WGS Extract v4 for the GUI-driven pieces.

Our view: WGS Extract is genuinely useful for two things, microarray-format exports and a fast look at haplogroups, and it bundles a reference genome library that saves you a painful download hunt. For everything else, run samtools and bcftools directly. You get flags, logs, exit codes, and the ability to script the whole thing. A GUI that hides the command it ran is a bad place to debug a reference mismatch.

1. Verify the file before you do anything to it

Providers ship truncated files more often than you would expect, usually from an interrupted download.

samtools quickcheck -v mygenome.cram && echo "EOF block OK"
md5sum mygenome.cram   # compare to the provider's manifest

quickcheck only reads the header and the EOF block, so it takes a second and catches truncation. It will not catch corruption in the middle. If the provider gave you an MD5, check it. If they did not, ask.

Then look at the header:

samtools view -H mygenome.cram | head -50
samtools view -H mygenome.cram | grep '^@PG'

The @PG lines tell you the aligner and, usually, the exact command line, including the path to the reference FASTA the provider used. That path is the single most useful string in the file.

2. Pin down the reference build exactly

“GRCh38” is not a specification. There are at least four widely used flavors and they are not interchangeable for CRAM decoding. Read the @SQ lines:

samtools view -H mygenome.cram | awk '$1=="@SQ"' | head -5

Decision table:

• SN:chr1 LN:248956422 means GRCh38 with UCSC-style names.
• SN:1 LN:249250621 means GRCh37 with Ensembl-style names (often hs37d5).
• SN:chr1 LN:249250621 means UCSC hg19.
• Presence of SN:chrUn_KI270302v1 or SN:HLA-A*01:01:01:01 means the full GRCh38 analysis set with alt contigs and HLA decoys.
• SN:chr14_KI270726v1_random without HLA lines usually means the no-alt analysis set.

Also check chrM length. GRCh38 and Ensembl GRCh37 both use rCRS at 16,569 bp. UCSC hg19 uses a different mitochondrial sequence at 16,571 bp. That difference breaks mtDNA haplogroup calling silently, and it is the single most common reason a Haplogrep result comes back nonsense.

The M5 tags in the @SQ lines are MD5 checksums of the actual sequence. Use them to confirm you have the right FASTA:

samtools dict GRCh38_full_analysis_set_plus_decoy_hla.fa | awk '$1=="@SQ"' | head -3

Compare the M5 values against the CRAM header. If they match, you have the exact reference. If they do not, you have a different file with the same name, and CRAM decoding will either fail loudly or fetch sequences from the EBI reference registry over the network, which is slow and depends on their uptime.

Set a local cache so it never goes to the network:

export REF_CACHE=$HOME/.cache/cram_refs/%2s/%2s/%s
export REF_PATH=$REF_CACHE
seq_cache_populate.pl -root $HOME/.cache/cram_refs GRCh38_full_analysis_set_plus_decoy_hla.fa

WGS Extract’s reference library covers the common builds and will do this lookup for you. If you are unsure what your provider used, opening the file in WGSE and reading its build detection is a reasonable five-minute shortcut. Then go back to the command line.

3. Decide whether to convert CRAM to BAM

Mostly: do not. CRAM at 30x is roughly 15 to 25 GB, BAM is 55 to 80 GB, and every modern tool reads CRAM given -T reference.fa. Converting costs you an hour and 50 GB for no analytical benefit.

Convert only when a specific tool refuses CRAM. If you must:

samtools view -@ 8 -b -T GRCh38_full_analysis_set_plus_decoy_hla.fa \
  -o mygenome.bam mygenome.cram
samtools index -@ 8 mygenome.bam

Expect 40 to 90 minutes at 8 threads on an SSD. CRAM decode is CPU-bound on the rANS/bzip2 codecs, so threads help nearly linearly up to about 8.

Going the other direction, to reclaim disk:

samtools view -@ 8 -C -T GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --output-fmt-option archive \
  -o mygenome.cram mygenome.bam

Lossless by default. Do not use --output-fmt-option lossy_names or quality binning unless you have decided you never want to recall variants.

If you need FASTQ (to realign to a different build, for example):

samtools collate -@ 8 -u -O mygenome.cram tmp_collate | \
  samtools fastq -@ 4 -1 R1.fq.gz -2 R2.fq.gz -0 /dev/null -s /dev/null -n -

collate is required. Running samtools fastq on a coordinate-sorted file without it produces interleaving that silently breaks read pairing. This is the failure mode we see most often in forum posts.

4. Get real coverage numbers

samtools coverage gives a per-chromosome summary in about two minutes:

samtools coverage --reference GRCh38_full_analysis_set_plus_decoy_hla.fa mygenome.cram

For distribution and callable-fraction numbers, use mosdepth:

mosdepth -n --fast-mode --by 1000 -t 4 \
  --fasta GRCh38_full_analysis_set_plus_decoy_hla.fa \
  sample mygenome.cram

Read sample.mosdepth.summary.txt and sample.mosdepth.global.dist.txt. What to look for on a “30x” genome:

• Mean autosomal depth of 28 to 34x. Below 25x, ask the provider what they promised. Providers quote raw yield, not aligned depth, and the gap is usually 8 to 12 percent after duplicate and unmapped reads.
• Fraction of the genome at ≥10x should be above 0.93 for PCR-free libraries. Below 0.88 suggests uneven coverage from a PCR-based prep.
• Duplicate rate under 10 percent. Check with samtools flagstat -@ 8.

Sex check from the same output: chrX depth near autosomal and chrY near zero outside pseudoautosomal regions indicates XX. chrX and chrY each near half autosomal depth indicates XY. Anything in between is worth a second look at the raw distribution before you conclude anything.

Contamination is the check people skip and should not:

VerifyBamID --SVDPrefix resource/1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --BamFile mygenome.cram --NumThread 8

FREEMIX above 0.02 means roughly 2 percent of reads come from another individual, which will inflate your heterozygous call count and produce false low-allele-fraction variants. Below 0.01 is clean.

5. Call variants yourself

Do not trust the provider’s VCF without knowing the caller and filter settings. Many DTC VCFs are produced with default GATK or a vendor pipeline tuned for a target other than yours, and some are filtered down to reported sites only.

DeepVariant is our default. It is a single container, it handles Illumina 30x PCR-free extremely well, and it produces calibrated quality scores without a separate VQSR step ¹.

docker run -v "$PWD":/data -w /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/mygenome.cram \
  --output_vcf=/data/mygenome.dv.vcf.gz \
  --output_gvcf=/data/mygenome.dv.g.vcf.gz \
  --num_shards=$(nproc) \
  --intermediate_results_dir=/data/dv_tmp

Budget 4 to 9 hours on 16 to 32 cores and about 150 GB of scratch in dv_tmp. Always emit the gVCF. The gVCF distinguishes “reference” from “no coverage here,” and that distinction matters every time you look up a specific position and find nothing.

Sanity-check the result before using it:

bcftools stats mygenome.dv.vcf.gz | grep -E '^SN|^TSTV'

Expect roughly 4 to 5 million variant sites against GRCh38, a transition/transversion ratio near 2.0 to 2.1 genome-wide, and a het/hom ratio around 1.5 to 1.6 for a non-African-ancestry sample (higher for African ancestry, reflecting greater heterozygosity). A Ti/Tv of 1.7 or below means false positives are leaking in. A variant count of 3 million or 6 million means something is wrong with the reference or the filtering.

Then normalize, because almost every downstream annotator assumes left-aligned, biallelic records:

bcftools norm -m -any -f GRCh38_full_analysis_set_plus_decoy_hla.fa \
  -Oz -o mygenome.norm.vcf.gz mygenome.dv.vcf.gz
bcftools index -t mygenome.norm.vcf.gz

On interpretation: a called variant is a measurement, and the distance between a measurement and a clinical conclusion is large. Most variants in any healthy genome are of uncertain significance, and the published experience with sequencing healthy adults is that incidental findings are common, frequently reclassified, and poorly served by self-interpretation ². Professional bodies have been explicit that whole-genome sequencing in a health context requires pre- and post-test counseling and a plan for unsolicited findings ³. If you find something that looks actionable, take it to a clinical geneticist and expect them to require confirmation in a CLIA-certified lab before anyone acts on it.

6. mtDNA and Y haplogroups

For mtDNA, extract chrM and call against rCRS:

samtools view -b -T GRCh38_full_analysis_set_plus_decoy_hla.fa \
  mygenome.cram chrM > chrM.bam
samtools index chrM.bam
bcftools mpileup -f GRCh38_full_analysis_set_plus_decoy_hla.fa -r chrM \
  -d 8000 -a AD chrM.bam | \
  bcftools call -mv --ploidy 1 -Ov -o chrM.vcf

The -d 8000 matters. Default mpileup depth caps at 250, and mtDNA typically sits at 1,000 to 5,000x because of copy number per cell. Capping depth biases heteroplasmy estimates and can drop real variants.

Upload chrM.vcf to Haplogrep 3, or run its CLI. Two caveats: if your build is UCSC hg19, the chrM sequence is not rCRS and you must realign those reads to rCRS first (WGS Extract does this conversion, which is a good reason to use it here). And heteroplasmy below about 10 percent is not reliably distinguishable from sequencing error at these depths without a caller designed for it.

For Y, yleaf reads the alignment directly:

yleaf -bam mygenome.bam -rg hg38 -o yleaf_out -r 1 -q 20 -t 8

It reports positions and a predicted haplogroup from the ISOGG/YFull marker sets. At 30x you will get a call several levels deeper than any microarray-based service, because arrays genotype a few thousand Y SNPs and your WGS covers the callable Y region at full depth.

7. Microarray-format export

This is where WGS Extract earns its place. Third-party sites (GEDmatch, DNA.Land-style tools, various genealogy uploads) accept 23andMe, AncestryDNA, FTDNA, and MyHeritage text formats, and not VCF. WGSE generates these from your alignment by genotyping the exact template positions each vendor uses.

Run the “Extract Data” tab, pick the output format, and let it work. It takes 20 to 60 minutes. Expect an output file with 600,000 to 950,000 rows depending on the template.

Three things to understand about the result:

1. Positions with insufficient coverage come out as --. A no-call rate of 1 to 3 percent is normal. Above 5 percent, revisit your coverage numbers from step 4.
2. Vendor formats encode genotypes on the plus strand of a specific build, usually GRCh37 for the older templates. WGSE handles the strand flip and coordinate conversion. Doing this by hand from a VCF is where most homemade converters go wrong, and the errors show up as systematically miscalled A/T and C/G sites.
3. These files carry no imputation. A real 23andMe file has had a genotyping-array-specific QC pass applied. Yours has not, so matching algorithms may behave slightly differently at the margins.

8. Lifting over between builds

Some tools only accept GRCh37. Convert the VCF rather than realigning:

CrossMap vcf hg38ToHg19.over.chain.gz mygenome.norm.vcf.gz \
  hg19.fa mygenome.hg19.vcf

Roughly 0.1 to 0.5 percent of variants will fail to map, mostly in regions restructured between builds. CrossMap writes these to a .unmap file. Read it. Do not liftover twice, and never liftover a lifted file back to the original build, because the round trip accumulates coordinate errors at every indel.

Common problems

CRAM decoding is slow or hangs. You have not set REF_CACHE and samtools is fetching reference sequences from EBI one contig at a time. Fix it with seq_cache_populate.pl.

samtools reports “reference mismatch” or garbage bases. Your FASTA is a different flavor of the same build. Compare M5 tags as in step 2. An hs37d5 FASTA will not decode a UCSC hg19 CRAM.

Variant count is 6 million or more. Usually contamination (check FREEMIX) or calling against alt contigs without alt-aware handling, which produces duplicate calls in the alt regions. Restrict to primary contigs: -r chr1,chr2,...,chrX,chrY,chrM.

Haplogroup comes back as a common root like H2a2a1 or R-M207 and goes no deeper. For mtDNA that is the rCRS-versus-hg19-chrM problem. For Y, check chrY depth first, since many pipelines mask or exclude chrY.

FASTQ extraction produces unpaired reads. You skipped samtools collate.

WGS Extract crashes on a large file. The v4 beta is memory-hungry during template genotyping. Close everything else, and if it still fails, run the microarray export on a BAM rather than a CRAM.

The provider’s VCF and yours disagree. They will, on roughly 1 to 2 percent of sites, mostly indels in homopolymer runs and low-complexity regions. Compare with bcftools isec and look at the discordant sites in IGV before deciding who is right. Neither caller is authoritative at those positions.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. 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. Noralane M. Lindor, Stephen N. Thibodeau, Wylie Burke. Whole-Genome Sequencing in Healthy People. Mayo Clinic Proceedings, 2017. https://doi.org/10.1016/j.mayocp.2016.10.019 ↩

3. Carla van El, Martina C. Cornel, Pascal Borry, et al. Whole-genome sequencing in health care. European Journal of Human Genetics, 2013. https://doi.org/10.1038/ejhg.2013.46 ↩