Skip to content

Don't use random FASTQ access during hcid check #115

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
BioWilko merged 2 commits into from
May 18, 2026
Merged

Don't use random FASTQ access during hcid check #115

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c4a187b
Don't use random FASTQ access during hcid check
BioWilko May 18, 2026
8ae1a42
stop iteration early if all needed reads already found
BioWilko May 18, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
25 changes: 20 additions & 5 deletions bin/check_hcid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -189,17 +189,32 @@ def report_findings(hcid_dict, read_file, prefix):
w.writerow({key: hcid_dict[taxid][key] for key in keys})

found = []
records = pyfastx.Fastq(read_file, full_name=False, build_index=True)
needed = {
name
for taxid in hcid_dict
if hcid_dict[taxid]["mapped_found"]
for name in hcid_dict[taxid]["mapped_read_ids"]
}
read_records = {}
for name, seq, qual in pyfastx.Fastq(read_file, full_name=False, build_index=False):

@Sam-Sims Sam-Sims May 18, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

very minor but could consider an early return here once all needed reads are found e.g all hcid reads were in first 10% of a 5gb fastq file then no need to continue parsing

as I type this though probably not worth it in our use case - so not blocking

@BioWilko BioWilko May 18, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah that's a good idea, let me have a look

if name in needed:
read_records[name] = (seq, qual)

for taxid in hcid_dict:
if hcid_dict[taxid]["mapped_found"]:
quals = []
lens = []
with open("%s.reads.fq" % taxid, "w") as f_reads:
for mapped_name in hcid_dict[taxid]["mapped_read_ids"]:
record = records[mapped_name]
f_reads.write(record.raw)
quals.append(stats.fmean(record.quali) if record.quali else 0)
lens.append(len(record.seq))
record = read_records.get(mapped_name)
if record is None:
print(f"WARNING: read {mapped_name} not found in FASTQ, skipping", file=sys.stderr)
continue
seq, qual = record
f_reads.write(f"@{mapped_name}\n{seq}\n+\n{qual}\n")
phred = [ord(c) - 33 for c in qual]
quals.append(stats.fmean(phred) if phred else 0)
lens.append(len(seq))
with open("%s.warning.json" % taxid, "w") as f_warn:
msg1 = f"WARNING: Found {hcid_dict[taxid]['classified_count']} classified reads ({hcid_dict[taxid]['mapped_count']} mapped reads) of {hcid_dict[taxid]['name']} and {hcid_dict[taxid]['classified_parent_count']} classified reads for the parent taxon.\n"
msg2 = f"Mapping details for required references (ref_accession:mapped_read_count:fraction_ref_covered) {hcid_dict[taxid]['mapped_required_details']}.\n"
Expand Down
2 changes: 1 addition & 1 deletion nextflow.config
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ manifest {
description = 'Classify metagenomic sequence data from human respiratory infections.'
mainScript = 'main.nf'
nextflowVersion = '>=20.10.0'
version = 'v2.2.0'
version = 'v2.2.1'
}

profiles {
Expand Down
Loading