CalcTrueQuality

Description

CalcTrueQuality analyzes mapped reads to calculate the actual accuracy of reported quality scores, then generates recalibration matrices to improve quality score accuracy. The tool performs statistical analysis on alignment data to determine how often bases with specific quality scores are actually correct or incorrect under various contextual conditions.

For diploid organisms with high heterozygosity rates, the tool may produce inaccurate recalibration at high quality scores unless SNP locations are masked or variations are called. For example, recalibrating human reads mapped to an unmasked human reference would generate a maximum Q-score of roughly 30 due to the human 1/1000 SNP rate.

If you have multiple sam/bam files demultiplexed from a single sequencing run, it is recommended to use all of them as input for increased statistical power. Once matrices are generated, recalibration can be done on mapped or unmapped reads. Better results may be obtained by recalibrating the fastq files and remapping the calibrated reads.

Basic Usage

calctruequality.sh in=<sam,sam,...sam> path=<directory>

Two-Step Process

Step 1: Generate matrices using CalcTrueQuality

calctruequality.sh in=mapped.sam path=calibration_dir

Step 2: Recalibrate reads using BBDuk

bbduk.sh in=reads.fq out=calibrated.fq recalibrate

Parameters

Parameters are organized by their function in the quality calibration process. Each parameter group serves a specific role in matrix generation and quality score recalibration.

Examples

calctruequality.sh in=aligned.bam path=calibration_data

calctruequality.sh in=sample1.bam,sample2.bam,sample3.bam path=merged_calibration

calctruequality.sh in=human.bam ref=human_genome.fa callvars=t ploidy=2 path=human_calibration

calctruequality.sh in=reads.bam passes=1 trackall=t path=quick_calibration

calctruequality.sh in=sample.bam vcf=known_variants.vcf path=variant_aware_calibration

Algorithm Details

Quality Score Recalibration Strategy

CalcTrueQuality implements a multi-dimensional matrix-based approach using separate "good" and "bad" count matrices (GBMatrixSet class) to track correct and incorrect base calls. The modify() method combines observed error rates with expected error rates using QualityTools.PROB_ERROR lookup tables and observation cutoff weighting to ensure statistical reliability when observation counts are low.

Matrix-Based Recalibration

The algorithm maintains parallel good/bad matrices with specific dimensional structures for different contextual factors. Each matrix type uses multi-dimensional long arrays with predefined size limits (QMAX2=42 for quality, BMAX=6 for bases, LENMAX=500 for position, TMAX=2048 for tiles):

• Base Context Matrices: qbpGoodMatrix[2][42][6][500] and qb123GoodMatrix[2][42][6][6][6] track accuracy using baseToNum[] lookup table conversion (A=0, C=1, G=2, T=3, U=3, E=4, other=5)
• Quality Context Matrices: q102GoodMatrix[2][42][42][42] considers quality scores of neighboring bases using Tools.mid(QMAX, quals[position], 0) boundary clamping
• Position Matrices: qapGoodMatrix[2][42][43][500] accounts for position-within-read effects using Tools.min(qpos, LENMAX-1) position capping
• Tile Matrices: qptGoodMatrix[2][42][500][2048] uses parseTile() method to extract tile numbers from Illumina read names via colon-separated field parsing

Two-Pass Recalibration

The two-pass system uses initializeMatrices() to load CountMatrixSet objects for each pass, then applies recalibrate() with pass-specific boolean flags (pass0, pass1) to iteratively refine quality scores. Pass 1 matrices are used by CountMatrixSet.recalibrate() to generate improved quality scores, which are then used as input for pass 2 matrix generation through the same GBMatrixSet accumulation process.

Error Probability Estimation

For each base position, the CountMatrixSet class implements four distinct error probability estimation methods selectable via USE_* boolean flags:

• Average Method (estimateErrorProbAvg): Sums error probabilities from all loaded matrix types and returns sum/(float)x where x is the number of contributing matrices
• Maximum Method (estimateErrorProbMax): Uses Tools.max(max, f) to track the highest error probability from any loaded matrix, returning the maximum value found
• Weighted Average (estimateErrorProbWeighted): Combines observation counts from all matrices using (bad+fakeBad)/(sum+fakeSum) where fakeBad and fakeSum provide regularization based on obs_cutoff and BAD_CUTOFF constants
• SMR Method (estimateErrorProbSMR): Implements a square-mean-root calculation combining matrix probabilities with geometric averaging

Variant Handling

The algorithm integrates with CallVariants class to generate VarMap and ScafMap objects. The AnalyzeVars.fixVars() method processes reads to identify variant positions, converting 'S' (substitution) symbols to 'V' (variant) in the match array when positions are found in the loaded variant map. This prevents true biological variants from being counted as sequencing errors during matrix accumulation.

Statistical Confidence

The modify() method implements Bayesian regularization using observation cutoffs (OBSERVATION_CUTOFF[pass]) to prevent overfitting. For low-observation entries, it computes sum2=sum+cutoff and bad2=bad+expected*cutoff where expected comes from QualityTools.PROB_ERROR[phred] lookup table. The final error probability is calculated as bad2/sum2, blending observed and expected error rates based on observation confidence.

Memory and Performance

The algorithm uses Worker threads (extending Thread class) with either ConcurrentReadInputStream or SamReadStreamer for parallel processing. Memory usage is determined by matrix dimensions: largest matrices like qb123GoodMatrix require 2×42×6×6×6×8 bytes (≈250KB each), while 5-dimensional matrices like q12b12GoodMatrix need 2×42×42×6×6×8 bytes (≈2.5MB each). The process_MT() method manages thread pools with Tools.mid(1, threads, 16) worker thread limits.

Support

For questions and support:

• Email: bbushnell@lbl.gov
• Documentation: bbmap.org