annotate blastxml_to_tabular.py @ 2:ab1a8640f817 draft

Uploaded v0.0.12 again, without extra path
author peterjc
date Thu, 23 Aug 2012 07:32:06 -0400
parents d375502056f1
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
d375502056f1 Uploaded
devteam
parents:
diff changeset
1 #!/usr/bin/env python
d375502056f1 Uploaded
devteam
parents:
diff changeset
2 """Convert a BLAST XML file to 12 column tabular output
d375502056f1 Uploaded
devteam
parents:
diff changeset
3
d375502056f1 Uploaded
devteam
parents:
diff changeset
4 Takes three command line options, input BLAST XML filename, output tabular
d375502056f1 Uploaded
devteam
parents:
diff changeset
5 BLAST filename, output format (std for standard 12 columns, or ext for the
d375502056f1 Uploaded
devteam
parents:
diff changeset
6 extended 24 columns offered in the BLAST+ wrappers).
d375502056f1 Uploaded
devteam
parents:
diff changeset
7
d375502056f1 Uploaded
devteam
parents:
diff changeset
8 The 12 columns output are 'qseqid sseqid pident length mismatch gapopen qstart
d375502056f1 Uploaded
devteam
parents:
diff changeset
9 qend sstart send evalue bitscore' or 'std' at the BLAST+ command line, which
d375502056f1 Uploaded
devteam
parents:
diff changeset
10 mean:
d375502056f1 Uploaded
devteam
parents:
diff changeset
11
d375502056f1 Uploaded
devteam
parents:
diff changeset
12 ====== ========= ============================================
d375502056f1 Uploaded
devteam
parents:
diff changeset
13 Column NCBI name Description
d375502056f1 Uploaded
devteam
parents:
diff changeset
14 ------ --------- --------------------------------------------
d375502056f1 Uploaded
devteam
parents:
diff changeset
15 1 qseqid Query Seq-id (ID of your sequence)
d375502056f1 Uploaded
devteam
parents:
diff changeset
16 2 sseqid Subject Seq-id (ID of the database hit)
d375502056f1 Uploaded
devteam
parents:
diff changeset
17 3 pident Percentage of identical matches
d375502056f1 Uploaded
devteam
parents:
diff changeset
18 4 length Alignment length
d375502056f1 Uploaded
devteam
parents:
diff changeset
19 5 mismatch Number of mismatches
d375502056f1 Uploaded
devteam
parents:
diff changeset
20 6 gapopen Number of gap openings
d375502056f1 Uploaded
devteam
parents:
diff changeset
21 7 qstart Start of alignment in query
d375502056f1 Uploaded
devteam
parents:
diff changeset
22 8 qend End of alignment in query
d375502056f1 Uploaded
devteam
parents:
diff changeset
23 9 sstart Start of alignment in subject (database hit)
d375502056f1 Uploaded
devteam
parents:
diff changeset
24 10 send End of alignment in subject (database hit)
d375502056f1 Uploaded
devteam
parents:
diff changeset
25 11 evalue Expectation value (E-value)
d375502056f1 Uploaded
devteam
parents:
diff changeset
26 12 bitscore Bit score
d375502056f1 Uploaded
devteam
parents:
diff changeset
27 ====== ========= ============================================
d375502056f1 Uploaded
devteam
parents:
diff changeset
28
d375502056f1 Uploaded
devteam
parents:
diff changeset
29 The additional columns offered in the Galaxy BLAST+ wrappers are:
d375502056f1 Uploaded
devteam
parents:
diff changeset
30
d375502056f1 Uploaded
devteam
parents:
diff changeset
31 ====== ============= ===========================================
d375502056f1 Uploaded
devteam
parents:
diff changeset
32 Column NCBI name Description
d375502056f1 Uploaded
devteam
parents:
diff changeset
33 ------ ------------- -------------------------------------------
d375502056f1 Uploaded
devteam
parents:
diff changeset
34 13 sallseqid All subject Seq-id(s), separated by a ';'
d375502056f1 Uploaded
devteam
parents:
diff changeset
35 14 score Raw score
d375502056f1 Uploaded
devteam
parents:
diff changeset
36 15 nident Number of identical matches
d375502056f1 Uploaded
devteam
parents:
diff changeset
37 16 positive Number of positive-scoring matches
d375502056f1 Uploaded
devteam
parents:
diff changeset
38 17 gaps Total number of gaps
d375502056f1 Uploaded
devteam
parents:
diff changeset
39 18 ppos Percentage of positive-scoring matches
d375502056f1 Uploaded
devteam
parents:
diff changeset
40 19 qframe Query frame
d375502056f1 Uploaded
devteam
parents:
diff changeset
41 20 sframe Subject frame
d375502056f1 Uploaded
devteam
parents:
diff changeset
42 21 qseq Aligned part of query sequence
d375502056f1 Uploaded
devteam
parents:
diff changeset
43 22 sseq Aligned part of subject sequence
d375502056f1 Uploaded
devteam
parents:
diff changeset
44 23 qlen Query sequence length
d375502056f1 Uploaded
devteam
parents:
diff changeset
45 24 slen Subject sequence length
d375502056f1 Uploaded
devteam
parents:
diff changeset
46 ====== ============= ===========================================
d375502056f1 Uploaded
devteam
parents:
diff changeset
47
d375502056f1 Uploaded
devteam
parents:
diff changeset
48 Most of these fields are given explicitly in the XML file, others some like
d375502056f1 Uploaded
devteam
parents:
diff changeset
49 the percentage identity and the number of gap openings must be calculated.
d375502056f1 Uploaded
devteam
parents:
diff changeset
50
d375502056f1 Uploaded
devteam
parents:
diff changeset
51 Be aware that the sequence in the extended tabular output or XML direct from
d375502056f1 Uploaded
devteam
parents:
diff changeset
52 BLAST+ may or may not use XXXX masking on regions of low complexity. This
d375502056f1 Uploaded
devteam
parents:
diff changeset
53 can throw the off the calculation of percentage identity and gap openings.
d375502056f1 Uploaded
devteam
parents:
diff changeset
54 [In fact, both BLAST 2.2.24+ and 2.2.25+ have a subtle bug in this regard,
d375502056f1 Uploaded
devteam
parents:
diff changeset
55 with these numbers changing depending on whether or not the low complexity
d375502056f1 Uploaded
devteam
parents:
diff changeset
56 filter is used.]
d375502056f1 Uploaded
devteam
parents:
diff changeset
57
d375502056f1 Uploaded
devteam
parents:
diff changeset
58 This script attempts to produce identical output to what BLAST+ would have done.
d375502056f1 Uploaded
devteam
parents:
diff changeset
59 However, check this with "diff -b ..." since BLAST+ sometimes includes an extra
d375502056f1 Uploaded
devteam
parents:
diff changeset
60 space character (probably a bug).
d375502056f1 Uploaded
devteam
parents:
diff changeset
61 """
d375502056f1 Uploaded
devteam
parents:
diff changeset
62 import sys
d375502056f1 Uploaded
devteam
parents:
diff changeset
63 import re
d375502056f1 Uploaded
devteam
parents:
diff changeset
64
d375502056f1 Uploaded
devteam
parents:
diff changeset
65 if sys.version_info[:2] >= ( 2, 5 ):
d375502056f1 Uploaded
devteam
parents:
diff changeset
66 import xml.etree.cElementTree as ElementTree
d375502056f1 Uploaded
devteam
parents:
diff changeset
67 else:
d375502056f1 Uploaded
devteam
parents:
diff changeset
68 from galaxy import eggs
d375502056f1 Uploaded
devteam
parents:
diff changeset
69 import pkg_resources; pkg_resources.require( "elementtree" )
d375502056f1 Uploaded
devteam
parents:
diff changeset
70 from elementtree import ElementTree
d375502056f1 Uploaded
devteam
parents:
diff changeset
71
d375502056f1 Uploaded
devteam
parents:
diff changeset
72 def stop_err( msg ):
d375502056f1 Uploaded
devteam
parents:
diff changeset
73 sys.stderr.write("%s\n" % msg)
d375502056f1 Uploaded
devteam
parents:
diff changeset
74 sys.exit(1)
d375502056f1 Uploaded
devteam
parents:
diff changeset
75
d375502056f1 Uploaded
devteam
parents:
diff changeset
76 #Parse Command Line
d375502056f1 Uploaded
devteam
parents:
diff changeset
77 try:
d375502056f1 Uploaded
devteam
parents:
diff changeset
78 in_file, out_file, out_fmt = sys.argv[1:]
d375502056f1 Uploaded
devteam
parents:
diff changeset
79 except:
d375502056f1 Uploaded
devteam
parents:
diff changeset
80 stop_err("Expect 3 arguments: input BLAST XML file, output tabular file, out format (std or ext)")
d375502056f1 Uploaded
devteam
parents:
diff changeset
81
d375502056f1 Uploaded
devteam
parents:
diff changeset
82 if out_fmt == "std":
d375502056f1 Uploaded
devteam
parents:
diff changeset
83 extended = False
d375502056f1 Uploaded
devteam
parents:
diff changeset
84 elif out_fmt == "x22":
d375502056f1 Uploaded
devteam
parents:
diff changeset
85 stop_err("Format argument x22 has been replaced with ext (extended 24 columns)")
d375502056f1 Uploaded
devteam
parents:
diff changeset
86 elif out_fmt == "ext":
d375502056f1 Uploaded
devteam
parents:
diff changeset
87 extended = True
d375502056f1 Uploaded
devteam
parents:
diff changeset
88 else:
d375502056f1 Uploaded
devteam
parents:
diff changeset
89 stop_err("Format argument should be std (12 column) or ext (extended 24 columns)")
d375502056f1 Uploaded
devteam
parents:
diff changeset
90
d375502056f1 Uploaded
devteam
parents:
diff changeset
91
d375502056f1 Uploaded
devteam
parents:
diff changeset
92 # get an iterable
d375502056f1 Uploaded
devteam
parents:
diff changeset
93 try:
d375502056f1 Uploaded
devteam
parents:
diff changeset
94 context = ElementTree.iterparse(in_file, events=("start", "end"))
d375502056f1 Uploaded
devteam
parents:
diff changeset
95 except:
d375502056f1 Uploaded
devteam
parents:
diff changeset
96 stop_err("Invalid data format.")
d375502056f1 Uploaded
devteam
parents:
diff changeset
97 # turn it into an iterator
d375502056f1 Uploaded
devteam
parents:
diff changeset
98 context = iter(context)
d375502056f1 Uploaded
devteam
parents:
diff changeset
99 # get the root element
d375502056f1 Uploaded
devteam
parents:
diff changeset
100 try:
d375502056f1 Uploaded
devteam
parents:
diff changeset
101 event, root = context.next()
d375502056f1 Uploaded
devteam
parents:
diff changeset
102 except:
d375502056f1 Uploaded
devteam
parents:
diff changeset
103 stop_err( "Invalid data format." )
d375502056f1 Uploaded
devteam
parents:
diff changeset
104
d375502056f1 Uploaded
devteam
parents:
diff changeset
105
d375502056f1 Uploaded
devteam
parents:
diff changeset
106 re_default_query_id = re.compile("^Query_\d+$")
d375502056f1 Uploaded
devteam
parents:
diff changeset
107 assert re_default_query_id.match("Query_101")
d375502056f1 Uploaded
devteam
parents:
diff changeset
108 assert not re_default_query_id.match("Query_101a")
d375502056f1 Uploaded
devteam
parents:
diff changeset
109 assert not re_default_query_id.match("MyQuery_101")
d375502056f1 Uploaded
devteam
parents:
diff changeset
110 re_default_subject_id = re.compile("^Subject_\d+$")
d375502056f1 Uploaded
devteam
parents:
diff changeset
111 assert re_default_subject_id.match("Subject_1")
d375502056f1 Uploaded
devteam
parents:
diff changeset
112 assert not re_default_subject_id.match("Subject_")
d375502056f1 Uploaded
devteam
parents:
diff changeset
113 assert not re_default_subject_id.match("Subject_12a")
d375502056f1 Uploaded
devteam
parents:
diff changeset
114 assert not re_default_subject_id.match("TheSubject_1")
d375502056f1 Uploaded
devteam
parents:
diff changeset
115
d375502056f1 Uploaded
devteam
parents:
diff changeset
116
d375502056f1 Uploaded
devteam
parents:
diff changeset
117 outfile = open(out_file, 'w')
d375502056f1 Uploaded
devteam
parents:
diff changeset
118 blast_program = None
d375502056f1 Uploaded
devteam
parents:
diff changeset
119 for event, elem in context:
d375502056f1 Uploaded
devteam
parents:
diff changeset
120 if event == "end" and elem.tag == "BlastOutput_program":
d375502056f1 Uploaded
devteam
parents:
diff changeset
121 blast_program = elem.text
d375502056f1 Uploaded
devteam
parents:
diff changeset
122 # for every <Iteration> tag
d375502056f1 Uploaded
devteam
parents:
diff changeset
123 if event == "end" and elem.tag == "Iteration":
d375502056f1 Uploaded
devteam
parents:
diff changeset
124 #Expecting either this, from BLAST 2.2.25+ using FASTA vs FASTA
d375502056f1 Uploaded
devteam
parents:
diff changeset
125 # <Iteration_query-ID>sp|Q9BS26|ERP44_HUMAN</Iteration_query-ID>
d375502056f1 Uploaded
devteam
parents:
diff changeset
126 # <Iteration_query-def>Endoplasmic reticulum resident protein 44 OS=Homo sapiens GN=ERP44 PE=1 SV=1</Iteration_query-def>
d375502056f1 Uploaded
devteam
parents:
diff changeset
127 # <Iteration_query-len>406</Iteration_query-len>
d375502056f1 Uploaded
devteam
parents:
diff changeset
128 # <Iteration_hits></Iteration_hits>
d375502056f1 Uploaded
devteam
parents:
diff changeset
129 #
d375502056f1 Uploaded
devteam
parents:
diff changeset
130 #Or, from BLAST 2.2.24+ run online
d375502056f1 Uploaded
devteam
parents:
diff changeset
131 # <Iteration_query-ID>Query_1</Iteration_query-ID>
d375502056f1 Uploaded
devteam
parents:
diff changeset
132 # <Iteration_query-def>Sample</Iteration_query-def>
d375502056f1 Uploaded
devteam
parents:
diff changeset
133 # <Iteration_query-len>516</Iteration_query-len>
d375502056f1 Uploaded
devteam
parents:
diff changeset
134 # <Iteration_hits>...
d375502056f1 Uploaded
devteam
parents:
diff changeset
135 qseqid = elem.findtext("Iteration_query-ID")
d375502056f1 Uploaded
devteam
parents:
diff changeset
136 if re_default_query_id.match(qseqid):
d375502056f1 Uploaded
devteam
parents:
diff changeset
137 #Place holder ID, take the first word of the query definition
d375502056f1 Uploaded
devteam
parents:
diff changeset
138 qseqid = elem.findtext("Iteration_query-def").split(None,1)[0]
d375502056f1 Uploaded
devteam
parents:
diff changeset
139 qlen = int(elem.findtext("Iteration_query-len"))
d375502056f1 Uploaded
devteam
parents:
diff changeset
140
d375502056f1 Uploaded
devteam
parents:
diff changeset
141 # for every <Hit> within <Iteration>
d375502056f1 Uploaded
devteam
parents:
diff changeset
142 for hit in elem.findall("Iteration_hits/Hit"):
d375502056f1 Uploaded
devteam
parents:
diff changeset
143 #Expecting either this,
d375502056f1 Uploaded
devteam
parents:
diff changeset
144 # <Hit_id>gi|3024260|sp|P56514.1|OPSD_BUFBU</Hit_id>
d375502056f1 Uploaded
devteam
parents:
diff changeset
145 # <Hit_def>RecName: Full=Rhodopsin</Hit_def>
d375502056f1 Uploaded
devteam
parents:
diff changeset
146 # <Hit_accession>P56514</Hit_accession>
d375502056f1 Uploaded
devteam
parents:
diff changeset
147 #or,
d375502056f1 Uploaded
devteam
parents:
diff changeset
148 # <Hit_id>Subject_1</Hit_id>
d375502056f1 Uploaded
devteam
parents:
diff changeset
149 # <Hit_def>gi|57163783|ref|NP_001009242.1| rhodopsin [Felis catus]</Hit_def>
d375502056f1 Uploaded
devteam
parents:
diff changeset
150 # <Hit_accession>Subject_1</Hit_accession>
d375502056f1 Uploaded
devteam
parents:
diff changeset
151 #
d375502056f1 Uploaded
devteam
parents:
diff changeset
152 #apparently depending on the parse_deflines switch
d375502056f1 Uploaded
devteam
parents:
diff changeset
153 sseqid = hit.findtext("Hit_id").split(None,1)[0]
d375502056f1 Uploaded
devteam
parents:
diff changeset
154 hit_def = sseqid + " " + hit.findtext("Hit_def")
d375502056f1 Uploaded
devteam
parents:
diff changeset
155 if re_default_subject_id.match(sseqid) \
d375502056f1 Uploaded
devteam
parents:
diff changeset
156 and sseqid == hit.findtext("Hit_accession"):
d375502056f1 Uploaded
devteam
parents:
diff changeset
157 #Place holder ID, take the first word of the subject definition
d375502056f1 Uploaded
devteam
parents:
diff changeset
158 hit_def = hit.findtext("Hit_def")
d375502056f1 Uploaded
devteam
parents:
diff changeset
159 sseqid = hit_def.split(None,1)[0]
d375502056f1 Uploaded
devteam
parents:
diff changeset
160 # for every <Hsp> within <Hit>
d375502056f1 Uploaded
devteam
parents:
diff changeset
161 for hsp in hit.findall("Hit_hsps/Hsp"):
d375502056f1 Uploaded
devteam
parents:
diff changeset
162 nident = hsp.findtext("Hsp_identity")
d375502056f1 Uploaded
devteam
parents:
diff changeset
163 length = hsp.findtext("Hsp_align-len")
d375502056f1 Uploaded
devteam
parents:
diff changeset
164 pident = "%0.2f" % (100*float(nident)/float(length))
d375502056f1 Uploaded
devteam
parents:
diff changeset
165
d375502056f1 Uploaded
devteam
parents:
diff changeset
166 q_seq = hsp.findtext("Hsp_qseq")
d375502056f1 Uploaded
devteam
parents:
diff changeset
167 h_seq = hsp.findtext("Hsp_hseq")
d375502056f1 Uploaded
devteam
parents:
diff changeset
168 m_seq = hsp.findtext("Hsp_midline")
d375502056f1 Uploaded
devteam
parents:
diff changeset
169 assert len(q_seq) == len(h_seq) == len(m_seq) == int(length)
d375502056f1 Uploaded
devteam
parents:
diff changeset
170 gapopen = str(len(q_seq.replace('-', ' ').split())-1 + \
d375502056f1 Uploaded
devteam
parents:
diff changeset
171 len(h_seq.replace('-', ' ').split())-1)
d375502056f1 Uploaded
devteam
parents:
diff changeset
172
d375502056f1 Uploaded
devteam
parents:
diff changeset
173 mismatch = m_seq.count(' ') + m_seq.count('+') \
d375502056f1 Uploaded
devteam
parents:
diff changeset
174 - q_seq.count('-') - h_seq.count('-')
d375502056f1 Uploaded
devteam
parents:
diff changeset
175 #TODO - Remove this alternative mismatch calculation and test
d375502056f1 Uploaded
devteam
parents:
diff changeset
176 #once satisifed there are no problems
d375502056f1 Uploaded
devteam
parents:
diff changeset
177 expected_mismatch = len(q_seq) \
d375502056f1 Uploaded
devteam
parents:
diff changeset
178 - sum(1 for q,h in zip(q_seq, h_seq) \
d375502056f1 Uploaded
devteam
parents:
diff changeset
179 if q == h or q == "-" or h == "-")
d375502056f1 Uploaded
devteam
parents:
diff changeset
180 xx = sum(1 for q,h in zip(q_seq, h_seq) if q=="X" and h=="X")
d375502056f1 Uploaded
devteam
parents:
diff changeset
181 if not (expected_mismatch - q_seq.count("X") <= int(mismatch) <= expected_mismatch + xx):
d375502056f1 Uploaded
devteam
parents:
diff changeset
182 stop_err("%s vs %s mismatches, expected %i <= %i <= %i" \
d375502056f1 Uploaded
devteam
parents:
diff changeset
183 % (qseqid, sseqid, expected_mismatch - q_seq.count("X"),
d375502056f1 Uploaded
devteam
parents:
diff changeset
184 int(mismatch), expected_mismatch))
d375502056f1 Uploaded
devteam
parents:
diff changeset
185
d375502056f1 Uploaded
devteam
parents:
diff changeset
186 #TODO - Remove this alternative identity calculation and test
d375502056f1 Uploaded
devteam
parents:
diff changeset
187 #once satisifed there are no problems
d375502056f1 Uploaded
devteam
parents:
diff changeset
188 expected_identity = sum(1 for q,h in zip(q_seq, h_seq) if q == h)
d375502056f1 Uploaded
devteam
parents:
diff changeset
189 if not (expected_identity - xx <= int(nident) <= expected_identity + q_seq.count("X")):
d375502056f1 Uploaded
devteam
parents:
diff changeset
190 stop_err("%s vs %s identities, expected %i <= %i <= %i" \
d375502056f1 Uploaded
devteam
parents:
diff changeset
191 % (qseqid, sseqid, expected_identity, int(nident),
d375502056f1 Uploaded
devteam
parents:
diff changeset
192 expected_identity + q_seq.count("X")))
d375502056f1 Uploaded
devteam
parents:
diff changeset
193
d375502056f1 Uploaded
devteam
parents:
diff changeset
194
d375502056f1 Uploaded
devteam
parents:
diff changeset
195 evalue = hsp.findtext("Hsp_evalue")
d375502056f1 Uploaded
devteam
parents:
diff changeset
196 if evalue == "0":
d375502056f1 Uploaded
devteam
parents:
diff changeset
197 evalue = "0.0"
d375502056f1 Uploaded
devteam
parents:
diff changeset
198 else:
d375502056f1 Uploaded
devteam
parents:
diff changeset
199 evalue = "%0.0e" % float(evalue)
d375502056f1 Uploaded
devteam
parents:
diff changeset
200
d375502056f1 Uploaded
devteam
parents:
diff changeset
201 bitscore = float(hsp.findtext("Hsp_bit-score"))
d375502056f1 Uploaded
devteam
parents:
diff changeset
202 if bitscore < 100:
d375502056f1 Uploaded
devteam
parents:
diff changeset
203 #Seems to show one decimal place for lower scores
d375502056f1 Uploaded
devteam
parents:
diff changeset
204 bitscore = "%0.1f" % bitscore
d375502056f1 Uploaded
devteam
parents:
diff changeset
205 else:
d375502056f1 Uploaded
devteam
parents:
diff changeset
206 #Note BLAST does not round to nearest int, it truncates
d375502056f1 Uploaded
devteam
parents:
diff changeset
207 bitscore = "%i" % bitscore
d375502056f1 Uploaded
devteam
parents:
diff changeset
208
d375502056f1 Uploaded
devteam
parents:
diff changeset
209 values = [qseqid,
d375502056f1 Uploaded
devteam
parents:
diff changeset
210 sseqid,
d375502056f1 Uploaded
devteam
parents:
diff changeset
211 pident,
d375502056f1 Uploaded
devteam
parents:
diff changeset
212 length, #hsp.findtext("Hsp_align-len")
d375502056f1 Uploaded
devteam
parents:
diff changeset
213 str(mismatch),
d375502056f1 Uploaded
devteam
parents:
diff changeset
214 gapopen,
d375502056f1 Uploaded
devteam
parents:
diff changeset
215 hsp.findtext("Hsp_query-from"), #qstart,
d375502056f1 Uploaded
devteam
parents:
diff changeset
216 hsp.findtext("Hsp_query-to"), #qend,
d375502056f1 Uploaded
devteam
parents:
diff changeset
217 hsp.findtext("Hsp_hit-from"), #sstart,
d375502056f1 Uploaded
devteam
parents:
diff changeset
218 hsp.findtext("Hsp_hit-to"), #send,
d375502056f1 Uploaded
devteam
parents:
diff changeset
219 evalue, #hsp.findtext("Hsp_evalue") in scientific notation
d375502056f1 Uploaded
devteam
parents:
diff changeset
220 bitscore, #hsp.findtext("Hsp_bit-score") rounded
d375502056f1 Uploaded
devteam
parents:
diff changeset
221 ]
d375502056f1 Uploaded
devteam
parents:
diff changeset
222
d375502056f1 Uploaded
devteam
parents:
diff changeset
223 if extended:
d375502056f1 Uploaded
devteam
parents:
diff changeset
224 sallseqid = ";".join(name.split(None,1)[0] for name in hit_def.split(">"))
d375502056f1 Uploaded
devteam
parents:
diff changeset
225 #print hit_def, "-->", sallseqid
d375502056f1 Uploaded
devteam
parents:
diff changeset
226 positive = hsp.findtext("Hsp_positive")
d375502056f1 Uploaded
devteam
parents:
diff changeset
227 ppos = "%0.2f" % (100*float(positive)/float(length))
d375502056f1 Uploaded
devteam
parents:
diff changeset
228 qframe = hsp.findtext("Hsp_query-frame")
d375502056f1 Uploaded
devteam
parents:
diff changeset
229 sframe = hsp.findtext("Hsp_hit-frame")
d375502056f1 Uploaded
devteam
parents:
diff changeset
230 if blast_program == "blastp":
d375502056f1 Uploaded
devteam
parents:
diff changeset
231 #Probably a bug in BLASTP that they use 0 or 1 depending on format
d375502056f1 Uploaded
devteam
parents:
diff changeset
232 if qframe == "0": qframe = "1"
d375502056f1 Uploaded
devteam
parents:
diff changeset
233 if sframe == "0": sframe = "1"
d375502056f1 Uploaded
devteam
parents:
diff changeset
234 slen = int(hit.findtext("Hit_len"))
d375502056f1 Uploaded
devteam
parents:
diff changeset
235 values.extend([sallseqid,
d375502056f1 Uploaded
devteam
parents:
diff changeset
236 hsp.findtext("Hsp_score"), #score,
d375502056f1 Uploaded
devteam
parents:
diff changeset
237 nident,
d375502056f1 Uploaded
devteam
parents:
diff changeset
238 positive,
d375502056f1 Uploaded
devteam
parents:
diff changeset
239 hsp.findtext("Hsp_gaps"), #gaps,
d375502056f1 Uploaded
devteam
parents:
diff changeset
240 ppos,
d375502056f1 Uploaded
devteam
parents:
diff changeset
241 qframe,
d375502056f1 Uploaded
devteam
parents:
diff changeset
242 sframe,
d375502056f1 Uploaded
devteam
parents:
diff changeset
243 #NOTE - for blastp, XML shows original seq, tabular uses XXX masking
d375502056f1 Uploaded
devteam
parents:
diff changeset
244 q_seq,
d375502056f1 Uploaded
devteam
parents:
diff changeset
245 h_seq,
d375502056f1 Uploaded
devteam
parents:
diff changeset
246 str(qlen),
d375502056f1 Uploaded
devteam
parents:
diff changeset
247 str(slen),
d375502056f1 Uploaded
devteam
parents:
diff changeset
248 ])
d375502056f1 Uploaded
devteam
parents:
diff changeset
249 #print "\t".join(values)
d375502056f1 Uploaded
devteam
parents:
diff changeset
250 outfile.write("\t".join(values) + "\n")
d375502056f1 Uploaded
devteam
parents:
diff changeset
251 # prevents ElementTree from growing large datastructure
d375502056f1 Uploaded
devteam
parents:
diff changeset
252 root.clear()
d375502056f1 Uploaded
devteam
parents:
diff changeset
253 elem.clear()
d375502056f1 Uploaded
devteam
parents:
diff changeset
254 outfile.close()