diff --git a/CHANGELOG.md b/CHANGELOG.md index 96912f8..4145a47 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,11 @@ * planned new feature: during import of long reads, (optionally) correct for short exon alignment issues. * separate new read import and classification of isoforms. +## [2.1.6] + +* new: `export_end_sequences` now includes strand in the fasta header (`chrom:start-end:strand`) +* new: `write_fasta` gained an `add_coord` option to include the genomic location in the fasta header -- the transcript's genomic span for transcript sequences, or the coding sequence (annotated CDS, or predicted ORF) used for translation for protein sequences + ## [2.1.5] * fixed: `add_hmmer_domains` crashed with `TypeError: Function call with ambiguous argument types` since pyhmmer 0.7.0, since `get_hmmer_sequences` returned a plain list instead of a `DigitalSequenceBlock`; also crashed with `AssertionError: expression should be a string` for the documented `query=True`/`ref_query=True` default ("include all transcripts") (#53) diff --git a/VERSION.txt b/VERSION.txt index cd57a8b..399088b 100644 --- a/VERSION.txt +++ b/VERSION.txt @@ -1 +1 @@ -2.1.5 +2.1.6 diff --git a/requirements.txt b/requirements.txt index 8741f13..61ff48c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,7 +6,7 @@ # pip install . # pip freeze | grep -v '^isotools==' > requirements.txt CPAT==3.0.5 -biopython==1.87 +biopython==1.88 certifi==2026.7.22 charset-normalizer==3.4.9 contourpy==1.3.3 @@ -21,7 +21,7 @@ matplotlib==3.11.1 narwhals==2.24.0 numba==0.66.0 numpy==2.4.6 -packaging==26.2 +packaging==26.3 pandas==3.0.5 patsy==1.0.2 pillow==12.3.0 diff --git a/src/isotools/_transcriptome_io.py b/src/isotools/_transcriptome_io.py index d1ac86d..3b76ed5 100644 --- a/src/isotools/_transcriptome_io.py +++ b/src/isotools/_transcriptome_io.py @@ -2093,13 +2093,17 @@ def export_end_sequences( - 1 # exclusive end -> last included base ) window_here = window if is_plus else window[::-1] - pos = (gene.chrom, center - window_here[0], center + window_here[1] + 1) + chr, loc_start, loc_end = ( + gene.chrom, + center - window_here[0], + center + window_here[1] + 1, + ) - if unique_loc and pos in known_positions[gene.chrom]: + if unique_loc and (chr, loc_start, loc_end) in known_positions[chr]: n_skipped_dup += 1 continue - seq = ref.fetch(*pos) + seq = ref.fetch(chr, loc_start, loc_end) if len(seq) != expected_len: logger.debug( @@ -2107,9 +2111,9 @@ def export_end_sequences( "fetched sequence length (%d) does not match expected length (%d)", transcript_id, gene.id, - pos[0], - pos[1], - pos[2], + chr, + loc_start, + loc_end, len(seq), expected_len, ) @@ -2119,9 +2123,9 @@ def export_end_sequences( if not is_plus: seq = reverse_complement(seq) fh.write( - f">{gene.id}\t{transcript_id}\t{pos[0]}:{pos[1]}-{pos[2]}\n{seq}\n" + f">{gene.id}\t{transcript_id}\t{chr}:{loc_start}-{loc_end}:{transcript['strand']}\n{seq}\n" ) - known_positions[gene.chrom].add(pos) + known_positions[chr].add((chr, loc_start, loc_end)) n_written += 1 logger.info( @@ -2615,6 +2619,7 @@ def write_fasta( reference=False, protein=False, coverage=None, + add_coord=False, **filter_args, ): """ @@ -2625,6 +2630,9 @@ def write_fasta( :param protein: Return protein sequences (ORF) instead of transcript sequences. :param coverage: By default, the coverage is not added to the header of the fasta. If set, the allowed values are: 'all', or 'sample'. 'all' - total coverage for all samples; 'sample' - coverage by sample. + :param add_coord: If set, include the genomic location "chr:start-end:strand" in the header. For transcript + sequences, this is the transcript's genomic span; for protein sequences, this is the coding sequence + (annotated CDS, or predicted ORF if not annotated) used for translation. :param fn: The filename to write the fasta. :param gzip: Compress the output as gzip. :param filter_args: Additional filter arguments (e.g. "region", "gois", "query") are passed to iter_transcripts. @@ -2646,15 +2654,35 @@ def write_fasta( tr_seqs = gene.get_sequence( genome_fn, transcript_ids, reference=reference, protein=protein ) + transcripts = gene.ref_transcripts if reference else gene.transcripts if len(tr_seqs) > 0: - f.write( - "\n".join( - f">{gene.id}_{k} gene={gene.name}" - f'{(" coverage=" + (str(gene.coverage[:, k].sum()) if coverage == "all" else str(gene.coverage[:, k])) if coverage else "")}\n{v}' - for k, v in tr_seqs.items() + lines = [] + for k, v in tr_seqs.items(): + transcript = transcripts[k] + coord_info = "" + if add_coord: + if protein: + cds = transcript.get("CDS", transcript.get("ORF")) + if cds: + coord_info = ( + f" {gene.chrom}:{cds[0]}-{cds[1]}:{gene.strand}" + ) + else: + coord_info = f" {gene.chrom}:{transcript['exons'][0][0]}-{transcript['exons'][-1][1]}:{gene.strand}" + coverage_info = ( + " coverage=" + + ( + str(gene.coverage[:, k].sum()) + if coverage == "all" + else str(gene.coverage[:, k]) + ) + if coverage + else "" ) - + "\n" - ) + lines.append( + f">{gene.id}_{k}{coord_info} gene={gene.name}{coverage_info}\n{v}" + ) + f.write("\n".join(lines) + "\n") def export_alternative_splicing( diff --git a/tests/export_test.py b/tests/export_test.py new file mode 100644 index 0000000..9352939 --- /dev/null +++ b/tests/export_test.py @@ -0,0 +1,72 @@ +from pysam import FastaFile +from isotools.transcriptome import Transcriptome + + +def _example_transcriptome(): + isoseq = Transcriptome.from_reference("tests/data/example.gff.gz") + for sa in ("CTL", "VPA"): + isoseq.add_sample_from_bam( + f"tests/data/example_1_{sa}.bam", + sample_name=sa, + group=sa, + platform="SequelII", + ) + return isoseq + + +def test_write_fasta_add_coord(tmp_path): + isoseq = _example_transcriptome() + with FastaFile("tests/data/example.fa") as genome_fh: + for gene in isoseq: + gene.add_orfs(genome_fh, reference=True) + gene.add_orfs(genome_fh, reference=False) + + gene = isoseq["FN1"] + ref_transcript = gene.ref_transcripts[0] + assert ( + "CDS" in ref_transcript + ), "expected FN1 reference transcript 0 to have an annotated CDS" + tr_start, tr_end = ref_transcript["exons"][0][0], ref_transcript["exons"][-1][1] + cds_start, cds_end = ref_transcript["CDS"] + + tr_fn = tmp_path / "tr.fa" + prot_fn = tmp_path / "prot.fa" + isoseq.write_fasta( + "tests/data/example.fa", + str(tr_fn), + reference=True, + protein=False, + add_coord=True, + gois=["FN1"], + ) + isoseq.write_fasta( + "tests/data/example.fa", + str(prot_fn), + reference=True, + protein=True, + add_coord=True, + gois=["FN1"], + ) + + tr_header = tr_fn.read_text().splitlines()[0] + prot_header = prot_fn.read_text().splitlines()[0] + + assert ( + tr_header + == f">{gene.id}_0 {gene.chrom}:{tr_start}-{tr_end}:{gene.strand} gene={gene.name}" + ) + assert ( + prot_header + == f">{gene.id}_0 {gene.chrom}:{cds_start}-{cds_end}:{gene.strand} gene={gene.name}" + ) + + +def test_write_fasta_no_coord_by_default(tmp_path): + # regression test: add_coord defaults to False, so the header must be + # unchanged from before add_coord was introduced. + isoseq = _example_transcriptome() + fn = tmp_path / "tr.fa" + isoseq.write_fasta("tests/data/example.fa", str(fn), reference=True, gois=["FN1"]) + header = fn.read_text().splitlines()[0] + gene = isoseq["FN1"] + assert header == f">{gene.id}_0 gene={gene.name}"