Module 4: Oncology File Formats, Phred Scores & Alignment Concepts: FASTQ, BAM & VCF
1. FASTQ Format & Phred Quality Score Mechanics
Next-Generation sequencing instruments generate raw read outputs in the standardized 4-line FASTQ format:
@NB501234:245:H7G75BGX2:1:11101:10023:1045 1:N:0:TCCGCGAA
ATGCTTCGATCGGATCGATCGATCGATCGATCGATCGATCGATCGATCGATC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFLine 4 encodes base-calling quality via ASCII characters mapped to Phred Quality Scores ($Q$) using the Sanger Phred+33 encoding:
$$Q = -10 \cdot \log_{10}(P_{\text{error}}) \implies P_{\text{error}} = 10^{-Q/10}$$| Phred Score ($Q$) | Probability of Incorrect Base Call | Base Call Accuracy | Clinical Interpretation |
|---|---|---|---|
| $Q_{10}$ | 1 in 10 ($10\%$) | 90.0% | Unacceptable clinical error rate; trimmed during QC. |
| $Q_{20}$ | 1 in 100 ($1.0\%$) | 99.0% | Minimum acceptable raw sequencing threshold. |
| $Q_{30}$ | 1 in 1,000 ($0.1\%$) | 99.9% | CLIA/CAP Clinical Benchmark: High-confidence base call. |
| $Q_{40}$ | 1 in 10,000 ($0.01\%$) | 99.99% | Ultra-high accuracy (typical of consensus UMI error-corrected reads). |
2. SAM, BAM & CRAM Alignment Architecture
Sequence reads are aligned against a reference human genome assembly (GRCh37/hg19 or GRCh38/hg38) using algorithms like BWA-MEM. The resulting alignments are stored in Sequence Alignment/Map (SAM) or compressed binary BAM / CRAM formats.
Decoding SAM Bitwise Flags
Column 2 of every BAM record contains an integer representing bitwise flags that describe read pairing and mapping properties:
# Common SAM Bitwise Flag Bitmask Values:
# 1 (0x1) -> Read paired in sequencing
# 2 (0x2) -> Read mapped in proper pair (correct distance and orientation)
# 4 (0x4) -> Read unmapped
# 16 (0x10) -> Read reverse strand
# 64 (0x40) -> First in pair (R1)
# 128 (0x80) -> Second in pair (R2)
# 256 (0x100)-> Secondary alignment
# 1024(0x400)-> PCR or optical duplicateFor example, a read with flag 99 ($64 + 32 + 2 + 1$) indicates: Read paired (1) + mapped in proper pair (2) + mate on reverse strand (32) + first in pair (64).
3. VCF (Variant Call Format) Architecture in Somatic Oncology
The VCF v4.2 standard records identified genomic alterations. Unlike germline VCFs that simply report genotypes (0/1, 1/1), somatic oncology VCFs embed depth and allele frequency metrics across INFO (variant-level metadata) and FORMAT / SAMPLE (sample-specific call data) columns:
| VCF Field | Clinical Representation | Real-World Example in Oncology VCF |
|---|---|---|
#CHROM | Chromosome identifier | chr7 |
POS | 1-based genomic position | 55259515 |
ID | Known database identifier (dbsnp / COSMIC) | COSM6224;rs121434568 |
REF | Reference base(s) | T |
ALT | Alternate (mutant) base(s) | G |
QUAL | Phred-scaled probability of variant existence | 852.77 |
FILTER | Filter status (PASS or reason for rejection) | PASS (or artifact_detection;low_vaf) |
INFO | Variant-level annotations (gene, consequence) | DP=642;ECNT=1;POP_AF=0.00001;TLOD=148.5 |
FORMAT | Data structure format string for sample column | GT:AD:AF:DP:F1R2:F2R1 |
SAMPLE | Tumor sample metrics corresponding to FORMAT | 0/1:412,230:0.358:642:200,110:212,120 |
In the sample column above, AD=412,230 and DP=642 reveal 412 reference reads and 230 mutant reads, yielding a VAF of 35.8% (0.358) at the classical EGFR L858R hotspot locus.
Common Pitfalls & Troubleshooting Checklist
- Java Heap Exhaustion (GATK / Picard): If processes terminate with
java.lang.OutOfMemoryError, set your environment variable explicitly:export _JAVA_OPTIONS="-Xmx16g -Djava.io.tmpdir=/tmp". - Index Timestamp Mismatch: If tools complain that a BAM index is missing or out of date, regenerate with
samtools index -@ 4 sample.bam. - Conda Channel Conflicts: Avoid ABI incompatibilities between bioconda and conda-forge by setting strict priority:
conda config --set channel_priority strict. - Chromosome Naming Contamination: Always verify whether reference FASTA headers use
chr1or1; divergent naming causes silent alignment drops.