Mercurial > repos > cpt > cpt_gff_rebase

comparison gff3_rebase.py @ 1:4f4b413056f6 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 94b0cd1fff0826c6db3e7dc0c91c0c5a8be8bb0c
author cpt
date Mon, 05 Jun 2023 02:44:12 +0000
parents
children
comparison
equal deleted inserted replaced
0:6e7e20cb1fc7 1:4f4b413056f6
1 #!/usr/bin/env python
2 import sys
3 import logging
4 import argparse
5 from gff3 import feature_lambda, feature_test_qual_value
6 from CPT_GFFParser import gffParse, gffWrite
7 from Bio.SeqFeature import FeatureLocation
8
9 log = logging.getLogger(__name__)
10 logging.basicConfig(level=logging.INFO)
11
12
13 def __get_features(child, interpro=False):
14 child_features = {}
15 for rec in gffParse(child):
16 log.info("Parsing %s", rec.id)
17 # Only top level
18 for feature in rec.features:
19 # Get the record id as parent_feature_id (since this is how it will be during remapping)
20 parent_feature_id = rec.id
21 # If it's an interpro specific gff3 file
22 if interpro:
23 # Then we ignore polypeptide features as they're useless
24 if feature.type == "polypeptide":
25 continue
26
27 try:
28 child_features[parent_feature_id].append(feature)
29 except KeyError:
30 child_features[parent_feature_id] = [feature]
31 # Keep a list of feature objects keyed by parent record id
32 return child_features
33
34
35 def __update_feature_location(feature, parent, protein2dna):
36 start = feature.location.start
37 end = feature.location.end
38 if protein2dna:
39 start *= 3
40 end *= 3
41
42 if parent.location.strand >= 0:
43 ns = parent.location.start + start
44 ne = parent.location.start + end
45 st = +1
46 else:
47 ns = parent.location.end - end
48 ne = parent.location.end - start
49 st = -1
50
51 # Don't let start/stops be less than zero.
52 #
53 # Instead, we'll replace with %3 to try and keep it in the same reading
54 # frame that it should be in.
55
56 if ns < 0:
57 ns %= 3
58 if ne < 0:
59 ne %= 3
60
61 feature.location = FeatureLocation(ns, ne, strand=st)
62
63 if hasattr(feature, "sub_features"):
64 for subfeature in feature.sub_features:
65 __update_feature_location(subfeature, parent, protein2dna)
66
67
68 def rebase(parent, child, interpro=False, protein2dna=False, map_by="ID"):
69 # get all of the features we will be re-mapping in a dictionary, keyed by parent feature ID
70 child_features = __get_features(child, interpro=interpro)
71
72 for rec in gffParse(parent):
73 replacement_features = []
74 # Horrifically slow I believe
75 for feature in feature_lambda(
76 rec.features,
77 # Filter features in the parent genome by those that are
78 # "interesting", i.e. have results in child_features array.
79 # Probably an unnecessary optimisation.
80 feature_test_qual_value,
81 {"qualifier": map_by, "attribute_list": child_features.keys()},
82 subfeatures=False,
83 ):
84
85 # Features which will be re-mapped
86 to_remap = child_features[feature.id]
87
88 fixed_features = []
89 for x in to_remap:
90 # Then update the location of the actual feature
91 __update_feature_location(x, feature, protein2dna)
92
93 if interpro:
94 for y in ("status", "Target"):
95 try:
96 del x.qualifiers[y]
97 except:
98 pass
99
100 fixed_features.append(x)
101 replacement_features.extend(fixed_features)
102 # We do this so we don't include the original set of features that we
103 # were rebasing against in our result.
104 rec.features = replacement_features
105 rec.annotations = {}
106 gffWrite([rec], sys.stdout)
107
108
109 if __name__ == "__main__":
110 parser = argparse.ArgumentParser(
111 description="rebase gff3 features against parent locations", epilog=""
112 )
113 parser.add_argument(
114 "parent", type=argparse.FileType("r"), help="Parent GFF3 annotations"
115 )
116 parser.add_argument(
117 "child",
118 type=argparse.FileType("r"),
119 help="Child GFF3 annotations to rebase against parent",
120 )
121 parser.add_argument(
122 "--interpro", action="store_true", help="Interpro specific modifications"
123 )
124 parser.add_argument(
125 "--protein2dna",
126 action="store_true",
127 help="Map protein translated results to original DNA data",
128 )
129 parser.add_argument("--map_by", help="Map by key", default="ID")
130 args = parser.parse_args()
131 rebase(**vars(args))