From FASTQ to Variants: The Bioinformatics ETL Pipeline

Under scientific review

This introductory page is under scientific review. Use it for orientation, not clinical decision-making.

In cancer genomics, turning raw physical DNA into clinical intelligence requires a multi-stage ETL (Extract, Transform, Load) pipeline. This pipeline ingests gigabytes of unaligned raw sequencer outputs (.fastq) and processes them into structured, annotated variant call files (.vcf / .maf).

For data engineers, this process is conceptually identical to a distributed log ingestion pipeline: raw unstructured logs $\rightarrow$ parsed & indexed logs $\rightarrow$ diff analysis $\rightarrow$ alert enrichment $\rightarrow$ filtering.

---

The Core Pipeline Architecture

[ FASTQ (Raw Reads) ] 
       │
       ▼ (FastQC / trimmomatic)
[ Clean FASTQ ] 
       │
       ▼ (BWA-MEM / STAR -> Reference GRCh38)
[ SAM / BAM (Aligned Reads) ] 
       │
       ▼ (Picard MarkDuplicates + GATK BQSR)
[ Recalibrated BAM ] 
       │
       ▼ (GATK Mutect2 / HaplotypeCaller)
[ Raw VCF (Variant Call Format) ] 
       │
       ▼ (VEP / ANNOVAR)
[ Annotated VCF ] 
       │
       ▼ (Filter: VAF, Depth, gnomAD)
[ Actionable Variants (MAF / Report) ]

---

1. FASTQ: The Raw Log File

A .fastq.gz file contains millions of raw sequencing reads. Each read consists of 4 standardized lines:

@A00142:315:H7YYFDSX2:1:1101:1245:1000 1:N:0:GATCAG
ATCGGCTAAGCTAGGCTAACGATCGATCGATCGATCGATCGATCGAT
+
FFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ

• Line 1 (@ID): Instrument ID, flow cell lane, tile coordinates, and sample barcode index.
• Line 2 (Sequence): The raw nucleotide base calls (A, C, T, G, or N for ambiguous).
• Line 3 (+): Separator.
• Line 4 (Phred Quality Scores): ASCII-encoded quality score ($Q$) for each base.

Phred Quality Scores Formula

$$Q = -10 \log_{10} (P_{\text{error}})$$

• Q20: $1%$ error probability ($99%$ accuracy).
• Q30: $0.1%$ error probability ($99.9%$ accuracy — standard clinical threshold).

---

2. Alignment: Indexing the Log Stream

Short reads (150 bp) carry no spatial context. Alignment maps each read to its exact coordinate on the 3-billion-base Human Reference Genome (GRCh38).

• Tools: BWA-MEM (DNA-seq), STAR (RNA-seq).
• Algorithm: Uses the Burrows-Wheeler Transform (BWT) and FM-index — a compressed full-text index that enables $O(m)$ string matching across billions of bases.
• Output Formats:
  • SAM (Sequence Alignment/Map): Human-readable plain text format.
  • BAM: Binary, block-compressed, indexed SAM file (.bam + .bai). Supports random access queries (like samtools view region chr17:7668402-7687550).
  • CRAM: Reference-based ultra-compressed format (~50% smaller than BAM).

---

3. Post-Alignment Cleaning (GATK Best Practices)

Raw BAM files contain technical noise that distorts variant calling:

1. Mark Duplicates (Picard / sambamba): Identifies PCR duplicates (identical DNA molecules amplified during library prep). Duplicate reads are flagged so they aren't counted multiple times as independent evidence.
2. Base Quality Score Recalibration (BQSR): Machine-learning model that adjusts raw Phred quality scores based on systematic optical or chemical hardware errors.

---

4. Variant Calling: The Genomic diff

Variant callers compare the recalibrated BAM file against the reference genome to identify differences:

• Germline Calling (GATK HaplotypeCaller): Evaluates inherited variants (present in 50% or 100% of reads, corresponding to heterozygous or homozygous alleles).
• Somatic Tumor-Normal Calling (Mutect2 / Strelka2): Compares Tumor.bam vs. Normal_Blood.bam from the same patient. This subtracts background germline variants to isolate acquired somatic mutations present at low Variant Allele Fractions (VAF) (e.g., 2% to 10% in heterogeneous tumors).

VCF Format (Variant Call Format)

#CHROM POS      ID        REF ALT QUAL FILTER INFO
chr7   140753336 rs121913364 A   T   .    PASS   GENE=BRAF;VAF=0.42;AA=V600E

---

5. Annotation & Filtering: Enriching the Alert

A raw VCF lists thousands of candidate variants. Annotation enriches each record with biological context:

• Functional Impact (VEP / SnpEff): Maps variant location to gene models $\rightarrow$ predicts consequence (Missense, Nonsense, Frameshift, Synonymous).
• Population Frequency (gnomAD): Attaches global allele frequencies. Variants present at $>1%$ in healthy populations are discarded as benign polymorphisms.
• Clinical Knowledgebases: Cross-references ClinVar, COSMIC, and Civic for known oncogenic pathogenicity and FDA drug associations.

---

6. Pipeline Execution: Production Workflows

In production clinical environments, these pipelines are orchestrated using containerized workflow managers:

• Nextflow (nf-core/sarek) or Cromwell (WDL).
• Running on cloud infrastructure (AWS Batch, GCP Life Sciences) or HPC slurm clusters.

---

References

1. Van der Auwera, G. A., & O'Connor, B. D. (2020). Genomics in the Cloud: Using Docker, GATK, and WDL in Terra. O'Reilly Media.
2. Li, H. (2013). Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv:1303.3997.
3. Cibulskis, K., et al. (2013). Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nature Biotechnology, 31(3), 213-219.