Mercurial > repos > galaxyp > blastxml_to_tabular_selectable
annotate blastxml_to_tabular_selectable.py @ 1:5da5dcc5e13a default tip
Update from 24 to 25 columns and update versions to 1.0.0
| author | Jim Johnson <jj@umn.edu> |
|---|---|
| date | Wed, 08 Oct 2014 18:57:39 -0500 |
| parents | 2bd0cbccb3c6 |
| children |
| rev | line source |
|---|---|
| 0 | 1 #!/usr/bin/env python |
| 2 """Convert a BLAST XML file to 12 column tabular output | |
| 3 | |
| 4 Takes three command line options, input BLAST XML filename, output tabular | |
| 5 BLAST filename, output format (std for standard 12 columns, or ext for the | |
| 6 extended 25 columns offered in the BLAST+ wrappers). | |
| 7 | |
| 8 The 12 columns output are 'qseqid sseqid pident length mismatch gapopen qstart | |
| 9 qend sstart send evalue bitscore' or 'std' at the BLAST+ command line, which | |
| 10 mean: | |
| 11 | |
| 12 ====== ========= ============================================ | |
| 13 Column NCBI name Description | |
| 14 ------ --------- -------------------------------------------- | |
| 15 1 qseqid Query Seq-id (ID of your sequence) | |
| 16 2 sseqid Subject Seq-id (ID of the database hit) | |
| 17 3 pident Percentage of identical matches | |
| 18 4 length Alignment length | |
| 19 5 mismatch Number of mismatches | |
| 20 6 gapopen Number of gap openings | |
| 21 7 qstart Start of alignment in query | |
| 22 8 qend End of alignment in query | |
| 23 9 sstart Start of alignment in subject (database hit) | |
| 24 10 send End of alignment in subject (database hit) | |
| 25 11 evalue Expectation value (E-value) | |
| 26 12 bitscore Bit score | |
| 27 ====== ========= ============================================ | |
| 28 | |
| 29 The additional columns offered in the Galaxy BLAST+ wrappers are: | |
| 30 | |
| 31 ====== ============= =========================================== | |
| 32 Column NCBI name Description | |
| 33 ------ ------------- ------------------------------------------- | |
| 34 13 sallseqid All subject Seq-id(s), separated by a ';' | |
| 35 14 score Raw score | |
| 36 15 nident Number of identical matches | |
| 37 16 positive Number of positive-scoring matches | |
| 38 17 gaps Total number of gaps | |
| 39 18 ppos Percentage of positive-scoring matches | |
| 40 19 qframe Query frame | |
| 41 20 sframe Subject frame | |
| 42 21 qseq Aligned part of query sequence | |
| 43 22 sseq Aligned part of subject sequence | |
| 44 23 qlen Query sequence length | |
| 45 24 slen Subject sequence length | |
| 46 25 salltitles All subject titles, separated by '<>' | |
| 47 ====== ============= =========================================== | |
| 48 | |
| 49 Most of these fields are given explicitly in the XML file, others some like | |
| 50 the percentage identity and the number of gap openings must be calculated. | |
| 51 | |
| 52 Be aware that the sequence in the extended tabular output or XML direct from | |
| 53 BLAST+ may or may not use XXXX masking on regions of low complexity. This | |
| 54 can throw the off the calculation of percentage identity and gap openings. | |
| 55 [In fact, both BLAST 2.2.24+ and 2.2.25+ have a subtle bug in this regard, | |
| 56 with these numbers changing depending on whether or not the low complexity | |
| 57 filter is used.] | |
| 58 | |
| 59 This script attempts to produce identical output to what BLAST+ would have done. | |
| 60 However, check this with "diff -b ..." since BLAST+ sometimes includes an extra | |
| 61 space character (probably a bug). | |
| 62 """ | |
| 63 import sys | |
| 64 import re | |
| 65 import os | |
| 66 from optparse import OptionParser | |
| 67 | |
| 68 if "-v" in sys.argv or "--version" in sys.argv: | |
|
1
5da5dcc5e13a
Update from 24 to 25 columns and update versions to 1.0.0
Jim Johnson <jj@umn.edu>
parents:
0
diff
changeset
|
69 print "v1.0.0" |
| 0 | 70 sys.exit(0) |
| 71 | |
| 72 if sys.version_info[:2] >= ( 2, 5 ): | |
| 73 try: | |
| 74 from xml.etree import cElementTree as ElementTree | |
| 75 except ImportError: | |
| 76 from xml.etree import ElementTree as ElementTree | |
| 77 else: | |
| 78 from galaxy import eggs | |
| 79 import pkg_resources; pkg_resources.require( "elementtree" ) | |
| 80 from elementtree import ElementTree | |
| 81 | |
| 82 def stop_err( msg ): | |
| 83 sys.stderr.write("%s\n" % msg) | |
| 84 sys.exit(1) | |
| 85 | |
| 86 usage = "usage: %prog [options] blastxml[,...]" | |
| 87 parser = OptionParser(usage=usage) | |
| 88 parser.add_option('-o','--output', dest='output', default = None, help='output file path', metavar="FILE") | |
| 89 parser.add_option("-c", "--columns", dest="columns", default='std', help="[std|ext|colname[,colname,...]] std: 12 column, ext: 25 column, or user specified columns") | |
| 90 parser.add_option("-a", "--allqueries", action="store_true", dest="allqueries", default=False, help="ouptut row for each query, including those with no hits") | |
| 91 parser.add_option("-d", "--qdef", action="store_true", dest="qdef", default=False, help="use Iteration_query-def for qseqid") | |
| 92 parser.add_option('-u', '--unmatched', dest='unmatched', default = None, help='unmatched queries file path', metavar="FILE") | |
| 93 parser.add_option('-m', '--maxhits', dest='maxhits', type="int", default = 1, help='maximum number of HITs for output for a Query') | |
| 94 parser.add_option('-p', '--maxhsps', dest='maxhsps', type="int", default = 1, help='maximum number of HSPs for output for a Hit') | |
| 95 (options, args) = parser.parse_args() | |
| 96 | |
| 97 colnames = 'qseqid,sseqid,pident,length,mismatch,gapopen,qstart,qend,sstart,send,evalue,bitscore,sallseqid,score,nident,positive,gaps,ppos,qframe,sframe,qseq,sseq,qlen,slen,salltitles'.split(',') | |
| 98 | |
| 99 if len(args) < 1: | |
| 100 parser.error("no BLASTXML input files") | |
| 101 | |
| 102 if options.columns == "std": | |
| 103 out_fmt = 'std' | |
| 104 extended = False | |
| 105 columns = colnames[0:12] | |
| 106 elif options.columns == "x22": | |
| 107 stop_err("Format argument x22 has been replaced with ext (extended 25 columns)") | |
| 108 elif options.columns == "ext": | |
| 109 out_fmt = 'ext' | |
| 110 extended = True | |
| 111 columns = colnames[0:25] | |
| 112 else: | |
| 113 out_fmt = 'cols' | |
| 114 extended = True | |
| 115 sel_cols = options.columns.split(',') | |
| 116 columns = [] | |
| 117 for col in sel_cols: | |
| 118 if col in colnames: | |
| 119 columns.append(col) | |
| 120 | |
| 121 for in_file in args: | |
| 122 if not os.path.isfile(in_file): | |
| 123 stop_err("Input BLAST XML file not found: %s" % in_file) | |
| 124 | |
| 125 unhitfile = open(options.unmatched,'w') if options.unmatched else None | |
| 126 | |
| 127 re_default_query_id = re.compile("^Query_\d+$") | |
| 128 assert re_default_query_id.match("Query_101") | |
| 129 assert not re_default_query_id.match("Query_101a") | |
| 130 assert not re_default_query_id.match("MyQuery_101") | |
| 131 re_default_subject_id = re.compile("^Subject_\d+$") | |
| 132 assert re_default_subject_id.match("Subject_1") | |
| 133 assert not re_default_subject_id.match("Subject_") | |
| 134 assert not re_default_subject_id.match("Subject_12a") | |
| 135 assert not re_default_subject_id.match("TheSubject_1") | |
| 136 | |
| 137 header = "#%s\n" % "\t".join(columns) | |
| 138 outfile = open(options.output, 'w') if options.output else sys.stdout | |
| 139 outfile.write(header) | |
| 140 | |
| 141 def handle_event(event, elem): | |
| 142 blast_program = None | |
| 143 if event == "end" and elem.tag == "BlastOutput_program": | |
| 144 blast_program = elem.text | |
| 145 # for every <Iteration> tag | |
| 146 if event == "end" and elem.tag == "Iteration": | |
| 147 #Expecting either this, from BLAST 2.2.25+ using FASTA vs FASTA | |
| 148 # <Iteration_query-ID>sp|Q9BS26|ERP44_HUMAN</Iteration_query-ID> | |
| 149 # <Iteration_query-def>Endoplasmic reticulum resident protein 44 OS=Homo sapiens GN=ERP44 PE=1 SV=1</Iteration_query-def> | |
| 150 # <Iteration_query-len>406</Iteration_query-len> | |
| 151 # <Iteration_hits></Iteration_hits> | |
| 152 # | |
| 153 #Or, from BLAST 2.2.24+ run online | |
| 154 # <Iteration_query-ID>Query_1</Iteration_query-ID> | |
| 155 # <Iteration_query-def>Sample</Iteration_query-def> | |
| 156 # <Iteration_query-len>516</Iteration_query-len> | |
| 157 # <Iteration_hits>... | |
| 158 qseqid = elem.findtext("Iteration_query-ID") | |
| 159 if options.qdef or re_default_query_id.match(qseqid): | |
| 160 #Place holder ID, take the first word of the query definition | |
| 161 qseqid = elem.findtext("Iteration_query-def").split(None,1)[0] | |
| 162 qlen = int(elem.findtext("Iteration_query-len")) | |
| 163 # If no hits in this iteration | |
| 164 if not elem.find("Iteration_hits/Hit"): | |
| 165 if unhitfile: | |
| 166 unhitfile.write("%s\n" % qseqid) | |
| 167 if options.allqueries: | |
| 168 if columns and len(columns) > 0: | |
| 169 v = [] | |
| 170 for name in columns: | |
| 171 if name == 'qseqid': | |
| 172 v.append(qseqid) | |
| 173 elif name == 'qlen': | |
| 174 v.append(str(qlen)) | |
| 175 else: | |
| 176 v.append('') | |
| 177 outfile.write("\t".join(v) + "\n") | |
| 178 # for every <Hit> within <Iteration> | |
| 179 for hitn, hit in enumerate(elem.findall("Iteration_hits/Hit")): | |
| 180 if hitn >= options.maxhits: | |
| 181 break | |
| 182 #Expecting either this, | |
| 183 # <Hit_id>gi|3024260|sp|P56514.1|OPSD_BUFBU</Hit_id> | |
| 184 # <Hit_def>RecName: Full=Rhodopsin</Hit_def> | |
| 185 # <Hit_accession>P56514</Hit_accession> | |
| 186 #or, | |
| 187 # <Hit_id>Subject_1</Hit_id> | |
| 188 # <Hit_def>gi|57163783|ref|NP_001009242.1| rhodopsin [Felis catus]</Hit_def> | |
| 189 # <Hit_accession>Subject_1</Hit_accession> | |
| 190 # | |
| 191 #apparently depending on the parse_deflines switch | |
| 192 # | |
| 193 #Or, with a local database not using -parse_seqids can get this, | |
| 194 # <Hit_id>gnl|BL_ORD_ID|2</Hit_id> | |
| 195 # <Hit_def>chrIII gi|240255695|ref|NC_003074.8| Arabidopsis thaliana chromosome 3, complete sequence</Hit_def> | |
| 196 # <Hit_accession>2</Hit_accession> | |
| 197 sseqid = hit.findtext("Hit_id").split(None,1)[0] | |
| 198 hit_def = sseqid + " " + hit.findtext("Hit_def") | |
| 199 if re_default_subject_id.match(sseqid) \ | |
| 200 and sseqid == hit.findtext("Hit_accession"): | |
| 201 #Place holder ID, take the first word of the subject definition | |
| 202 hit_def = hit.findtext("Hit_def") | |
| 203 sseqid = hit_def.split(None,1)[0] | |
| 204 # for every <Hsp> within <Hit> | |
| 205 for hspn,hsp in enumerate(hit.findall("Hit_hsps/Hsp")): | |
| 206 if hspn >= options.maxhsps: | |
| 207 break | |
| 208 nident = hsp.findtext("Hsp_identity") | |
| 209 length = hsp.findtext("Hsp_align-len") | |
| 210 pident = "%0.2f" % (100*float(nident)/float(length)) | |
| 211 | |
| 212 q_seq = hsp.findtext("Hsp_qseq") | |
| 213 h_seq = hsp.findtext("Hsp_hseq") | |
| 214 m_seq = hsp.findtext("Hsp_midline") | |
| 215 assert len(q_seq) == len(h_seq) == len(m_seq) == int(length) | |
| 216 gapopen = str(len(q_seq.replace('-', ' ').split())-1 + \ | |
| 217 len(h_seq.replace('-', ' ').split())-1) | |
| 218 | |
| 219 mismatch = m_seq.count(' ') + m_seq.count('+') \ | |
| 220 - q_seq.count('-') - h_seq.count('-') | |
| 221 #TODO - Remove this alternative mismatch calculation and test | |
| 222 #once satisifed there are no problems | |
| 223 expected_mismatch = len(q_seq) \ | |
| 224 - sum(1 for q,h in zip(q_seq, h_seq) \ | |
| 225 if q == h or q == "-" or h == "-") | |
| 226 xx = sum(1 for q,h in zip(q_seq, h_seq) if q=="X" and h=="X") | |
| 227 if not (expected_mismatch - q_seq.count("X") <= int(mismatch) <= expected_mismatch + xx): | |
| 228 stop_err("%s vs %s mismatches, expected %i <= %i <= %i" \ | |
| 229 % (qseqid, sseqid, expected_mismatch - q_seq.count("X"), | |
| 230 int(mismatch), expected_mismatch)) | |
| 231 | |
| 232 #TODO - Remove this alternative identity calculation and test | |
| 233 #once satisifed there are no problems | |
| 234 expected_identity = sum(1 for q,h in zip(q_seq, h_seq) if q == h) | |
| 235 if not (expected_identity - xx <= int(nident) <= expected_identity + q_seq.count("X")): | |
| 236 stop_err("%s vs %s identities, expected %i <= %i <= %i" \ | |
| 237 % (qseqid, sseqid, expected_identity, int(nident), | |
| 238 expected_identity + q_seq.count("X"))) | |
| 239 | |
| 240 | |
| 241 evalue = hsp.findtext("Hsp_evalue") | |
| 242 if evalue == "0": | |
| 243 evalue = "0.0" | |
| 244 else: | |
| 245 evalue = "%0.0e" % float(evalue) | |
| 246 | |
| 247 bitscore = float(hsp.findtext("Hsp_bit-score")) | |
| 248 if bitscore < 100: | |
| 249 #Seems to show one decimal place for lower scores | |
| 250 bitscore = "%0.1f" % bitscore | |
| 251 else: | |
| 252 #Note BLAST does not round to nearest int, it truncates | |
| 253 bitscore = "%i" % bitscore | |
| 254 | |
| 255 values = [qseqid, | |
| 256 sseqid, | |
| 257 pident, | |
| 258 length, #hsp.findtext("Hsp_align-len") | |
| 259 str(mismatch), | |
| 260 gapopen, | |
| 261 hsp.findtext("Hsp_query-from"), #qstart, | |
| 262 hsp.findtext("Hsp_query-to"), #qend, | |
| 263 hsp.findtext("Hsp_hit-from"), #sstart, | |
| 264 hsp.findtext("Hsp_hit-to"), #send, | |
| 265 evalue, #hsp.findtext("Hsp_evalue") in scientific notation | |
| 266 bitscore, #hsp.findtext("Hsp_bit-score") rounded | |
| 267 ] | |
| 268 | |
| 269 if extended: | |
| 270 sallseqid = ";".join(name.split(None,1)[0] for name in hit_def.split(">")) | |
| 271 salltitles = "<>".join(name.split(None,1)[1] for name in hit_def.split(" >")) | |
| 272 #print hit_def, "-->", sallseqid | |
| 273 positive = hsp.findtext("Hsp_positive") | |
| 274 ppos = "%0.2f" % (100*float(positive)/float(length)) | |
| 275 qframe = hsp.findtext("Hsp_query-frame") | |
| 276 sframe = hsp.findtext("Hsp_hit-frame") | |
| 277 if blast_program == "blastp": | |
| 278 #Probably a bug in BLASTP that they use 0 or 1 depending on format | |
| 279 if qframe == "0": qframe = "1" | |
| 280 if sframe == "0": sframe = "1" | |
| 281 slen = int(hit.findtext("Hit_len")) | |
| 282 values.extend([sallseqid, | |
| 283 hsp.findtext("Hsp_score"), #score, | |
| 284 nident, | |
| 285 positive, | |
| 286 hsp.findtext("Hsp_gaps"), #gaps, | |
| 287 ppos, | |
| 288 qframe, | |
| 289 sframe, | |
| 290 #NOTE - for blastp, XML shows original seq, tabular uses XXX masking | |
| 291 q_seq, | |
| 292 h_seq, | |
| 293 str(qlen), | |
| 294 str(slen), | |
| 295 salltitles, | |
| 296 ]) | |
| 297 if out_fmt == 'cols': | |
| 298 if columns and len(columns) > 0: | |
| 299 v = [] | |
| 300 for name in columns: | |
| 301 v.append(values[colnames.index(name)]) | |
| 302 values = v | |
| 303 #print "\t".join(values) | |
| 304 outfile.write("\t".join(values) + "\n") | |
| 305 # prevents ElementTree from growing large datastructure | |
| 306 root.clear() | |
| 307 elem.clear() | |
| 308 | |
| 309 | |
| 310 for in_file in args: | |
| 311 # get an iterable | |
| 312 try: | |
| 313 context = ElementTree.iterparse(in_file, events=("start", "end")) | |
| 314 except: | |
| 315 stop_err("Invalid data format.") | |
| 316 # turn it into an iterator | |
| 317 context = iter(context) | |
| 318 # get the root element | |
| 319 try: | |
| 320 event, root = context.next() | |
| 321 except: | |
| 322 stop_err( "Invalid data format." ) | |
| 323 for event, elem in context: | |
| 324 handle_event(event, elem) | |
| 325 | |
| 326 if unhitfile: | |
| 327 unhitfile.close() | |
| 328 if options.output: | |
| 329 outfile.close() |