bio-fastq-quality

## Version Compatibility

Reference examples tested with: BioPython 1.83+

Before using code patterns, verify installed versions match. If versions differ:
- Python: `pip show <package>` then `help(module.function)` to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

# FASTQ Quality Scores

**"Filter my FASTQ reads by quality score"** → Access, analyze, and filter Phred quality scores, trim low-quality bases, and generate per-position quality profiles.
- Python: `SeqIO.parse()` with `letter_annotations['phred_quality']` (BioPython)

Analyze and manipulate FASTQ quality scores using Biopython.

## Required Imports

```python
from Bio import SeqIO
from Bio.Seq import Seq
```

## Accessing Quality Scores

Quality scores are stored in `letter_annotations['phred_quality']` as a list of integers.

```python
for record in SeqIO.parse('reads.fastq', 'fastq'):
    qualities = record.letter_annotations['phred_quality']
    print(record.id, qualities[:10])  # First 10 quality scores
```

## Quality Score Basics

| Phred Score | Error Probability | Accuracy |
|-------------|-------------------|----------|
| 10 | 1 in 10 | 90% |
| 20 | 1 in 100 | 99% |
| 30 | 1 in 1000 | 99.9% |
| 40 | 1 in 10000 | 99.99% |

## Code Patterns

### Calculate Average Quality per Read
```python
for record in SeqIO.parse('reads.fastq', 'fastq'):
    quals = record.letter_annotations['phred_quality']
    avg_qual = sum(quals) / len(quals)
    print(f'{record.id}: {avg_qual:.1f}')
```

### Filter Reads by Mean Quality
```python
def high_quality_reads(records, min_avg_qual=20):
    for record in records:
        quals = record.letter_annotations['phred_quality']
        if sum(quals) / len(quals) >= min_avg_qual:
            yield record

records = SeqIO.parse('reads.fastq', 'fastq')
good_reads = high_quality_reads(records, min_avg_qual=25)
SeqIO.write(good_reads, 'filtered.fastq', 'fastq')
```

### Filter by Minimum Quality at Any Position
```python
def all_bases_above(records, min_qual=20):
    for record in records:
        if min(record.letter_annotations['phred_quality']) >= min_qual:
            yield record
```

### Trim Low-Quality Ends (3' Trimming)
```python
def trim_low_quality(record, min_qual=20):
    quals = record.letter_annotations['phred_quality']
    trim_pos = len(quals)
    for i in range(len(quals) - 1, -1, -1):
        if quals[i] >= min_qual:
            trim_pos = i + 1
            break
    return record[:trim_pos]

records = SeqIO.parse('reads.fastq', 'fastq')
trimmed = (trim_low_quality(rec) for rec in records)
SeqIO.write(trimmed, 'trimmed.fastq', 'fastq')
```

### Sliding Window Quality Trim

**Goal:** Trim a read at the first position where average quality in a sliding window drops below threshold.

**Approach:** Slide a fixed-size window across quality scores; when the window average falls below the cutoff, truncate the record at that position.

**Reference (BioPython 1.83+):**
```python
def sliding_window_trim(record, window_size=5, min_avg_qual=20):
    quals = record.letter_annotations['phred_quality']
    for i in range(len(quals) - window_size + 1):
        window = quals[i:i + window_size]
        if sum(window) / window_size < min_avg_qual:
            return record[:i] if i > 0 else None
    return record
```

### Quality Statistics Summary
```python
import statistics

all_quals = []
for record in SeqIO.parse('reads.fastq', 'fastq'):
    all_quals.extend(record.letter_annotations['phred_quality'])

print(f'Mean quality: {statistics.mean(all_quals):.1f}')
print(f'Median quality: {statistics.median(all_quals):.1f}')
print(f'Min quality: {min(all_quals)}')
print(f'Max quality: {max(all_quals)}')
```

### Per-Position Quality Profile

**Goal:** Compute mean quality at each read position to identify systematic quality drops (e.g., read-end degradation).

**Approach:** Accumulate quality scores by position across all reads, then compute per-position means.

**Reference (BioPython 1.83+):**
```python
from collections import defaultdict

position_quals = defaultdict(list)
for record in SeqIO.parse('reads.fastq', 'fastq'):
    for i, q in enumerate(record.letter_annotations['phred_quality']):
        position_quals[i].append(q)

for pos in sorted(position_quals.keys())[:20]:
    quals = position_quals[pos]
    print(f'Position {pos}: mean={sum(quals)/len(quals):.1f}')
```

