#!/usr/bin/python3 ################################################################ # copyright (c) 2017,2018 by William R. Pearson and The Rector & # Visitors of the University of Virginia */ ################################################################ # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. # You may obtain a copy of the License at # # http://www.apache.org/licenses/LICENSE-2.0 # # Unless required by applicable law or agreed to in writing, # software distributed under this License is distributed on an "AS # IS" BASIS, WITHOUT WRRANTIES OR CONDITIONS OF ANY KIND, either # express or implied. See the License for the specific language # governing permissions and limitations under the License. ################################################################ ################################################################ # annot_blast_btop4.py --query query.file --ann_script ann_pfam_www.pl --include_doms blast_tab_btop_file ################################################################ # annot_blast_btop4.py associates domain annotation information and # subalignment scores with a blast tabular (-outfmt 6 or -outfmt 7) # file that contains the raw score and the BTOP alignment encoding # This file can be generated from "blastp/n" or "blast_formatter" # using the command: # blast_formatter -archive blast_output.asn -outfmt '7 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore score btop' > blast_output.tab_annot # # If the BTOP field or query_file is not available, the script # produces domain content without sub-alignment scores. ################################################################ ## 2-Dec-2019 # added --have_qslen, --raw_score/--no_raw_score # made more robust to multiple HSPs when using --ann_file # ################################################################ ## 4-Nov-2018 # add --include_doms, which adds a new field with the coordinates of # the domains in the protein (independent of alignment) # ################################################################ ## 21-July-2018 # include sequence length (actually alignment end) to produce NODOM's (no NODOM's without length). # ################################################################ ## 13-Jan-2017 # modified to provide query/subject coordinates and identities if no # query sequence -- does not decrement for reverse-complement fastx/blastx DNA ################################################################ ## 16-Nov-2015 # modify to allow multi-query blast searches ################################################################ ## 19-Dec-2015 # add -q_annot_script to annotate query sequence # import argparse import fileinput import sys import re import shutil import subprocess from math import log # read lines of the form: # gi|121694|sp|P20432.1|GSTT1_DROME gi|121694|sp|P20432|GSTT1_DROME 100.00 209 0 0 1 209 1 209 6e-156 433 1113 209 # gi|121694|sp|P20432.1|GSTT1_DROME gi|1170090|sp|P04907|GSTF3_MAIZE 26.77 198 123 7 4 185 6 197 2e-08 51.2 121 FL1YG ... 1NKRA1YW1 # gi|121694|sp|P20432.1|GSTT1_DROME gi|81174731|sp|P0ACA5|SSPA_ECO57 39.66 58 32 2 43 100 49 103 8e-06 43.9 102 EDFLLI ... V-I-NEQS3FM # gi|121694|sp|P20432.1|GSTT1_DROME gi|121695|sp|P12653|GSTF1_MAIZE 27.62 181 107 7 32 203 34 199 9e-05 40.8 94 LI1LF ... N-1AS1CLLM1 # and report the domain content ala -m 8CC def init_blosum62(): # ncbi_blaa -- list of amino acids ncbi_blaa = "A R N D C Q E G H I L K M F P S T W Y V B Z X *".split(' ') # blosum62: 2D dict of scoring matrix values blosum62 = {} blosum62['A'] = dict(zip(ncbi_blaa,[ 4,-1,-2,-2, 0,-1,-1, 0,-2,-1,-1,-1,-1,-2,-1, 1, 0,-3,-2, 0,-2,-1, 0,-4])) blosum62['R'] = dict(zip(ncbi_blaa,[-1, 5, 0,-2,-3, 1, 0,-2, 0,-3,-2, 2,-1,-3,-2,-1,-1,-3,-2,-3,-1, 0,-1,-4])) blosum62['N'] = dict(zip(ncbi_blaa,[-2, 0, 6, 1,-3, 0, 0, 0, 1,-3,-3, 0,-2,-3,-2, 1, 0,-4,-2,-3, 3, 0,-1,-4])) blosum62['D'] = dict(zip(ncbi_blaa,[-2,-2, 1, 6,-3, 0, 2,-1,-1,-3,-4,-1,-3,-3,-1,0,-1,-4,-3,-3,4,1,-1,-4])) blosum62['C'] = dict(zip(ncbi_blaa,[ 0,-3,-3,-3, 9,-3,-4,-3,-3,-1,-1,-3,-1,-2,-3,-1,-1,-2,-2,-1,-3,-3,-2,-4])) blosum62['Q'] = dict(zip(ncbi_blaa,[-1, 1, 0, 0,-3, 5, 2,-2, 0,-3,-2,1,0,-3,-1,0,-1,-2,-1,-2,0,3,-1,-4])) blosum62['E'] = dict(zip(ncbi_blaa,[-1, 0, 0, 2,-4, 2, 5,-2, 0,-3,-3,1,-2,-3,-1,0,-1,-3,-2,-2,1,4,-1,-4])) blosum62['G'] = dict(zip(ncbi_blaa,[ 0,-2, 0,-1,-3,-2,-2, 6,-2,-4,-4,-2,-3,-3,-2,0,-2,-2,-3,-3,-1,-2,-1,-4])) blosum62['H'] = dict(zip(ncbi_blaa,[-2, 0, 1,-1,-3, 0, 0,-2, 8,-3,-3,-1,-2,-1,-2,-1,-2,-2,2,-3,0,0,-1,-4])) blosum62['I'] = dict(zip(ncbi_blaa,[-1,-3,-3,-3,-1,-3,-3,-4,-3,4,2,-3,1,0,-3,-2,-1,-3,-1,3,-3,-3,-1,-4])) blosum62['L'] = dict(zip(ncbi_blaa,[-1,-2,-3,-4,-1,-2,-3,-4,-3,2,4,-2,2,0,-3,-2,-1,-2,-1,1,-4,-3,-1,-4])) blosum62['K'] = dict(zip(ncbi_blaa,[-1, 2, 0,-1,-3, 1, 1,-2,-1,-3,-2,5,-1,-3,-1,0,-1,-3,-2,-2,0,1,-1,-4])) blosum62['M'] = dict(zip(ncbi_blaa,[-1,-1,-2,-3,-1, 0,-2,-3,-2,1,2,-1,5,0,-2,-1,-1,-1,-1,1,-3,-1,-1,-4])) blosum62['F'] = dict(zip(ncbi_blaa,[-2,-3,-3,-3,-2,-3,-3,-3,-1,0,0,-3,0,6,-4,-2,-2,1,3,-1,-3,-3,-1,-4])) blosum62['P'] = dict(zip(ncbi_blaa,[-1,-2,-2,-1,-3,-1,-1,-2,-2,-3,-3,-1,-2,-4,7,-1,-1,-4,-3,-2,-2,-1,-2,-4])) blosum62['S'] = dict(zip(ncbi_blaa,[ 1,-1, 1, 0,-1, 0, 0, 0,-1,-2,-2,0,-1,-2,-1,4,1,-3,-2,-2,0,0,0,-4])) blosum62['T'] = dict(zip(ncbi_blaa,[ 0,-1, 0,-1,-1,-1,-1,-2,-2,-1,-1,-1,-1,-2,-1,1,5,-2,-2,0,-1,-1,0,-4])) blosum62['W'] = dict(zip(ncbi_blaa,[-3 -3,-4,-4,-2,-2,-3,-2,-2,-3,-2,-3,-1,1,-4,-3,-2,11,2,-3,-4,-3,-2,-4])) blosum62['Y'] = dict(zip(ncbi_blaa,[-2,-2,-2,-3,-2,-1,-2,-3, 2,-1,-1,-2,-1,3,-3,-2,-2,2,7,-1,-3,-2,-1,-4])) blosum62['V'] = dict(zip(ncbi_blaa,[ 0,-3,-3,-3,-1,-2,-2,-3,-3,3,1,-2,1,-1,-2,-2,0,-3,-1,4,-3,-2,-1,-4])) blosum62['B'] = dict(zip(ncbi_blaa,[-2,-1, 3, 4,-3, 0, 1,-1, 0,-3,-4,0,-3,-3,-2,0,-1,-4,-3,-3, 4, 1,-1,-4])) blosum62['Z'] = dict(zip(ncbi_blaa,[-1, 0, 0, 1,-3, 3, 4,-2, 0,-3,-3,1,-1,-3,-1,0,-1,-3,-2,-2,1,4,-1,-4])) blosum62['X'] = dict(zip(ncbi_blaa,[ 0,-1,-1,-1,-2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-2,0,0,-2,-1,-1,-1,-1,-1,-4])) blosum62['*'] = dict(zip(ncbi_blaa,[-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,1])) if (len(blosum62.keys()) != len(ncbi_blaa)): sys.stderr.write(" blosum62 length mismatch %d != %d\n" %(len(blosum62), len(ncbi_blaa))) print(' '.join(ncbi_blaa),file=sys.stderr) print(' '.join(blosum62.keys()),file=sys.stderr) exit(1) blosum62_diag = {x:blosum62[x][x] for x in ncbi_blaa} return (blosum62, blosum62_diag, -11, -1) ################ # read_annots (\@hit_list) # input: hit_entry['s_seq_id, etc'], target # output: modified $hit_entry['domains'] # modified $hit_entry['sites'] # # extend to make robust to multiple hits on the same subject def read_annots(Reader): target_set = {} current_domain = "" hit_ix = 0 seq_domains = [] seq_sites = [] subj_domains = {} for line in Reader: if (line[0]=='='): continue line = line.strip("\n") # check for header if (line[0] == '>'): if (current_domain): # previous domains/sites have already been found and parsed if (current_domain not in target_set): target_set[current_domain] = {} target_set[current_domain]['domains'] = [ d for d in seq_domains ] # previous domains target_set[current_domain]['sites'] = [ s for s in seq_sites ] # previous sites else: sys.stderr.write("*** phase error: %s duplicate\n"%(current_domain)) seq_domains = []; # current domains seq_sites = []; # current sites current_domain = line.split(' ')[0][1:] else: # check for data a_fields = line.split('\t') a_fields[0]=int(a_fields[0]) if (a_fields[1] == '-'): a_fields[2]=int(a_fields[2]) annot_info = dict(zip(('d_pos','type','d_end','descr'), a_fields)) re_df=re.compile(r' :(\d+)$') annot_info['descr'] = re_df.sub(r'~\1',annot_info['descr']) seq_domains.append(annot_info) else: annot_info = dict(zip(('d_pos','type', 'd_val', 'descr', a_fields))) annot_info['d_end'] = annot_info['d_pos'] seq_sites.append(annot_info) Reader.close() # get the last one if (current_domain): # previous domains/sites have already been found and parsed if (current_domain not in target_set): target_set[current_domain] = {} target_set[current_domain]['domains'] = [ d for d in seq_domains ] # previous domains target_set[current_domain]['sites'] = [ s for s in seq_sites ] # previous sites # else: # sys.stderr.write("*** phase error: %s duplicate\n"%(current_domain)) return target_set ################ # merge_annots(hit_r): # # take different annotations in hit_r and put them in one list # def merge_annots(hit_r): merged_annots = [] if ('q_aligned_domains' in hit_r): for annot in hit_r['q_aligned_domains']: annot['target']=1 merged_annots.append(annot) if ('aligned_domains' in hit_r): for annot in hit_r['aligned_domains']: annot['target']=0 merged_annots.append(annot) merged_annots = sorted(merged_annots, key=lambda x: x['qa_start']) return(merged_annots) ################ # get_file_annots(file_name) # def get_file_annots(file_name, hit_list): with open(file_name,'r') as Reader: ann_set = read_annots(Reader) return ann_set ################ # get_script_annots(script_name, hit_list) # # set up stdin/stdout pipe to send in hit list info and read results # def get_script_annots(script_name, hit_list, key_list): seq_set = {} proc = subprocess.Popen(script_name, stdin=subprocess.PIPE, stdout=subprocess.PIPE, shell=True, encoding='utf-8') for hit in hit_list: (seq_id, seq_len) = (hit[key_list[0]],hit[key_list[1]]) if (seq_id not in seq_set): proc.stdin.write(("%s\t%s\n"%(seq_id,seq_len))) proc.stdin.close() while (proc.returncode is None): proc.poll() return read_annots(proc.stdout) ################ # # link_annots(hit_list, annot_set) # # put 'domains' and 'sites' into each hit in the hit list # def link_annots(hit_list, annot_set): for hit in hit_list: seqid = hit['s_seq_id'] if (seqid in annot_set): if ('domains' in annot_set[seqid]): hit['domains']=annot_set[seqid]['domains'] if ('sites' in annot_set[seqid]): hit['sites']=annot_set[seqid]['sites'] # input: a blast BTOP string of the form: "1VA160TS7KG10RK27" # returns a list_ref of tokens: (1, "VA", 60, "TS", 7, "KG, 10, "RK", 27) # def decode_btop(btop_str): tokens = re.split(r'(\d+)',btop_str) # split with capture returns both strings between and separator (\d+) if not tokens[0]: tokens = tokens[1:] out_tokens = [] for token in tokens: if re.match(r'\d+',token): out_tokens.append(token) else: mis_tokens = re.split(r'(..)',token) # split with capture for mis in mis_tokens: if (mis): out_tokens.append(mis) return out_tokens def parse_query_lib(query_file): query_seqs = {} with open(query_file,"r") as qfd: header='' seq_data='' for line in qfd: line = line.strip("\n") if (line[0]=='>'): if (header): # save existing sequence seq_data = '' + seq_data query_seqs[header]=seq_data header = line[1:].split(' ')[0] else: line = re.sub(r'[^A-Za-z]','',line) seq_data += line.upper() # save last entry if (header): query_seqs[header]=seq_data return query_seqs # given: (1) a query sequence; (2) an encoded alignment; (3) a scoring matrix # calculate: # (1) the overall score # (2) a per residue dictionary of scores and mappings from query -> subject # (2) a per residue dictionary of scores and mappings from subject -> query # def alignment_score(query_r, hit, matrix_2d, matrix_diag, g_open, g_ext): query_start, subj_start = (int(hit['q_start']),int(hit['s_start'])) btop_align_r = decode_btop(hit['BTOP']) hit['btop_align'] = btop_align_r q_map = [] s_map = [] gap0, gap1 = (0 ,0) q_ix = query_start - 1 # start from zero s_ix = subj_start - 1 score, m_score = (0, 0) seq0, seq1 = ("","") for btop in btop_align_r: if (re.search(r'^\d+$',btop)): # matching query sequence, add it up for i in range(0,int(btop)): res = query_r[q_ix] score += matrix_diag[res] q_map.append({'s':score, 'y_ix':s_ix, 'res':res}) s_map.append({'s':score, 'y_ix':q_ix, 'res':res}) q_ix += 1 s_ix += 1 else: seq0, seq1 = (btop[0],btop[1]) if (re.search(r'\-',btop)): # is there a gap? if (seq0 == '-'): # is it in query? if gap0: # are we in a gap? score += g_ext else: score += g_open+g_ext gap0 = True # 'y_ix':-1 indicates alignment to gap s_map.append({'s':score, 'y_ix':-1, 'res':seq1}) s_ix += 1 else: # gap is in subject if gap1: score += g_ext else: score += g_open+g_ext gap1 = True # 'y_ix':-1 indicates alignment to gap q_map.append({'s':score, 'y_ix':-1, 'res':seq0}) q_ix += 1 else: # mismatch, not gap score += matrix_2d[seq0][seq1] gap1=gap0 = False q_map.append({'s':score, 'y_ix':s_ix, 'res':seq0}) s_map.append({'s':score, 'y_ix':q_ix, 'res':seq1}) q_ix += 1 s_ix += 1 return score, q_map, s_map ################################################################ # sub_alignment_stats -- calculate stats for ONE domain entry # given x_map, xa_start, xa_end where x=q/s depending on target # domain_r : domain boundaries # # calculate a score, identity and boundaries in both sequences and return values # def one_sub_alignment_stats(domain_r, x_map, y_map, xa_start, xa_end, ya_start, ya_end): td_start, td_end = (domain_r['d_pos'],domain_r['d_end']) if (td_end < xa_start or td_start > xa_end): return 0 if (td_start < xa_start): td_start = xa_start if (td_end > xa_end): td_end = xa_end td_start -= xa_start td_end -= xa_start left_score = 0 if (td_start>0) : left_score = x_map[td_start-1]['s'] score = x_map[td_end]['s'] - left_score # map[] coordinates are 0-based # ya_start = x_map[td_start]['y_ix']+1 # ya_end = x_map[td_end]['y_ix']+1 #### identity calculation: n_len = 0 n_id = 0 for xi in range(td_start, td_end+1): this_x = x_map[xi] x_res=this_x['res'] if (this_x['y_ix'] >= 0): n_len += 1 y_res = y_map[this_x['y_ix']-ya_start+1]['res'] if (x_res.upper() == y_res.upper()): n_id += 1 ident = float(n_id)/float(n_len) return score, ident, td_start+xa_start-1, td_end+xa_start-1, x_map[td_start]['y_ix'], x_map[td_end]['y_ix'] ################ # get domain scores, idents, boundaries for list of domains # def do_sub_alignment_stats(domain_list, x_map, y_map, xa_start, xa_end, ya_start, ya_end, keys_str): aligned_doms = [] for domain in domain_list: subalign_data = one_sub_alignment_stats(domain, x_map, y_map, xa_start, xa_end, ya_start, ya_end) if (subalign_data and len(subalign_data)==6): sub_data = dict(zip(keys_str,subalign_data)) for k in ('type','descr'): sub_data[k] = domain[k] aligned_doms.append(sub_data) return(aligned_doms) #### # print raw domain info: # |DX:%d-%d;C=dom_info|XD:%d-%d:C=dom_info # def format_dom_info(q_dom_r, dom_r): dom_str = "" for dom in q_dom_r: dom_str += "|DX:%d-%d;C=%s"%(dom['d_pos'],dom['d_end'], dom['descr']) for dom in dom_r: dom_str += "|XD:%d-%d;C=%s"%(dom['d_pos'],dom['d_end'], dom['descr']) return dom_str def format_annot_info(annot_list_r, hit): annot_str = ""; # two types of annotations, domains and sites. score_scale = hit['score']/hit['raw_score'] for annot_r in (annot_list_r ): if (annot_r['type'] == '-'): fsub_score = annot_r['score']/hit['raw_score'] ns_score, s_bit = (int(annot_r['score'] * score_scale + 0.5), fsub_score * hit['bits']) qval = 0.0 if (hit['evalue'] == 0.0): if (s_bit > 50.0): qval = 3000.0 else: qval = -10.0 * (2.0*log(400.0) + s_bit)/log(10.0) else: qval = -10.0*log(hit['evalue'])*fsub_score/log(10.0) if qval < 0.0: qval = 0.0 rx_str = 'XR' if (annot_r['target']): rx_str = "RX" annot_str += ';'.join(("|%s:%d-%d:%d-%d:s=%d"%(rx_str, annot_r['qa_start']+1,annot_r['qa_end']+1, annot_r['sa_start']+1,annot_r['sa_end']+1,ns_score), "b=%.1f"%(s_bit),"I=%.3f"%(annot_r['ident']), "Q=%.1f"%(qval),"C=%s"%(annot_r['descr']))) else: # site annotation ann_type = annot_r['type']; site_str = "|%cX"%(ann_type) if (annot_r['target'] == 1): site_str = "|X%c"%(ann_type) elif (annot_r['target'] == 2): site_str = "|%c%c"%(ann_type, ann_type) annot_str += "%s:"%(site_str) annot_str += "%d%s%s%d%s"%(annot_r['qa_pos'], annot_r['q_res'], annot_r['m_symb'], annot_r['sa_pos'], annot_r['s_res']) return annot_str def main(args): blosum62, blosum62_diag, g_open, g_ext = init_blosum62() if (args.query_file): # query_lib_r has a set of query sequences query_lib_r = parse_query_lib(args.query_file) else: sys.stderr.write("--query required\n") exit(1) tab_fields = "q_seqid s_seqid percid alen mismatch gopen q_start q_end s_start s_end evalue bits BTOP".split(' ') int_fields = "alen mismatch gopen q_start q_end s_start s_end".split(' ') float_fields = "percid evalue bits score".split(' ') if (args.have_qslen): tab_fields = "q_seqid q_len s_seqid s_len percid alen mismatch gopen q_start q_end s_start s_end evalue bits BTOP".split(' ') int_fields = "q_len s_len alen mismatch gopen q_start q_end s_start s_end".split(' ') # the fields that are displayed are listed here. By default, all fields except score and BTOP are displayed. out_tab_fields = tab_fields[0:-1] in_tab_fields = tab_fields[0:-1] if (args.raw_out): out_tab_fields.append("raw_score") if (args.raw_in): in_tab_fields.append("score") ## always add BTOP in_tab_fields.append("BTOP") tab_fields = in_tab_fields if (args.out_fields): out_tab_fields = out_fields.split(" ") header_lines = [] next_line = "" have_data = False hit_list = [] q_hit_list = [] for line in fileinput.input(args.files): if (line[0] == '#'): if (have_data): next_line = line have_data = False break else: header_lines.append(line) continue have_data = True line = line.strip('\n') if (line): this_data = dict(zip(tab_fields, line.split("\t"))) for k in this_data.keys(): if (k in int_fields): this_data[k] = int(this_data[k]) if (k in float_fields): this_data[k] = float(this_data[k]) hit_list.append(this_data) # get the query annotations q_hit_list = [] if (args.q_ann_file): q_seqid = hit_list[0]['q_seqid'] q_hit_list.append({'s_seq_id':q_seqid, 's_end':len(query_lib_r[q_seqid])}) q_annots = get_file_annots(args.q_ann_file, q_hit_list) link_annots(q_hit_list, q_annots) elif (args.q_ann_script): args.q_ann_script = re.sub(r'\+',' ',args.q_ann_script) if (args.q_ann_script and shutil.which(args.q_ann_script.split(" ")[0])): q_seqid = hit_list[0]['q_seqid'] q_hit_list.append({'s_seq_id':q_seqid, 's_end':len(query_lib_r[q_seqid])}) q_annots = get_script_annots(args.q_ann_script, q_hit_list, ['s_seq_id','s_end']) link_annots(q_hit_list, q_annots) # get the subject annotations # first set up the list with sequence lengths if (args.ann_file or args.ann_script): s_len = 100000 for hit in hit_list: hit['s_seq_id']=hit['s_seqid'] if (not args.have_qslen): hit['s_end']=s_len if (args.ann_file): s_annots = get_file_annots(args.ann_file, hit_list) link_annots(hit_list, s_annots) elif (args.ann_script): args.ann_script = re.sub(r'\+',' ',args.ann_script) if (shutil.which(args.ann_script.split(" ")[0])): s_annots = get_script_annots(args.ann_script, hit_list,['s_seq_id','s_end']) link_annots(hit_list, s_annots) for line in header_lines: print(line, end='') header_lines = [next_line] # now get query annotation if available for hit in hit_list: list_covered = [] # If I have an encoded aligment {BTOP} and a query sequence query_lib_r && query_lib_r[hit['q_seqid']] # then I can calculate sub-alignment scores if ('BTOP' in hit and query_lib_r and hit['q_seqid'] in query_lib_r): # calculate raw_score and mappings hit['raw_score'], q_map, s_map = alignment_score(query_lib_r[hit['q_seqid']], hit,blosum62, blosum62_diag, g_open, g_ext) if ('score' not in hit): hit['score'] = hit['raw_score'] # calculate sub-alignment scores in subject/library coordinates if ('domains' in hit and len(hit['domains'])>0): hit['aligned_domains'] = do_sub_alignment_stats(hit['domains'], s_map, q_map, hit['s_start'],hit['s_end'],hit['q_start'],hit['q_end'], ('score','ident','sa_start', 'sa_end', 'qa_start', 'qa_end')) # calculate sub-alignment scores in query coordinates if (len(q_hit_list) > 0 and 'domains' in q_hit_list[0] and len(q_hit_list[0]['domains'])>0): hit['q_aligned_domains'] = do_sub_alignment_stats(q_hit_list[0]['domains'], q_map, s_map, hit['q_start'],hit['q_end'],hit['s_start'],hit['s_end'], ('score','ident','qa_start', 'qa_end', 'sa_start', 'sa_end')) ################ ## final output display print("\t".join([str(hit[x]) for x in out_tab_fields]),end='') # show fields from original blast tabular file merged_annots_r = merge_annots(hit) # merge the four possible annotation lists into one. if (len(merged_annots_r)>0): print("\t"+format_annot_info(merged_annots_r, hit),end='') if (args.dom_info): if (len(q_hit_list) > 0 and 'domains' in q_hit_list[0]): print("\t"+format_dom_info(q_hit_list[0]['domains'], hit['domains']),end='') else: print("\t"+format_dom_info([], hit['domains']),end='') elif (len(list_covered)>0): print("\t" + ";".join(list_covered)) if (args.dom_info): print("\t"+format_dom_info(q_hit_list[0]['domains'], hit['domains']),end='') print() for line in header_lines: print(line,end="") if __name__ == '__main__': print('# ' + ' '.join(sys.argv)) parser=argparse.ArgumentParser(description='annot_blast_btop4.py : annotate blast tabular format with BTOP ') # not implemented # parser.add_argument('--matrix', help='scoring matrix',dest='matrix',action='store',default='BL62') parser.add_argument('--ann_script', help='script for subject annotations',dest='ann_script',action='store') parser.add_argument('--q_ann_script', help='script for query annotations',dest='q_ann_script',action='store') parser.add_argument('--ann_file', help='subject annotation file',dest='ann_file',action='store') parser.add_argument('--q_ann_file', help='query annotation file',dest='q_ann_file',action='store') parser.add_argument('--have_qslen', help='query/subject lenghts in tab file',dest='have_qslen',action='store_true',default=False) parser.add_argument('--dom_info', help='show unaligned domain coordinates',dest='dom_info',action='store_true',default=False) parser.add_argument('--sub2query', help='get query annots from self-subject',dest='sub_query',action='store_true',default=False) parser.add_argument('--query', help='file of query sequences',dest='query_file',action='store') parser.add_argument('--out_fields', help='names/order of output fields',dest='out_fields',action='store') parser.add_argument('--raw_score', help='raw score after bit score',dest='raw_in',action='store_true',default=True) parser.add_argument('--no_raw_score', help='raw score after bit score',dest='raw_in',action='store_false', default=True) parser.add_argument('--no-raw_score', help='raw score after bit score',dest='raw_in',action='store_false', default=True) parser.add_argument('--raw_score_out', help='display raw score',dest='raw_out',action='store_true',default=False) parser.add_argument('files', metavar='FILE', help='Blast tabular BTOP files to read', nargs='*') args=parser.parse_args() main(args)