Add genomic coordinates to export_end_sequences and write_fasta headers (2.1.6)

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)

VERSION.txt

@@ -1 +1 @@
-2.1.5
+2.1.6

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

src/isotools/_transcriptome_io.py

@@ -2093,23 +2093,27 @@ 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(
                         "Skipping transcript %s of gene %s at %s:%d-%d: "
                         "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(

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}"