Analysing My 30× Whole-Genome Sequence

A whole-genome sequencing dataset is unusual personal data: compact enough to keep on a workstation, large enough to require a genuine processing pipeline and complicated enough that producing an answer is much easier than producing a trustworthy one.

I recently received the raw data from my own 30× whole-genome sequencing run and began building a reproducible workflow around it.

The immediate objective is not to search for interesting traits or medical conclusions. It is to understand the dataset, verify its quality and establish a reliable foundation for later analysis.

What 30× Sequencing Means

The dataset was produced using paired-end whole-genome sequencing at approximately 30× coverage.

Thirty-fold coverage means that each position in the genome has been read around thirty times on average. Coverage is not distributed perfectly evenly, but this depth is generally sufficient for reliable detection of many single-nucleotide variants and small insertions or deletions.

My delivery contained two compressed FASTQ files and a supplied VCF:

ulredcbg2673_1.fq.gz     59.33 GB
ulredcbg2673_2.fq.gz     51.25 GB
ulredcbg2673.vcf.gz     168.60 MB

The FASTQ files contain the original sequencing reads and quality scores. The VCF contains variants already identified by the sequencing provider.

That distinction matters.

The VCF is convenient for exploration, but it represents the output of somebody else’s reference genome, alignment process, variant caller, filters and assumptions. Keeping the raw reads means I can inspect those decisions and eventually reproduce the analysis independently.

Starting with the Raw Data

Before alignment or variant calling, I began with an inventory and basic quality-control stage.

Each FASTQ contained approximately 480 million reads of 150 bases:

reads per file          480,029,091
bases per file       72,004,363,650
read length                     150
gc content                  ~40.44%

The reported base quality was high:

q20 bases     99.12–99.19%
q30 bases     95.52–96.62%

A Q30 score corresponds to an estimated base-call error probability of one in one thousand. It does not prove that every region of the genome is represented correctly, but it is a useful first indication that the sequencing run itself produced high-quality reads.

These checks are deliberately unremarkable. At this stage, a boring result is a good result.

The important questions are whether both read files are complete, whether their identifiers match, whether the read counts agree and whether there are obvious quality, adapter or composition problems before the much more expensive processing stages begin.

Inspecting the Supplied Variant File

I also profiled the supplied VCF to understand what had already been called.

It contained approximately 4.53 million variants marked as passing the provider’s filters, including:

pass variants              4,531,106
biallelic snps             3,721,122
biallelic indels             734,335
heterozygous variants      2,698,912
homozygous-alt variants    1,832,194
multiallelic positions        76,194

The genotype distribution also included reference calls, missing genotypes and a smaller number of multiallelic genotypes.

These counts are useful as consistency checks, but they are not conclusions by themselves. A plausible number of variants does not establish that each individual call is correct.

I therefore also looked at genotype quality and read depth. Among non-reference variants marked as passing:

gq ≥ 20     4,105,123
gq ≥ 30     3,840,750

dp ≥ 20     4,103,506
dp ≥ 30     3,693,802

This gives me a baseline against which I can later compare variants generated by my own pipeline.

Rebuilding the Pipeline

The next stage is to process the reads independently against GRCh38.

The planned workflow is broadly:

compressed fastq
    → read quality control
    → reference preparation
    → alignment
    → sorting and duplicate marking
    → alignment quality control
    → germline variant calling
    → gvcf and joint-ready output
    → variant normalisation and filtering
    → annotation and analysis

The intention is not to reproduce the provider’s VCF byte for byte. Different aligners, callers, reference builds and filtering policies will naturally produce some differences.

Instead, I want to understand where the results agree, where they diverge and whether those differences can be explained.

That includes comparing:

• chromosome and region coverage;
• transition-to-transversion ratios;
• heterozygous-to-homozygous ratios;
• SNP and indel counts;
• genotype concordance;
• quality and depth distributions;
• variants unique to either call set.

