diff --git a/bin/check_hcid.py b/bin/check_hcid.py index 3e34a71a..0f0c2b7b 100755 --- a/bin/check_hcid.py +++ b/bin/check_hcid.py @@ -189,17 +189,39 @@ 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): + if name in needed: + read_records[name] = (seq, qual) + + if needed - set(read_records.keys()) == set(): + break + 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" diff --git a/nextflow.config b/nextflow.config index af207bfc..e6a33b27 100644 --- a/nextflow.config +++ b/nextflow.config @@ -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 {