### Count Reads by Quality Threshold
```python
thresholds = [20, 25, 30, 35]
counts = {t: 0 for t in thresholds}

for record in SeqIO.parse('reads.fastq', 'fastq'):
    avg = sum(record.letter_annotations['phred_quality']) / len(record.seq)
    for t in thresholds:
        if avg >= t:
            counts[t] += 1

for t, c in counts.items():
    print(f'Q>={t}: {c} reads')
```

### Remove N Bases and Low Quality Together
```python
def clean_read(record, min_qual=20):
    seq = str(record.seq)
    quals = record.letter_annotations['phred_quality']
    keep = [(s, q) for s, q in zip(seq, quals) if s != 'N' and q >= min_qual]
    if not keep:
        return None
    new_seq, new_quals = zip(*keep)
    new_record = record[:0]  # Empty copy with same metadata
    new_record.seq = Seq(''.join(new_seq))
    new_record.letter_annotations['phred_quality'] = list(new_quals)
    return new_record
```

## FASTQ Format Variants

| Variant | Format String | Quality Encoding | ASCII Range |
|---------|---------------|------------------|-------------|
| Sanger/Illumina 1.8+ | `'fastq'` | Phred+33 (standard) | 33-126 |
| Solexa | `'fastq-solexa'` | Solexa+64 | 59-126 |
| Illumina 1.3-1.7 | `'fastq-illumina'` | Phred+64 | 64-126 |

Most modern data uses standard `'fastq'` (Sanger/Illumina 1.8+).

## Quality Score Conversion

For legacy data using different quality encodings:

```python
from Bio.SeqIO.QualityIO import phred_quality_from_solexa, solexa_quality_from_phred
```

### Convert Solexa to Phred

```python
from Bio.SeqIO.QualityIO import phred_quality_from_solexa

# Convert single score
solexa_score = 10
phred_score = phred_quality_from_solexa(solexa_score)

# Convert list of scores
solexa_scores = [10, 20, 30, 40]
phred_scores = [phred_quality_from_solexa(s) for s in solexa_scores]
```

### Convert Phred to Solexa

```python
from Bio.SeqIO.QualityIO import solexa_quality_from_phred

phred_score = 30
solexa_score = solexa_quality_from_phred(phred_score)
```

### Convert Between FASTQ Variants

```python
from Bio import SeqIO

# Read old Illumina format, write standard format
records = SeqIO.parse('old_reads.fastq', 'fastq-illumina')
SeqIO.write(records, 'standard_reads.fastq', 'fastq')

# Read Solexa format, write standard format
records = SeqIO.parse('solexa_reads.fastq', 'fastq-solexa')
SeqIO.write(records, 'standard_reads.fastq', 'fastq')
```

### Auto-Detect Quality Encoding

**Goal:** Determine which FASTQ quality encoding (Sanger, Solexa, or Illumina 1.3+) a file uses.

**Approach:** Sample quality lines, find the minimum ASCII value, and compare against known offset ranges.

**Reference (BioPython 1.83+):**
```python
def detect_quality_encoding(filepath, sample_size=1000):
    min_qual = 126
    max_qual = 0
    count = 0

    with open(filepath) as f:
        for i, line in enumerate(f):
            if i % 4 == 3:  # Quality line
                for char in line.strip():
                    min_qual = min(min_qual, ord(char))
                    max_qual = max(max_qual, ord(char))
                count += 1
                if count >= sample_size:
                    break

    if min_qual < 59:
        return 'fastq'  # Sanger/Illumina 1.8+ (Phred+33)
    elif min_qual < 64:
        return 'fastq-solexa'  # Solexa+64
    else:
        return 'fastq-illumina'  # Illumina 1.3+ (Phred+64)
```

## Common Errors

| Error | Cause | Solution |
|-------|-------|----------|
| `KeyError: 'phred_quality'` | Not FASTQ or wrong variant | Check format, try 'fastq-illumina' |
| Quality scores all 0 | Wrong encoding assumed | Try different FASTQ variant |
| Trimmed reads empty | Too aggressive trimming | Lower quality threshold |

## Related Skills

- read-sequences - Parse FASTQ files
- filter-sequences - Filter reads by other criteria (length, content)
- paired-end-fastq - Handle R1/R2 paired quality filtering
- sequence-statistics - Generate summary statistics including quality
- alignment-files - After filtering, align reads with bwa/bowtie2; quality scores in BAM