Mercurial > repos > estrain > seqsero2_v1_0_1

comparison SeqSero2_package.py @ 13:66277925a124

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author estrain
date Sat, 07 Sep 2019 11:12:04 -0400
parents a3accb97716a
children eca5802b1d6e
comparison
equal deleted inserted replaced
12:b43b7ef0df6a 13:66277925a124
10 import pickle 10 import pickle
11 import argparse 11 import argparse
12 import itertools 12 import itertools
13 from distutils.version import LooseVersion 13 from distutils.version import LooseVersion
14 from distutils.spawn import find_executable 14 from distutils.spawn import find_executable
15 import math
15 16
16 try: 17 try:
17 from .version import SeqSero2_version 18 from .version import SeqSero2_version
18 except Exception: #ImportError 19 except Exception: #ImportError
19 from version import SeqSero2_version 20 from version import SeqSero2_version
516 cover_region=set() #fixed problem that repeated calculation leading percentage > 1 517 cover_region=set() #fixed problem that repeated calculation leading percentage > 1
517 List.append(handle[i].query.strip()+"___"+handle[i].alignments[j].hit_def) 518 List.append(handle[i].query.strip()+"___"+handle[i].alignments[j].hit_def)
518 for z in range(len(handle[i].alignments[j].hsps)): 519 for z in range(len(handle[i].alignments[j].hsps)):
519 hsp=handle[i].alignments[j].hsps[z] 520 hsp=handle[i].alignments[j].hsps[z]
520 temp=set(range(hsp.query_start,hsp.query_end)) 521 temp=set(range(hsp.query_start,hsp.query_end))
522 # Bit score is not parsed correctly, calculate it from raw score
523 est_bit=(handle[i].ka_params[0]*hsp.score-math.log(handle[i].ka_params[1]))/(math.log(2))
521 if len(cover_region)==0: 524 if len(cover_region)==0:
522 cover_region=cover_region|temp 525 cover_region=cover_region|temp
523 fraction=1 526 fraction=1
524 else: 527 else:
525 fraction=1-len(cover_region&temp)/float(len(temp)) 528 fraction=1-len(cover_region&temp)/float(len(temp))
526 cover_region=cover_region|temp 529 cover_region=cover_region|temp
527 if "last" in handle[i].query or "first" in handle[i].query: 530 if "last" in handle[i].query or "first" in handle[i].query:
528 score+=hsp.bits*fraction 531 #score+=hsp.bits*fraction
532 score+=est_bit*fraction
529 ids+=float(hsp.identities)/handle[i].query_length*fraction 533 ids+=float(hsp.identities)/handle[i].query_length*fraction
530 else: 534 else:
531 score+=hsp.bits*fraction 535 #score+=hsp.bits*fraction
536 score+=est_bit*fraction
532 ids+=float(hsp.identities)/handle[i].query_length*fraction 537 ids+=float(hsp.identities)/handle[i].query_length*fraction
533 List_score.append(score) 538 List_score.append(score)
534 List_ids.append(ids) 539 List_ids.append(ids)
535 List_query_region.append(cover_region) 540 List_query_region.append(cover_region)
536 temp=zip(List,List_score,List_ids,List_query_region) 541 temp=zip(List,List_score,List_ids,List_query_region)