Handle paired-end FASTQ files (R1/R2) using Biopython while keeping mates synchronized. Use when working with Illumina paired reads, synchronizing pairs, filtering both mates together with orphan routing, interleaving/deinterleaving, or matching mates by read name.
Handle paired-end FASTQ files (R1/R2) using Biopython while keeping mates synchronized. Use when working with Illumina paired reads, synchronizing pairs, filtering both mates together with orphan routing, interleaving/deinterleaving, or matching mates by read name.
tool_type
python
primary_tool
Bio.SeqIO
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).
The text before the space is IDENTICAL for both mates; the 1:/2: lives only after the space
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.
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.
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
defiterate_pairs(r1_file, r2_file, format='fastq'):
r1_records = SeqIO.parse(r1_file, format)
r2_records = SeqIO.parse(r2_file, format)
for r1, r2 inzip(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+):
deffilter_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}
withopen(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 inzip(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'] += 1elif r1_ok:
SeqIO.write(r1, o1, 'fastq')
counts['r1_orphan'] += 1elif r2_ok:
SeqIO.write(r2, o2, 'fastq')
counts['r2_orphan'] += 1return 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.
Create Interleaved File
Goal: Merge synchronized R1/R2 files into one interleaved file.
Approach: Zip both iterators and yield alternating records through a generator so nothing is materialized in memory.
Goal: Recover separate R1/R2 files from an interleaved file, streaming to avoid loading everything.
Approach: Parse once; route even-indexed records to R1, odd-indexed to R2.
Reference (BioPython 1.83+):
defdeinterleave_streaming(interleaved_file, r1_file, r2_file, format='fastq'):
records = SeqIO.parse(interleaved_file, format)
pairs = 0withopen(r1_file, 'w') as r1_h, open(r2_file, 'w') as r2_h:
for i, record inenumerate(records):
if i % 2 == 0:
SeqIO.write(record, r1_h, format)
else:
SeqIO.write(record, r2_h, format)
pairs += 1return pairs
Even/odd splitting only stays correct if the interleaved file has perfect alternation. If an upstream per-read filter left an orphan in the file, every record after it lands in the wrong output. Guard by checking mate_key equality between each even/odd pair after splitting, or deinterleave with a name check.
Paired Statistics
Count and Verify Pairs
defpaired_stats(r1_file, r2_file):
r1_count = sum(1for _ in SeqIO.parse(r1_file, 'fastq'))
r2_count = sum(1for _ in SeqIO.parse(r2_file, 'fastq'))
if r1_count != r2_count:
print(f'WARNING: Unequal counts! R1={r1_count}, R2={r2_count} -> files are desynced')
else:
print(f'Pairs: {r1_count}, total reads: {r1_count * 2}')
return r1_count, r2_count
Paired Quality Summary
defpaired_quality_summary(r1_file, r2_file):
r1_quals, r2_quals = [], []
for r1, r2 inzip(SeqIO.parse(r1_file, 'fastq'), SeqIO.parse(r2_file, 'fastq')):
r1_quals.append(sum(r1.letter_annotations['phred_quality']) / len(r1.seq))
r2_quals.append(sum(r2.letter_annotations['phred_quality']) / len(r2.seq))
print(f'R1 mean quality: {sum(r1_quals)/len(r1_quals):.1f}')
print(f'R2 mean quality: {sum(r2_quals)/len(r2_quals):.1f}')
R2 commonly shows lower mean quality than R1 (the reverse read is sequenced later in the run); a modest R1/R2 gap is expected, not a defect.
Find Paired Files
Auto-Detect R2 from R1
from pathlib import Path
deffind_r2(r1_path):
r1_path = Path(r1_path)
name = r1_path.name
patterns = [('_R1', '_R2'), ('_R1_', '_R2_'), ('.R1.', '.R2.'), ('_1', '_2')]
for p1, p2 in patterns:
if p1 in name:
r2_path = r1_path.parent / name.replace(p1, p2, 1)
if r2_path.exists():
return r2_path
returnNone
Order matters: test the specific _R1 patterns before the bare _1, otherwise sample_R1.fastq would match _1 and produce sample_R2.fastq only by luck. replace(..., 1) replaces the first occurrence only, so a sample named sample_R1_lane_R1.fastq swaps just the first token.
Compressed Paired Files
import gzip
defiterate_gzipped_pairs(r1_gz, r2_gz):
with gzip.open(r1_gz, 'rt') as r1_h, gzip.open(r2_gz, 'rt') as r2_h:
for r1, r2 inzip(SeqIO.parse(r1_h, 'fastq'), SeqIO.parse(r2_h, 'fastq')):
yield r1, r2
Use text mode 'rt', not 'rb', when handing a gzip handle to SeqIO.parse; the parser expects decoded text.
Common Errors
Symptom
Cause
Fix
Aligner reports "mismatched read names" or "unpaired reads"
R1 and R2 desynced by filtering/trimming one mate independently
Always filter both mates together; route orphans to separate files (see synchronized filter)
Silent mismapping, nonsensical insert sizes, no error
A per-mate operation dropped reads from one file -> i-th records no longer mates
Re-pair from source; never trust outputs from independent per-mate filtering
Mates never recognized as pairs
Mixing pre-1.8 /1/2 data with 1.8+ space-format ids, or comparing full descriptions instead of the pre-space ID
Match on mate_key (ID up to first space, /1//2 stripped), not the whole header
Deinterleave produces shifted/wrong pairs
An orphan in the interleaved file broke the strict R1,R2 alternation
Remove orphans before interleaving; verify each even/odd pair with mate_key after splitting
zip loop processes fewer reads than expected
R1 and R2 have unequal counts; zip stops at the shorter and silently drops the tail
Run paired_stats first; counts must be equal
Memory error on large files
list(SeqIO.parse(...)) materializes every record
Stream with generators; for random access use SeqIO.index/index_db
Related Skills
read-sequences - Parse individual FASTQ files and choose parse vs index
fastq-quality - Phred encoding and quality interpretation before paired filtering
filter-sequences - Single-file filtering criteria (apply to both mates here)
compressed-files - gzip vs BGZF handling for paired files
