Bio Paired End Fastq

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.

Paired-End FASTQ

"Work with my paired-end FASTQ files" -> Iterate R1/R2 pairs in sync, filter both mates together (routing orphans out), interleave/deinterleave files, and match mates by read name.

• Python: SeqIO.parse() with zip() iteration (BioPython)

The Governing Principle: R1 and R2 Are Parallel Streams

Aligners (bwa mem, bowtie2, STAR) consume R1 and R2 as two parallel streams and pair the i-th record of each file: same order, same count. They assume the k-th read in R1 is the mate of the k-th read in R2.

This makes independent per-mate processing the #1 paired-end correctness trap. Filtering or trimming ONE mate without the other DESYNCS the files:

• Best case: the aligner detects a read-name mismatch and crashes loudly.
• Worst case: it silently pairs the wrong R1 with the wrong R2 -> mismapping, corrupt insert sizes, no error at all.

Governing rule: never filter, trim, sort, or subsample one mate independently. Process both mates as a unit. When a read fails but its mate passes, route the survivor to a separate singleton/orphan file rather than leaving a gap that desyncs the stream. Proper paired trimmers (Trimmomatic PE, fastp, cutadapt -p) do exactly this: synchronized paired output plus separate orphan files.

Required Import

from Bio import SeqIO

Read-Name Conventions: How Mates Are Matched

Two distinct naming layers exist. File naming (which file is R1 vs R2) is separate from read naming (how a tool decides two records are mates).

File naming patterns

• sample_R1.fastq / sample_R2.fastq
• sample_1.fastq / sample_2.fastq
• sample_R1_001.fastq / sample_R2_001.fastq (Illumina bcl2fastq)
• sample.R1.fastq.gz / sample.R2.fastq.gz

Read naming: mate matched by shared ID up to the first whitespace

In CASAVA 1.8+ the second field is <read>:<is_filtered>:<control>:<index> -> read=1 or 2 (mate number), is_filtered=Y (failed chastity) or N, control=0 normally, index=barcode.

Most tools (and Biopython) take the read ID as everything up to the first whitespace. Biopython puts that token in record.id and the full header line in record.description. So for 1.8+ data, r1.id == r2.id directly; the 1:/2: distinction is only visible in record.description. For pre-1.8 data, strip the /1//2 suffix before comparing.

A normalizer that handles both eras:

def mate_key(record):
    return record.id.rsplit('/', 1)[0]

record.id already excludes anything after the first space, so this single rsplit covers both the space-format (1.8+) and the slash-suffix (pre-1.8) conventions.

Iterate Pairs Together

Basic Paired Iteration

r1_records = SeqIO.parse('reads_R1.fastq', 'fastq')
r2_records = SeqIO.parse('reads_R2.fastq', 'fastq')

for r1, r2 in zip(r1_records, r2_records):
    print(f'R1: {r1.id}, R2: {r2.id}')
    print(f'Lengths: {len(r1.seq)}, {len(r2.seq)}')

zip stops at the shorter iterator. If R1 has 1000 reads and R2 has 998, zip silently processes 998 and drops the tail with no warning. Verify counts match (see Paired Statistics) before trusting a zip loop on files of unknown provenance.

Verify Pair Matching

def iterate_pairs(r1_file, r2_file, format='fastq'):
    r1_records = SeqIO.parse(r1_file, format)
    r2_records = SeqIO.parse(r2_file, format)

    for r1, r2 in zip(r1_records, r2_records):
        if mate_key(r1) != mate_key(r2):
            raise ValueError(f'Pair mismatch: {r1.id} vs {r2.id}')
        yield r1, r2

for r1, r2 in iterate_pairs('reads_R1.fastq', 'reads_R2.fastq'):
    process_pair(r1, r2)

Filter Pairs Together (Synchronized, With Orphan Routing)

Goal: Quality-filter paired reads so that R1 and R2 stay in lockstep, and reads whose mate was discarded are routed to orphan files instead of silently desyncing the stream.

Approach: Stream both files together with zip. Evaluate both mates. If both pass, write to the paired outputs. If exactly one passes, write the survivor to its orphan file. This mirrors the four-output behavior of Trimmomatic PE (paired R1, paired R2, orphan R1, orphan R2).

Reference (BioPython 1.83+):

def filter_pairs_synced(r1_in, r2_in, r1_out, r2_out, r1_orphan, r2_orphan, min_qual=25):
    '''Keep a pair only if both mates pass; route lone survivors to orphan files.'''
    r1_records = SeqIO.parse(r1_in, 'fastq')
    r2_records = SeqIO.parse(r2_in, 'fastq')

    counts = {'paired': 0, 'r1_orphan': 0, 'r2_orphan': 0}
    with open(r1_out, 'w') as p1, open(r2_out, 'w') as p2, \
         open(r1_orphan, 'w') as o1, open(r2_orphan, 'w') as o2:
        for r1, r2 in zip(r1_records, r2_records):
            r1_ok = sum(r1.letter_annotations['phred_quality']) / len(r1.seq) >= min_qual
            r2_ok = sum(r2.letter_annotations['phred_quality']) / len(r2.seq) >= min_qual
            if r1_ok and r2_ok:
                SeqIO.write(r1, p1, 'fastq')
                SeqIO.write(r2, p2, 'fastq')
                counts['paired'] += 1
            elif r1_ok:
                SeqIO.write(r1, o1, 'fastq')
                counts['r1_orphan'] += 1
            elif r2_ok:
                SeqIO.write(r2, o2, 'fastq')
                counts['r2_orphan'] += 1
    return counts

The paired outputs stay synchronized because every pass-both write goes to BOTH files in the same iteration. Mean quality is one criterion; swap the test for a length threshold, adapter check, or any predicate, but always apply it to both mates and route orphans the same way. min_qual=25 is a common Q-score cutoff (Phred 25 ~= 0.3% error); tune per experiment.

Interleave Pairs

Interleaved FASTQ holds both mates in one file alternating R1, R2, R1, R2 (record 2k = forward, 2k+1 = reverse). bwa mem -p reads this format. The strict alternation IS the pairing, so a single mate-less read shifts every downstream record by one and desyncs everything. Only interleave files that are known to be synchronized, and pull orphans out first.