A Data-Engineering Problem

Much of this project resembles a conventional data platform more than an interactive scientific notebook.

The raw input is over 110 GB. Alignment files can become considerably larger, and intermediate stages are expensive enough that accidentally repeating work is undesirable.

The workflow therefore needs:

• deterministic inputs and reference files;
• versioned container images;
• recorded command lines and tool versions;
• checksums for important artefacts;
• restartable processing stages;
• explicit intermediate and derived-data directories;
• enough logging to explain how an output was produced.

I initially considered adopting a larger workflow framework, but chose a simpler Docker-and-script-based design for the current scope.

The pipeline is divided into explicit numbered stages with completion markers and provenance information. This is sufficient for one sample, keeps the execution model transparent and avoids introducing workflow-engine complexity before it provides a clear benefit.

Reusable processing logic lives outside the notebooks. Marimo provides the interactive interface for reviewing results, visualising distributions and documenting findings, rather than becoming the place where the full pipeline is implemented.

Compute and Storage

I am running the work as part of my Daystrom lab environment, using local NVMe storage and containerised tools.

The DGX Spark is available for parts of the project, although genomics does not automatically become a GPU workload simply because a GPU is present.

Many established bioinformatics tools remain primarily CPU-, memory- or storage-bound. Alignment and variant calling also have different bottlenecks from large-model inference. Fast sequential and random storage access, sufficient working space, careful parallelism and predictable restart behaviour matter at least as much as accelerator availability.

The Spark is nevertheless useful as a flexible system for later experiments involving machine learning, embedding scientific literature or running local language models alongside the analysis.

The broader lesson is similar to my experience with local AI inference:

Having enough compute to run a workload does not remove the need to understand its bottlenecks.

Interpretation Comes Later

A personal genome invites premature interpretation.

It is easy to search a variant database, find an association and treat that association as an explanation. The underlying evidence is usually more complicated.

A variant may have a small effect, apply only to a particular population, appear through linkage with another variant or come from a study that has not replicated well. Complex traits are generally influenced by many variants together with environmental and developmental factors.

There are also classes of variation that a conventional short-read small-variant pipeline handles less reliably, including:

• structural variants;
• copy-number changes;
• repeat expansions;
• low-level mosaicism;
• mitochondrial heteroplasmy;
• difficult or highly homologous genomic regions.

For those reasons, I am treating scientific evidence, genome build, allele orientation, population context and study quality as part of the data model rather than as details to resolve afterwards.

This is an engineering and research project, not a clinical diagnostic process.

What I Want to Explore

Once the independently generated call set and QC work are complete, I want to explore several areas.

The first is concordance: how closely my own processing agrees with the supplied VCF and what explains the differences.

The second is ancestry and population structure using suitable public reference datasets. I have previously worked with genotype matrices and 1000 Genomes data, which provides a useful basis for building this carefully rather than relying on opaque consumer-genetics summaries.

The third is association research around complex traits, including dyslexia. That work needs particular restraint. The objective is to review published evidence and understand how polygenic associations are constructed, not to infer a deterministic explanation from a handful of variants.

I am also interested in using the dataset to study techniques such as:

• variant annotation;
• genotype storage and indexing;
• population-scale comparison;
• polygenic-score methodology;
• structural-variant detection;
• reproducible scientific reporting.

Why Analyse It Myself

Commercial reports are useful, but they compress a large amount of processing and interpretation into a small number of conclusions.

I am more interested in the path from raw reads to those conclusions.

Working from FASTQ data forces me to confront the reference genome, alignment choices, sequencing uncertainty, variant representation and evidence quality. It turns the genome from a report into a dataset whose provenance can be examined.

The goal is not to replace geneticists or clinical laboratories.

It is to understand what the data can support, what it cannot support and how much engineering sits between sequencing a sample and making a defensible claim about it.

For now, success means producing a clean, reproducible and well-understood variant dataset.

The interesting questions can come afterwards.