Sequencing data quality assessment

Practical 8: Genomics, systems biology and multi-omics

Learning objectives

Decode Phred+33 qualities; calculate per-read expected errors; distinguish read, library and biological replication; apply an explicit QC rule.

Biological question and data

How do low-quality positions influence a small read-quality summary?

The fictional FASTQ records use ASCII Phred+33 encoding. The error probabilities are interpreted as calibrated per-base probabilities solely for this example. The practical does not assess mapping or a real sequencing run.

No biological download is required. Numerical inputs are hypothetical, not observed animal data. The runnable examples use base R. Count and abundance units are specified where used.

1. Decode quality strings

Confirm the sequence and quality strings have matching lengths.

reads <- c(read1="ACGTACGT",read2="ACGTACGT",read3="ACGTNNGT")
quality <- c(read1="IIIIIIII",read2="IIII!!II",read3="55555555")
stopifnot(identical(names(reads),names(quality)),all(nchar(reads)==nchar(quality)))
Q <- t(vapply(quality,function(q) utf8ToInt(q)-33,numeric(8)))
Q

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
read1   40   40   40   40   40   40   40   40
read2   40   40   40   40    0    0   40   40
read3   20   20   20   20   20   20   20   20

stopifnot(Q[1,1]==40,Q[2,5]==0,Q[3,1]==20)

I encodes Q40, ! Q0 and 5 Q20 under the stated convention. The convention must not be guessed from appearance.

2. Summarise errors and unknown calls

Sum encoded error probabilities and count unknown bases separately.

expected_errors <- rowSums(10^(-Q/10))
unknown <- vapply(reads,function(s) sum(strsplit(s,"")[[1]]=="N"),integer(1))
qc <- data.frame(read=names(reads),mean_Q=rowMeans(Q),expected_errors,unknown)
qc

       read mean_Q expected_errors unknown
read1 read1     40          0.0008       0
read2 read2     30          2.0006       0
read3 read3     20          0.0800       2

stopifnot(abs(expected_errors[1]-.0008)<1e-12,expected_errors[2]>2)

Two Q0 positions dominate the second read’s expected-error total. Mean quality alone compresses their location and severity.

3. Apply and report a specified rule

The supplied rule retains reads with expected errors below .1 and no unknown calls; it is illustrative, not a universal threshold.

qc$keep <- qc$expected_errors<.1 & qc$unknown==0
qc

       read mean_Q expected_errors unknown  keep
read1 read1     40          0.0008       0  TRUE
read2 read2     30          2.0006       0 FALSE
read3 read3     20          0.0800       2 FALSE

c(input_reads=nrow(qc),retained=sum(qc$keep))

input_reads    retained 
          3           1 

stopifnot(sum(qc$keep)==1)

The third read meets the encoded-error rule but fails the unknown-base rule. Report the exclusions and their reason.

4. Compare average quality with expected errors

The Phred transformation is nonlinear. Transforming mean quality to an error probability is different from averaging per-base error probabilities.

quality <- c(10,40)
c(error_from_mean_Q=10^(-mean(quality)/10),
 mean_error_probability=mean(10^(-quality/10)),
 expected_errors=sum(10^(-quality/10)))

     error_from_mean_Q mean_error_probability        expected_errors 
           0.003162278            0.050050000            0.100100000 

The low-quality base contributes most of the expected errors. Expected error sums use linearity of expectation and do not require independent errors; the probability of no errors would require additional dependence assumptions.

5. Evaluate consequences of filtering

A stricter threshold can improve retained-read summaries while discarding evidence unevenly between samples. Report retained counts and proportions per library, not only the pooled fraction. Quality filtering is one part of a workflow: mapping, duplicate evidence, contamination and biological sample integrity still need checks.

The small quality strings demonstrate encoding and decisions. They do not characterize an instrument run or establish the best threshold for real sequencing data. Choose a rule based on the intended downstream use and inspect its effect on group balance.

Interpretation challenge

Two libraries retain 90% of reads overall, but one treatment group loses twice as many reads. What should be inspected?

Suggested reasoning

Group-specific quality and processing, library sizes, acquisition/batch structure and downstream exposure. The pooled retention fraction can hide an imbalanced measurement process.

Interpretation and limits

Quality encoding, adapter content, length distributions, contaminants, mapping, coverage and sample identity matter in a complete sequencing workflow. A threshold chosen after looking at a desired result can bias downstream analysis. More retained reads do not create more biological replicates.

Exercises

1. Why can mean Q hide a problematic read?
2. Does summing error probabilities require independent errors?
3. What information would mapping QC add?

Suggested answers

1. A few extremely low-quality positions can be diluted in the average.
2. No; expectation is additive. Independence matters for some probability calculations, not this sum.
3. Reference alignment, multimapping, coverage and mapping-rate patterns, each requiring interpretation.

Continue learning

B05 Genomics and B06 Transcriptomics.

Continue with RNA-seq. Return to the practical index or the unified map.