Skip to content

PairSeq mismatch: .gz + index=True ignores key_func` #108

Description

@NathaneilKnight

Dear developers,
Thanks for creating such a amazing tool!

Summary

When using .gz FASTQ with index=True, readSeqFile ignores key_func and uses rec.id directly as keys.
This causes PairSeq to fail for --coord presto because rec.id contains metadata (e.g. |CONSCOUNT=...) and does not match getCoordKey(...) results.

Steps to Reproduce

  1. Use PairSeq.py on presto Consensus files in .fastq.gz:
PairSeq.py -1 /path/to/*.R1_consensus-pass.fastq.gz \
          -2 /path/to/*.R2_consensus-pass.fastq.gz \
          --coord presto --outdir /path/to/out --outname sample
  1. Observe that no pairs are found (PASS=0).

Expected Behavior

--coord presto should match reads by stripping annotations after |, so paired reads should be found.

Actual Behavior

No pairs found (PASS=0), because .gz indexing uses rec.id directly.

Example:

rec.id = GTAGCATATG|CONSCOUNT=1|PRCONS=IgA_5-primer_Leading_sequence|PRFREQ=1.0
getCoordKey(..., coord_type="presto") -> GTAGCATATG

Mismatch caused by .gz indexing ignoring key_func.


Root Cause (in code)

Current .gz path (index=True):

seq_records = {rec.id: rec for rec in SeqIO.parse(handle, seq_type)}

This bypasses key_func.


Proposed Fix

Use key_func for .gz indexing as well:

if key_func is None:
    seq_records = {rec.id: rec for rec in SeqIO.parse(handle, seq_type)}
else:
    seq_records = {key_func(rec.id): rec for rec in SeqIO.parse(handle, seq_type)}

This makes .gz behavior consistent with non-.gz indexing.

Best regards,
Chenyu Pei

Metadata

Metadata

Assignees

Labels

bugSomething isn't working

Type

No type

Fields

No fields configured for issues without a type.

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions