Repository: cogent3/cogent3
Pull request #1037 (2022-04-23 12:16:38)

BUG: fixed error in applying off annotations to a sequence


Patch:
diff --git a/src/cogent3/core/location.py b/src/cogent3/core/location.py
index 9f151a24d..a09092574 100644
--- a/src/cogent3/core/location.py
+++ b/src/cogent3/core/location.py
@@ -575,10 +575,6 @@ class Map(object):
                     raise RuntimeError(
                         f"located outside sequence: {str((start, end, parent_length))}"
                     )
-                elif max(start, end) < 0:
-                    diff = min(start, end)
-                    start = [start, 0][start < 0]
-                    end = [end, 0][end < 0]
                 elif min(start, end) > parent_length:
                     diff = max(start, end) - parent_length
                     start = [start, parent_length][start > parent_length]
diff --git a/src/cogent3/core/sequence.py b/src/cogent3/core/sequence.py
index a7be415a2..55bb3d65f 100644
--- a/src/cogent3/core/sequence.py
+++ b/src/cogent3/core/sequence.py
@@ -872,12 +872,24 @@ class Sequence(_Annotatable, SequenceI):
                 annot.copy_annotations_to(new)
         return new
 
+    def _get_feature_start(self, feature):
+        """returns feature offset relative to parent feature(s)"""
+        start = feature.map.start
+        offset = 0
+        while feature.parent is not self:
+            feature = feature.parent
+            if feature.map.start - start:
+                offset += feature.map.start
+
+        return offset + start
+
     def annotate_from_gff(self, f, pre_parsed=False):
         """annotates a Sequence from a gff file where each entry has the same SeqID"""
         # only features with parent features included in the 'features' dict
         gff_contents = f if pre_parsed else gff.gff_parser(f)
         top_level = defaultdict(list)
         grouped = defaultdict(list)
+        num_no_id = 0
         for gff_dict in gff_contents:
             if gff_dict["SeqID"] != self.name:
                 # we can only handle features for this sequence
@@ -885,13 +897,19 @@ class Sequence(_Annotatable, SequenceI):
 
             id_ = gff_dict["Attributes"]["ID"]
             parents = gff_dict["Attributes"].get("Parent", None)
-            if parents is None:
+            if parents is None and id_:
                 assert id_ not in top_level, f"non-unique id {id_}"
                 top_level[id_].append(
                     self.add_feature(
                         gff_dict["Type"], id_, [(gff_dict["Start"], gff_dict["End"])]
                     )
                 )
+            elif parents is None:
+                id_ = f"no-id-{num_no_id}"
+                num_no_id += 1
+                self.add_feature(
+                    gff_dict["Type"], id_, [(gff_dict["Start"], gff_dict["End"])]
+                )
             else:
                 for parent in parents:
                     grouped[parent].append(gff_dict)
@@ -904,33 +922,18 @@ class Sequence(_Annotatable, SequenceI):
                     break
 
             for feature in features:
-                feature_start = feature.map.start
+                feature_start = self._get_feature_start(feature)
                 for gff_dict in child_features:
