diff --git a/gfftk/gff.py b/gfftk/gff.py index 5b47936..301d2a3 100755 --- a/gfftk/gff.py +++ b/gfftk/gff.py @@ -2137,7 +2137,7 @@ def _longest_orf(annot, fadict, minlen=50, table=1): CDS.append((start_pos, f)) cov += f - start_pos + 1 elif len(CDS) > 0: - if cov <= lenOrf: + if cov < lenOrf: remainder = lenOrf - cov if (f - s) < remainder: CDS.append((s, f)) @@ -2155,13 +2155,21 @@ def _longest_orf(annot, fadict, minlen=50, table=1): v["CDS"] = [sorted(CDS, key=lambda tup: tup[0])] else: v["CDS"] = [sorted(CDS, key=lambda tup: tup[0], reverse=True)] - v["phase"] = ["?"] + v["phase"] = [["?"] * len(CDS)] + v["partialStart"] = [False] + v["partialStop"] = [False] Clean[k] = v except AssertionError: + v["partialStart"] = [False] + v["partialStop"] = [False] Clean[k] = v else: + v["partialStart"] = [False] + v["partialStop"] = [False] Clean[k] = v else: # did not find orfs so add as is + v["partialStart"] = [False] + v["partialStop"] = [False] Clean[k] = v return Clean diff --git a/tests/unit/test_gff.py b/tests/unit/test_gff.py index 444bf89..ef5a0b9 100644 --- a/tests/unit/test_gff.py +++ b/tests/unit/test_gff.py @@ -72,6 +72,8 @@ def test_simplify_go(self): assert result == [] + + class TestCombinedGFFFormat: """Tests for combined GFF3+FASTA format functionality.""" @@ -349,3 +351,106 @@ def test_sgd_gene_structure(self): finally: os.unlink(gff_name) os.unlink(fasta_name) + + def test_longest_orf_gapmm2_alignment(self): + """Test _longest_orf properly handles partialStart/Stop and exact boundary conditions.""" + # Create a fasta file + fd, fasta_file = tempfile.mkstemp(suffix=".fasta") + os.close(fd) + + # Create a mock gapmm2 alignment gff3 file + fd2, gff_file = tempfile.mkstemp(suffix=".gff3") + os.close(fd2) + + try: + # We want an exact matching coverage length + # The bug was cov <= lenOrf triggering a zero/negative length CDS + # mrna_seq = 102 bases, which covers exactly exon 1 (51) and exon 2 (51) + # The third exon should not be included in the CDS. + mrna_seq = "ATG" + "A" * 96 + "TAG" + seq = ["N"] * 200 + seq[9:60] = list(mrna_seq[0:51]) + seq[69:120] = list(mrna_seq[51:102]) + seq[129:186] = ["A"] * 57 + fasta_str = "".join(seq) + + with open(fasta_file, "w") as f: + f.write(">chr1\n") + f.write(fasta_str + "\n") + + with open(gff_file, "w") as f: + # 3 exons: (10, 60), (70, 120), (130, 186) + f.write("chr1\tgapmm2\tcDNA_match\t10\t60\t100\t+\t.\tID=gapmm2_4;Target=transcript1 1 51 +\n") + f.write("chr1\tgapmm2\tcDNA_match\t70\t120\t100\t+\t.\tID=gapmm2_4;Target=transcript1 52 102 +\n") + f.write("chr1\tgapmm2\tcDNA_match\t130\t186\t100\t+\t.\tID=gapmm2_4;Target=transcript1 103 159 +\n") + + # Use gff_format="alignment" which triggers the _longest_orf function + parsed_dict = gff2dict(gff_file, fasta_file, gff_format="alignment", debug=False) + + assert "gapmm2_4" in parsed_dict + gene_data = parsed_dict["gapmm2_4"] + + # Since longestORF length is 102, which is NOT > minlen*3 (150) + # It will fall to the else block. Let's make sure it handles it without crashing. + # partialStart and partialStop should be properly initialized + assert gene_data["partialStart"] in ([False], [None]) + assert gene_data["partialStop"] in ([False], [None]) + + finally: + if os.path.exists(fasta_file): + os.remove(fasta_file) + if os.path.exists(gff_file): + os.remove(gff_file) + + def test_longest_orf_exact_boundary(self): + """Test _longest_orf correctly processes an ORF that exactly ends on an exon boundary.""" + # Create a fasta file + fd, fasta_file = tempfile.mkstemp(suffix=".fasta") + os.close(fd) + + # Create a mock gapmm2 alignment gff3 file + fd2, gff_file = tempfile.mkstemp(suffix=".gff3") + os.close(fd2) + + try: + # mrna_seq = 156 bases (52 codons) -> > 150 (minlen 50) + # 156 bases exactly covers exon 1 (51), exon 2 (51), exon 3 (54) + # wait, let's make it cover exon 1 (51) and exon 2 (105) = 156. + # Then an extra exon 3. + mrna_seq = "ATG" + "A" * 150 + "TAG" + seq = ["N"] * 400 + seq[9:60] = list(mrna_seq[0:51]) + seq[69:174] = list(mrna_seq[51:156]) + seq[199:250] = ["A"] * 51 + fasta_str = "".join(seq) + + with open(fasta_file, "w") as f: + f.write(">chr1\n") + f.write(fasta_str + "\n") + + with open(gff_file, "w") as f: + f.write("chr1\tgapmm2\tcDNA_match\t10\t60\t100\t+\t.\tID=gapmm2_long;Target=transcript1 1 51 +\n") + f.write("chr1\tgapmm2\tcDNA_match\t70\t174\t100\t+\t.\tID=gapmm2_long;Target=transcript1 52 156 +\n") + f.write("chr1\tgapmm2\tcDNA_match\t200\t250\t100\t+\t.\tID=gapmm2_long;Target=transcript1 157 207 +\n") + + parsed_dict = gff2dict(gff_file, fasta_file, gff_format="alignment", debug=False) + + assert "gapmm2_long" in parsed_dict + gene_data = parsed_dict["gapmm2_long"] + + # Since longestORF length is 156, it will successfully find the ORF. + # partialStart and partialStop should be properly initialized. + assert len(gene_data["partialStart"]) == 1 + assert len(gene_data["partialStop"]) == 1 + + # Ensure the CDS has only 2 elements, not 3 (where the third is reversed) + assert len(gene_data["CDS"][0]) == 2 + assert gene_data["CDS"][0] == [(10, 60), (70, 174)] + + finally: + if os.path.exists(fasta_file): + os.remove(fasta_file) + if os.path.exists(gff_file): + os.remove(gff_file) + +