annotate cpt_blastn_to_gff/blast_to_gff3.py @ 0:54c3aabcb3e7 draft

Uploaded
author cpt
date Fri, 13 May 2022 04:42:45 +0000
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
1 #!/usr/bin/env python
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
2 import argparse
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
3 import copy
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
4 import logging
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
5 import re
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
6 import sys
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
7 from CPT_GFFParser import gffParse, gffWrite, gffSeqFeature
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
8 from Bio.Blast import NCBIXML
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
9 from Bio.Seq import Seq
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
10 from Bio.SeqRecord import SeqRecord
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
11 from Bio.SeqFeature import SeqFeature, FeatureLocation
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
12
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
13 logging.basicConfig(level=logging.INFO)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
14 log = logging.getLogger(name="blast2gff3")
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
15
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
16 __doc__ = """
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
17 Convert BlastXML or Blast 25 Column Table output into GFF3
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
18 """
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
19
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
20 # note for all FeatureLocations, Biopython saves in zero index and Blast provides one indexed locations, thus a Blast Location of (123,500) should be saved as (122, 500)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
21 def blast2gff3(blast, blastxml=False, blasttab=False, include_seq=False):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
22 # Call correct function based on xml or tabular file input, raise error if neither or both are provided
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
23 if blastxml and blasttab:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
24 raise Exception("Cannot provide both blast xml and tabular flag")
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
25
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
26 if blastxml:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
27 return blastxml2gff3(blast, include_seq)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
28 elif blasttab:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
29 return blasttsv2gff3(blast, include_seq)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
30 else:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
31 raise Exception("Must provide either blast xml or tabular flag")
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
32
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
33
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
34 def check_bounds(ps, pe, qs, qe):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
35 # simplify the constant boundary checking used in subfeature generation
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
36 if qs < ps:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
37 ps = qs
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
38 if qe > pe:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
39 pe = qe
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
40 if ps <= 0:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
41 ps = 1
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
42 return (min(ps, pe), max(ps, pe))
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
43
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
44
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
45 def clean_string(s):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
46 clean_str = re.sub("\|", "_", s) # Replace any \ or | with _
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
47 clean_str = re.sub(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
48 "[^A-Za-z0-9_\ .-]", "", clean_str
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
49 ) # Remove any non-alphanumeric or _.- chars
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
50 return clean_str
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
51
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
52
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
53 def clean_slist(l):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
54 cleaned_list = []
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
55 for s in l:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
56 cleaned_list.append(clean_string(s))
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
57 return cleaned_list
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
58
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
59
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
60 def blastxml2gff3(blastxml, include_seq=False):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
61
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
62 blast_records = NCBIXML.parse(blastxml)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
63 for idx_record, record in enumerate(blast_records):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
64 # http://www.sequenceontology.org/browser/release_2.4/term/SO:0000343
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
65 # match_type = { # Currently we can only handle BLASTN, BLASTP
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
66 # "BLASTN": "nucleotide_match",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
67 # "BLASTP": "protein_match",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
68 # }.get(record.application, "match")
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
69 match_type = "match"
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
70 collected_records = []
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
71
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
72 recid = record.query
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
73 if " " in recid:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
74 recid = clean_string(recid[0 : recid.index(" ")])
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
75
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
76 for idx_hit, hit in enumerate(record.alignments):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
77 # gotta check all hsps in a hit to see boundaries
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
78 rec = SeqRecord("", id=recid)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
79 parent_match_start = 0
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
80 parent_match_end = 0
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
81 hit_qualifiers = {
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
82 "ID": "b2g.%s.%s" % (idx_record, idx_hit),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
83 "source": "blast",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
84 "accession": hit.accession,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
85 "hit_id": clean_string(hit.hit_id),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
86 "score": None,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
87 "length": hit.length,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
88 "hit_titles": clean_slist(hit.title.split(" >")),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
89 "hsp_count": len(hit.hsps),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
90 }
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
91 desc = hit.title.split(" >")[0]
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
92 hit_qualifiers["Name"] = desc
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
93 sub_features = []
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
94 for idx_hsp, hsp in enumerate(hit.hsps):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
95 if idx_hsp == 0:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
96 # -2 and +1 for start/end to convert 0 index of python to 1 index of people, -2 on start because feature location saving issue
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
97 parent_match_start = hsp.query_start
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
98 parent_match_end = hsp.query_end
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
99 hit_qualifiers["score"] = hsp.expect
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
100 # generate qualifiers to be added to gff3 feature
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
101 hit_qualifiers["score"] = min(hit_qualifiers["score"], hsp.expect)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
102 hsp_qualifiers = {
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
103 "ID": "b2g.%s.%s.hsp%s" % (idx_record, idx_hit, idx_hsp),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
104 "source": "blast",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
105 "score": hsp.expect,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
106 "accession": hit.accession,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
107 "hit_id": clean_string(hit.hit_id),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
108 "length": hit.length,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
109 "hit_titles": clean_slist(hit.title.split(" >")),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
110 }
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
111 if include_seq:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
112 if (
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
113 "blast_qseq",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
114 "blast_sseq",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
115 "blast_mseq",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
116 ) in hit_qualifiers.keys():
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
117 hit_qualifiers.update(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
118 {
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
119 "blast_qseq": hit_qualifiers["blast_qseq"] + hsp.query,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
120 "blast_sseq": hit_qualifiers["blast_sseq"] + hsp.sbjct,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
121 "blast_mseq": hit_qualifiers["blast_mseq"] + hsp.match,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
122 }
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
123 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
124 else:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
125 hit_qualifiers.update(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
126 {
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
127 "blast_qseq": hsp.query,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
128 "blast_sseq": hsp.sbjct,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
129 "blast_mseq": hsp.match,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
130 }
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
131 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
132 for prop in (
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
133 "score",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
134 "bits",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
135 "identities",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
136 "positives",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
137 "gaps",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
138 "align_length",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
139 "strand",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
140 "frame",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
141 "query_start",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
142 "query_end",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
143 "sbjct_start",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
144 "sbjct_end",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
145 ):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
146 hsp_qualifiers["blast_" + prop] = getattr(hsp, prop, None)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
147
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
148 # check if parent boundary needs to increase to envelope hsp
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
149 # if hsp.query_start < parent_match_start:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
150 # parent_match_start = hsp.query_start - 1
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
151 # if hsp.query_end > parent_match_end:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
152 # parent_match_end = hsp.query_end + 1
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
153
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
154 parent_match_start, parent_match_end = check_bounds(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
155 parent_match_start, parent_match_end, hsp.query_start, hsp.query_end
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
156 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
157
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
158 # add hsp to the gff3 feature as a "match_part"
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
159 sub_features.append(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
160 gffSeqFeature(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
161 FeatureLocation(hsp.query_start - 1, hsp.query_end),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
162 type="match_part",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
163 strand=0,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
164 qualifiers=copy.deepcopy(hsp_qualifiers),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
165 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
166 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
167
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
168 # Build the top level seq feature for the hit
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
169 hit_qualifiers["description"] = "Residue %s..%s hit to %s" % (parent_match_start, parent_match_end, desc,)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
170 top_feature = gffSeqFeature(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
171 FeatureLocation(parent_match_start - 1, parent_match_end),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
172 type=match_type,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
173 strand=0,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
174 qualifiers=hit_qualifiers,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
175 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
176 # add the generated subfeature hsp match_parts to the hit feature
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
177 top_feature.sub_features = copy.deepcopy(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
178 sorted(sub_features, key=lambda x: int(x.location.start))
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
179 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
180 # Add the hit feature to the record
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
181 rec.features.append(top_feature)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
182 rec.annotations = {}
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
183 collected_records.append(rec)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
184
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
185 if not len(collected_records):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
186 print("##gff-version 3\n##sequence-region null 1 4")
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
187
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
188 for rec in collected_records:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
189 yield rec
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
190
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
191
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
192 def combine_records(records):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
193 # Go through each record and identify those records with
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
194 cleaned_records = {}
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
195 for rec in records:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
196 combo_id = (
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
197 rec.features[0].qualifiers["target"]
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
198 + rec.features[0].qualifiers["accession"]
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
199 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
200 if combo_id not in cleaned_records.keys():
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
201 # First instance of a query ID + subject ID combination
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
202 # Save this record as it's only item
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
203 newid = rec.features[0].qualifiers["ID"] + ".0"
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
204 rec.features[0].qualifiers["ID"] = newid
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
205 rec.features[0].sub_features[0].qualifiers["ID"] = newid + ".hsp0"
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
206 cleaned_records[combo_id] = rec
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
207 else:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
208 # Query ID + Subject ID has appeared before
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
209 # Combine the Match Parts as subfeatures
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
210 sub_features = copy.deepcopy(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
211 cleaned_records[combo_id].features[0].sub_features
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
212 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
213 addtnl_features = rec.features[0].sub_features
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
214 # add the current records sub features to the ones previous
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
215 for feat in addtnl_features:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
216 sub_features.append(feat)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
217 cleaned_records[combo_id].features[0].subfeatures = copy.deepcopy(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
218 sub_features
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
219 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
220 cleaned_records[combo_id].features[0].qualifiers["score"] = min(cleaned_records[combo_id].features[0].qualifiers["score"], rec.features[0].qualifiers["score"])
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
221 # now we need to update the IDs for the features when combined
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
222 # sort them into the proper order, then apply new ids
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
223 # and also ensure the parent record boundaries fit the whole span of subfeatures
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
224 sub_features = sorted(sub_features, key=lambda x: int(x.location.start))
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
225 new_parent_start = cleaned_records[combo_id].features[0].location.start + 1
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
226 new_parent_end = cleaned_records[combo_id].features[0].location.end
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
227 for idx, feat in enumerate(sub_features):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
228 feat.qualifiers["ID"] = "%s.hsp%s" % (
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
229 cleaned_records[combo_id].features[0].qualifiers["ID"],
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
230 idx,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
231 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
232 new_parent_start, new_parent_end = check_bounds(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
233 new_parent_start,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
234 new_parent_end,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
235 feat.location.start + 1,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
236 feat.location.end,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
237 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
238 cleaned_records[combo_id].features[0].qualifiers["score"] = min(cleaned_records[combo_id].features[0].qualifiers["score"], feat.qualifiers["blast_score"])
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
239 # if feat.location.start < new_parent_start:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
240 # new_parent_start = feat.location.start - 1
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
241 # if feat.location.end > new_parent_end:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
242 # new_parent_end = feat.location.end + 1
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
243 cleaned_records[combo_id].features[0].location = FeatureLocation(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
244 new_parent_start - 1, new_parent_end
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
245 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
246 cleaned_records[combo_id].features[0].qualifiers[
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
247 "description"
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
248 ] = "Residue %s..%s hit to %s" % (
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
249 new_parent_start,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
250 new_parent_end,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
251 cleaned_records[combo_id].features[0].qualifiers["Name"],
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
252 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
253 # save the renamed and ordered feature list to record
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
254 cleaned_records[combo_id].features[0].sub_features = copy.deepcopy(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
255 sub_features
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
256 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
257 return sorted(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
258 cleaned_records.values(), key=lambda x: int(x.features[0].location.start)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
259 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
260
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
261
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
262 def blasttsv2gff3(blasttsv, include_seq=False):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
263
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
264 # http://www.sequenceontology.org/browser/release_2.4/term/SO:0000343
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
265 # match_type = { # Currently we can only handle BLASTN, BLASTP
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
266 # "BLASTN": "nucleotide_match",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
267 # "BLASTP": "protein_match",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
268 # }.get(type, "match")
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
269 match_type = "match"
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
270
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
271 columns = [
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
272 "qseqid", # 01 Query Seq-id (ID of your sequence)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
273 "sseqid", # 02 Subject Seq-id (ID of the database hit)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
274 "pident", # 03 Percentage of identical matches
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
275 "length", # 04 Alignment length
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
276 "mismatch", # 05 Number of mismatches
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
277 "gapopen", # 06 Number of gap openings
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
278 "qstart", # 07 Start of alignment in query
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
279 "qend", # 08 End of alignment in query
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
280 "sstart", # 09 Start of alignment in subject (database hit)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
281 "send", # 10 End of alignment in subject (database hit)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
282 "evalue", # 11 Expectation value (E-value)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
283 "bitscore", # 12 Bit score
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
284 "sallseqid", # 13 All subject Seq-id(s), separated by a ';'
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
285 "score", # 14 Raw score
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
286 "nident", # 15 Number of identical matches
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
287 "positive", # 16 Number of positive-scoring matches
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
288 "gaps", # 17 Total number of gaps
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
289 "ppos", # 18 Percentage of positive-scoring matches
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
290 "qframe", # 19 Query frame
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
291 "sframe", # 20 Subject frame
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
292 "qseq", # 21 Aligned part of query sequence
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
293 "sseq", # 22 Aligned part of subject sequence
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
294 "qlen", # 23 Query sequence length
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
295 "slen", # 24 Subject sequence length
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
296 "salltitles", # 25 All subject title(s), separated by a '<>'
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
297 ]
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
298 collected_records = []
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
299 for record_idx, record in enumerate(blasttsv):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
300 if record.startswith("#"):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
301 continue
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
302
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
303 dc = {k: v for (k, v) in zip(columns, (x.strip() for x in record.split("\t")))}
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
304
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
305 rec = SeqRecord("", id=dc["qseqid"])
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
306
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
307 feature_id = "b2g.%s" % (record_idx)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
308 hit_qualifiers = {
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
309 "ID": feature_id,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
310 "Name": (dc["salltitles"].split("<>")[0]),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
311 "description": "Residue {sstart}..{send} hit to {x}".format(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
312 x=dc["salltitles"].split("<>")[0], **dc
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
313 ),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
314 "source": "blast",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
315 "score": dc["evalue"],
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
316 "accession": clean_string(dc["sseqid"]),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
317 "length": dc["qlen"],
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
318 "hit_titles": clean_slist(dc["salltitles"].split("<>")),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
319 "target": clean_string(dc["qseqid"]),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
320 }
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
321 hsp_qualifiers = {"source": "blast"}
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
322 for key in dc.keys():
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
323 # Add the remaining BLAST info to the GFF qualifiers
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
324 if key in ("salltitles", "sallseqid", "sseqid", "qseqid", "qseq", "sseq",):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
325 continue
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
326 hsp_qualifiers["blast_%s" % key] = clean_string(dc[key])
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
327
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
328 # Below numbers stored as strings, convert to proper form
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
329 for (
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
330 integer_numerical_key
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
331 ) in "gapopen gaps length mismatch nident positive qend qframe qlen qstart score send sframe slen sstart".split(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
332 " "
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
333 ):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
334 dc[integer_numerical_key] = int(dc[integer_numerical_key])
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
335
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
336 for float_numerical_key in "bitscore evalue pident ppos".split(" "):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
337 dc[float_numerical_key] = float(dc[float_numerical_key])
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
338
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
339 parent_match_start = dc["qstart"]
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
340 parent_match_end = dc["qend"]
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
341
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
342 parent_match_start, parent_match_end = check_bounds(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
343 parent_match_start, parent_match_end, dc["qstart"], dc["qend"]
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
344 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
345
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
346 # The ``match`` feature will hold one or more ``match_part``s
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
347 top_feature = gffSeqFeature(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
348 FeatureLocation(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
349 min(parent_match_start, parent_match_end) - 1,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
350 max(parent_match_start, parent_match_end),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
351 ),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
352 type=match_type,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
353 strand=0,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
354 qualifiers=hit_qualifiers,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
355 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
356 top_feature.sub_features = []
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
357 # There is a possibility of multiple lines containing the HSPS
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
358 # for the same hit.
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
359 # Unlike the parent feature, ``match_part``s have sources.
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
360 hsp_qualifiers["ID"] = clean_string(dc["sseqid"])
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
361 match_part_start = dc["qstart"]
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
362 match_part_end = dc["qend"]
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
363
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
364 top_feature.sub_features.append(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
365 gffSeqFeature(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
366 FeatureLocation(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
367 min(match_part_start, match_part_end) - 1,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
368 max(match_part_start, match_part_end),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
369 ),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
370 type="match_part",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
371 strand=0,
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
372 qualifiers=copy.deepcopy(hsp_qualifiers),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
373 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
374 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
375 top_feature.sub_features = sorted(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
376 top_feature.sub_features, key=lambda x: int(x.location.start)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
377 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
378 rec.features = [top_feature]
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
379 rec.annotations = {}
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
380 collected_records.append(rec)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
381
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
382 collected_records = combine_records(collected_records)
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
383 if not len(collected_records):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
384 print("##gff-version 3\n##sequence-region null 1 4")
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
385 for rec in collected_records:
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
386 yield rec
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
387
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
388
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
389 if __name__ == "__main__":
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
390 parser = argparse.ArgumentParser(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
391 description="Convert BlastP or BlastN output to GFF3, must provide XML or Tabular output",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
392 epilog="",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
393 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
394 parser.add_argument(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
395 "blast",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
396 type=argparse.FileType("r"),
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
397 help="Blast XML or 25 Column Tabular Output file",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
398 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
399 parser.add_argument(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
400 "--blastxml", action="store_true", help="Process file as Blast XML Output"
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
401 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
402 parser.add_argument(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
403 "--blasttab",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
404 action="store_true",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
405 help="Process file as Blast 25 Column Tabular Output",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
406 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
407 parser.add_argument(
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
408 "--include_seq",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
409 action="store_true",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
410 help="Include sequence, only used for Blast XML",
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
411 )
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
412 args = parser.parse_args()
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
413
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
414 for rec in blast2gff3(**vars(args)):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
415 if len(rec.features):
54c3aabcb3e7 Uploaded
cpt
parents:
diff changeset
416 gffWrite([rec], sys.stdout)