-                    top_level[gff_dict["Attributes"]["ID"]].append(
-                        feature.add_feature(
-                            gff_dict["Type"],
-                            id_,
-                            [
-                                (
-                                    gff_dict["Start"] - feature_start,
-                                    gff_dict["End"] - feature_start,
-                                )
-                            ],
-                        )
+                    id_ = gff_dict["Attributes"]["ID"]
+                    b = gff_dict["Start"]
+                    e = gff_dict["End"]
+                    type_ = gff_dict["Type"]
+                    sub_feat = feature.add_feature(
+                        type_,
+                        id_,
+                        [(b - feature_start, e - feature_start)],
                     )
-
-    def _sort_parents(self, parents, ordered, key):
-        """returns a list of feature id's with parents before children"""
-        keys = parents.keys()
-        if key in keys:
-            for parent in parents[key]:
-                if parent in keys:
-                    return self._sort_parents(parents, ordered, parent)
-        ordered.append(key)
-        parents.pop(key)
-        if not parents:
-            return ordered
-        return self._sort_parents(parents, ordered, next(iter(keys)))
+                    top_level[gff_dict["Attributes"]["ID"]].append(sub_feat)
 
     def with_masked_annotations(
         self, annot_types, mask_char=None, shadow=False, extend_query=False
@@ -972,8 +975,7 @@ class Sequence(_Annotatable, SequenceI):
         i = 0
         segments = []
         for b, e in region.get_coordinates():
-            segments.append(self._seq[i:b])
-            segments.append(mask_char * (e - b))
+            segments.extend((self._seq[i:b], mask_char * (e - b)))
             i = e
         segments.append(self._seq[i:])
 
@@ -1001,8 +1003,7 @@ class Sequence(_Annotatable, SequenceI):
 
     def gapped_by_map_motif_iter(self, map):
         for segment in self.gapped_by_map_segment_iter(map):
-            for motif in segment:
-                yield motif
+            yield from segment
 
     def gapped_by_map(self, map, recode_gaps=False):
         segments = self.gapped_by_map_segment_iter(map, True, recode_gaps)
@@ -1022,14 +1023,13 @@ class Sequence(_Annotatable, SequenceI):
         myclass = f"{self.__class__.__name__}"
         myclass = myclass.split(".")[-1]
         if len(self) > 10:
-            seq = str(self._seq[:7]) + f"... {len(self)}"
+            seq = f"{str(self._seq[:7])}... {len(self)}"
         else:
             seq = str(self._seq)
         return f"{myclass}({seq})"
 
     def get_name(self):
         """Return the sequence name -- should just use name instead."""
-
         return self.name
 
     def __len__(self):
@@ -1085,17 +1085,17 @@ class Sequence(_Annotatable, SequenceI):
         seq = self._seq
         if motif_length == 1:
             return seq
-        else:
-            length = len(seq)
-            remainder = length % motif_length
-            if remainder and log_warnings:
-                warnings.warn(
-                    f'Dropped remainder "{seq[-remainder:]}" from end of sequence'
-                )
-            return [
-                seq[i : i + motif_length]
-                for i in range(0, length - remainder, motif_length)
-            ]
+
+        length = len(seq)
+        remainder = length % motif_length
+        if remainder and log_warnings:
+            warnings.warn(
+                f'Dropped remainder "{seq[-remainder:]}" from end of sequence'
+            )
+        return [
+            seq[i : i + motif_length]
+            for i in range(0, length - remainder, motif_length)
+        ]
 
     def parse_out_gaps(self):
         gapless = []
diff --git a/src/cogent3/parse/gff.py b/src/cogent3/parse/gff.py
index b8bdedf4c..6af34d767 100644
--- a/src/cogent3/parse/gff.py
+++ b/src/cogent3/parse/gff.py
@@ -1,4 +1,8 @@
 #!/usr/bin/env python
+from pathlib import Path
+
+from cogent3.util.io import open_
+
 
 __author__ = "Peter Maxwell"
 __copyright__ = "Copyright 2007-2022, The Cogent Project"
@@ -14,10 +18,6 @@ __maintainer__ = "Peter Maxwell"
 __email__ = "pm67nz@gmail.com"
 __status__ = "Production"
 
-from pathlib import Path
-
-from cogent3.util.io import open_
-
 
 def gff_parser(f):
     """parses a gff file
@@ -51,6 +51,11 @@ def _gff_parser(f):
         gff3 = gff3_header in f.readline()
         f.seek(0)
 
+    if gff3:
+        attribute_parser = parse_attributes_gff3
+    else:
+        attribute_parser = parse_attributes_gff2
+
     for line in f:
         # comments and blank lines
         if "#" in line:
@@ -66,10 +71,10 @@ def _gff_parser(f):
         if len(cols) == 8:
             cols.append("")
         assert len(cols) == 9, len(line)
-        (seqid, source, type_, start, end, score, strand, phase, attributes) = cols
+        seqid, source, type_, start, end, score, strand, phase, attributes = cols
 
         # adjust for 0-based indexing
-        (start, end) = (int(start) - 1, int(end))
+        start, end = int(start) - 1, int(end)
         # start is always meant to be less than end in GFF
         # features that extend beyond sequence have negative indices
         if start < 0 or end < 0:
@@ -81,10 +86,6 @@ def _gff_parser(f):
             (start, end) = (end, start)
 
         # all attributes have an "ID" but this may not be unique
-        if gff3:
-            attribute_parser = parse_attributes_gff3
-        else:
-            attribute_parser = parse_attributes_gff2
         attributes = attribute_parser(attributes, (start, end))
 
         rtn = {
@@ -114,16 +115,12 @@ def parse_attributes_gff3(attributes, span):
     """Returns a dictionary containing all the attributes"""
     attributes = attributes.strip(";")
     attributes = attributes.split(";")
-    if attributes[0]:
-        attributes = dict(t.split("=") for t in attributes)
-    else:
-        attributes = {}
-    if "Parent" in attributes.keys():
+    attributes = dict(t.split("=") for t in attributes) if attributes[0] else {}
+    if "Parent" in attributes:
         # There may be multiple parents
         if "," in attributes["Parent"]:
             attributes["Parent"] = attributes["Parent"].split(",")
         else:
             attributes["Parent"] = [attributes["Parent"]]
-    if "ID" not in attributes.keys():
-        attributes["ID"] = ""
+    attributes["ID"] = attributes.get("ID", "")
     return attributes


Test patch:
diff --git a/tests/data/ensembl_sample.gff3 b/tests/data/ensembl_sample.gff3
new file mode 100644
index 000000000..55a95aeea
--- /dev/null
+++ b/tests/data/ensembl_sample.gff3
@@ -0,0 +1,17 @@
+##gff-version 3
+##sequence-region   22 1 40
+#!genome-build Genome Reference Consortium GRCh38.p13
+#!genome-version GRCh38
+#!genome-date 2013-12
+#!genome-build-accession GCA_000001405.28
+#!genebuild-last-updated 2021-11
+22	GRCh38	chromosome	1	40	.	.	.	ID=chromosome:22;Alias=CM000684.2,chr22,NC_000022.11
+###
+22	.	biological_region	37	38	2.94e+03	.	.	external_name=oe %3D 0.82;logic_name=cpg
+22	.	biological_region	39	39	0.999	-	.	logic_name=eponine
+###
+22	.	biological_region	2	5	2.94e+03	.	.	external_name=oe %3D 0.82;logic_name=cpg
+22	ensembl_havana	gene	6	36	.	+	.	ID=gene:ENSG00000215568;Name=GAB4;biotype=protein_coding;description=GRB2 associated binding protein family member 4 [Source:HGNC Symbol%3BAcc:HGNC:18325];gene_id=ENSG00000215568;logic_name=ensembl_havana_gene_homo_sapiens;version=9
+22	havana	mRNA	8	30	.	+	.	ID=transcript:ENST00000651146;Parent=gene:ENSG00000215568;Name=GAB4-206;biotype=nonsense_mediated_decay;transcript_id=ENST00000651146;version=1
+22	havana	exon	12	20	.	+	.	Parent=transcript:ENST00000651146;Name=ENSE00003848117;constitutive=0;ensembl_end_phase=-1;ensembl_phase=-1;exon_id=ENSE00003848117;rank=9;version=1
+22	havana	exon	23	27	.	+	.	Parent=transcript:ENST00000651146;Name=ENSE00003513343;constitutive=0;ensembl_end_phase=-1;ensembl_phase=-1;exon_id=ENSE00003513343;rank=8;version=1
diff --git a/tests/test_core/test_sequence.py b/tests/test_core/test_sequence.py
index b9a32e50f..8e19c560f 100644
--- a/tests/test_core/test_sequence.py
+++ b/tests/test_core/test_sequence.py
@@ -248,6 +248,48 @@ class SequenceTests(TestCase):
         # 13 features with one having 2 parents, so 14 instances should be found
         self.assertEqual(len(matches), 14)
 
+    def test_annotate_gff_nested_features(self):
+        """correctly annotate a sequence with nested features"""
+        # the synthetic example
+        #          1111111111222222222333333333334
+        # 1234567890123456789012345678901234567890
+        #  **** biological_region
+        #                                     ** biological_region
+        #                                       * biological_region
+        #      *******************************  gene
+        #         *********************   mRNA
+        #            *********            exon
+        #                       *****     exon
+        # ACCCCGGAAAATTTTTTTTTAAGGGGGAAAAAAAAACCCCCCC...
+        seq = DNA.make_seq("ACCCCGGAAAATTTTTTTTTAAGGGGGAAAAAAAAACCCCCCC", name="22")
+        gff3_path = os.path.join("data/ensembl_sample.gff3")
+        seq.annotate_from_gff(gff3_path)
+        # we have 1 "full chromosome" annotation, 3 generic regions and 1 gene
+        self.assertEqual(len(seq.annotations), 5)
+
+        # get the gene and check it has a single annotation and that
+        # its slice is correct
+        ann = seq.get_annotations_matching("gene")
+        self.assertEqual(len(ann), 1)
+        self.assertEqual(len(ann[0].annotations), 1)
+        seq = ann[0].get_slice()
+        self.assertEqual(str(seq), "GGAAAATTTTTTTTTAAGGGGGAAAAAAAAA")
+
+        # the gene has 1 transcript
+        ann = seq.get_annotations_matching("mRNA", extend_query=True)
+        self.assertEqual(len(ann), 1)
+        self.assertEqual(len(ann[0].annotations), 2)  # 2 exons
+        seq = ann[0].get_slice()
+        self.assertEqual(str(seq), "AAAATTTTTTTTTAAGGGGGAAA")
+
+        # the transcript has 2 exons
+        ann = seq.get_annotations_matching("exon", extend_query=True)
+        self.assertEqual(len(ann), 2)
+        exon_seqs = ("TTTTTTTTT", "GGGGG")
+        for x in ann:
+            self.assertEqual(len(x.annotations), 0)
+            self.assertTrue(str(x.get_slice()) in exon_seqs, msg=x.get_slice())
+
     def test_strip_degenerate(self):
         """Sequence strip_degenerate should remove any degenerate bases"""
         self.assertEqual(self.RNA("UCAG-").strip_degenerate(), "UCAG-")