|
| 1 | +import pandas as pd |
| 2 | +import numpy as np |
| 3 | +import pytest |
| 4 | + |
| 5 | +from thoraxe.subexons.tidy import get_tidy_table |
| 6 | + |
| 7 | + |
| 8 | +def _mk_row(gene="G", trx="T", subexon_rank=15, strand=-1, start=0, end=None, label="X"): |
| 9 | + # End defaults to one base before start on the negative strand or one base after on positive |
| 10 | + if end is None: |
| 11 | + end = start - 3 if strand == -1 else start + 3 |
| 12 | + |
| 13 | + # Keep protein length consistent with S_exon_Lengths ("21") |
| 14 | + subexon_prot = "A" * 21 |
| 15 | + subexon_nt = "ATG" * 7 # 21 aa -> 63 nt; not strictly required but tidy to keep |
| 16 | + |
| 17 | + # NOTE: get_tidy_table reads S_exon_Sequences (plural) from input rows |
| 18 | + s_exon_seq = "A" * 21 |
| 19 | + |
| 20 | + return { |
| 21 | + "GeneID": gene, |
| 22 | + "TranscriptIDCluster": trx, |
| 23 | + "SubexonRank": subexon_rank, |
| 24 | + "Strand": strand, |
| 25 | + "SubexonCodingStart": int(start), |
| 26 | + "SubexonCodingEnd": int(end), |
| 27 | + |
| 28 | + # Columns consumed by tidy.get_tidy_table() |
| 29 | + "SubexonProteinSequence": subexon_prot, |
| 30 | + "SubexonSequence": subexon_nt, |
| 31 | + |
| 32 | + # One s-exon per row in these tests; order is what we care about |
| 33 | + "S_exons": label, |
| 34 | + "S_exon_Lengths": "21", |
| 35 | + "S_exon_Sequences": s_exon_seq, |
| 36 | + } |
| 37 | + |
| 38 | + |
| 39 | +@pytest.mark.parametrize( |
| 40 | + "species,trx,strand,starts", |
| 41 | + [ |
| 42 | + # Platypus (ornithorhynchus_anatinus) transcripts |
| 43 | + ("ornithorhynchus_anatinus", "ENSOANT00000050045", -1, |
| 44 | + {"17_2": 29149002, "1_0": 29143989, "1_1": 29143910}), |
| 45 | + ("ornithorhynchus_anatinus", "ENSOANT00000062336", -1, |
| 46 | + {"17_2": 29149002, "1_0": 29143989, "1_1": 29143910}), |
| 47 | + # Xenopus (xenopus_tropicalis) transcript |
| 48 | + ("xenopus_tropicalis", "ENSXETT00000106354", -1, |
| 49 | + {"17_2": 3169580, "1_0": 3164877, "1_1": 3164798}), |
| 50 | + ], |
| 51 | +) |
| 52 | +def test_real_negative_strand_order(species, trx, strand, starts): |
| 53 | + # Create three rows with the same SubexonRank to trigger tie-breaking |
| 54 | + # Intentionally insert rows in an order that is NOT genomic order |
| 55 | + # so the test fails if get_tidy_table does not break ties using |
| 56 | + # genomic start positions (via the signed-start sort key). |
| 57 | + insertion_order = ["1_1", "1_0", "17_2"] |
| 58 | + rows = [ |
| 59 | + _mk_row( |
| 60 | + gene="G", |
| 61 | + trx=trx, |
| 62 | + subexon_rank=15, |
| 63 | + strand=strand, |
| 64 | + start=starts[label], |
| 65 | + end=starts[label] - 1, |
| 66 | + label=label, |
| 67 | + ) |
| 68 | + for label in insertion_order |
| 69 | + ] |
| 70 | + pre = pd.DataFrame(rows) |
| 71 | + |
| 72 | + # gene2species mapping is required by get_tidy_table in your ThorAxe branch |
| 73 | + tidy = get_tidy_table(pre, gene2species={"G": species}) |
| 74 | + |
| 75 | + grp = tidy[tidy["TranscriptIDCluster"] == trx].sort_values("S_exon_Rank") |
| 76 | + got = grp["S_exonID"].tolist() |
| 77 | + |
| 78 | + # Expected order: for negative strand, 5′→3′ is decreasing start. |
| 79 | + # The tidy table exposes `SubexonCodingStart` (not `S_exon_CodingStart`). |
| 80 | + signed = np.where( |
| 81 | + grp["Strand"].to_numpy() == 1, |
| 82 | + grp["SubexonCodingStart"].to_numpy(), |
| 83 | + -grp["SubexonCodingStart"].to_numpy(), |
| 84 | + ) |
| 85 | + exp = grp.iloc[np.argsort(signed, kind="stable")]["S_exonID"].tolist() |
| 86 | + |
| 87 | + assert got == exp, f"{species}:{trx} expected {exp} by genomic order, got {got}" |
0 commit comments