Mercurial > repos > estrain > seqsero_v1

comparison SeqSero/libs/special_gene_test_assemblies.py @ 0:c577b57b7c74 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author estrain
date Wed, 06 Dec 2017 15:59:29 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:c577b57b7c74
1 #just an possible use, we can use it to replace H**.py, treat fliC and fljB as the target genes?
2
3 from __future__ import division
4 import sys
5 import os
6 from Bio.Blast import NCBIXML
7 from Bio import SeqIO
8
9 def special_gene(target_fie,database,gene_list):
10 database=database.split("/")[-1]##########1/27/2015
11 os.system('makeblastdb -in database/'+database+' -out '+database+'_db -dbtype nucl')##########1/28/2015
12 os.system('blastn -query '+target_file+' -db '+database+'_db -out '+database+'_vs_'+target_file+'.xml '+'-outfmt 5')##########1/28/2015
13 xml_file=database+'_vs_'+target_file+'.xml'
14 result_handle=open(xml_file)
15 blast_record=NCBIXML.parse(result_handle)
16 records=list(blast_record)
17 E_thresh=1e-10
18 for x in gene_list:
19 handle=SeqIO.parse("database/"+database,"fasta")##########1/28/2015
20 length_list=[]
21 for y in handle:
22 if x in y.description:
23 length_x=len(y.seq)
24 length_list.append(length_x)
25 aver_len=float(sum(length_list))/len(length_list)
26 hspbit=[]
27 alignmentlist=[]
28 for record in records:
29 for alignment in record.alignments:
30 if x in alignment.hit_def: #multi gene database, so...
31 print x,"got a hit, evaluating the hit quality..."
32 score=0
33 for hsp in alignment.hsps:
34 if hsp.expect<E_thresh:
35 score+=hsp.bits
36 alignment=alignment.hit_def+':'+str(score)
37 hspbit.append(score)
38 alignmentlist.append(alignment)
39 scorelist=dict(zip(alignmentlist,hspbit))
40 score=0
41 for Htype in scorelist:
42 if scorelist[Htype]>score:
43 First_Choice=Htype
44 score=scorelist[Htype]
45 if float(score)>=0.1*aver_len:
46 print "$$$",First_Choice,"got a hit, score:",score
47 else:
48 print "$$$No ",x,"exists"
49 os.system("rm "+database+"_db.*")##########1/28/2015
50 os.system("rm "+xml_file)##########1/28/2015
51
52 database=sys.argv[1]
53 target_file=sys.argv[2]
54 gene_list=[]
55 a=1
56 i=3
57 while a==1:
58 try:
59 gene_list.append(sys.argv[i])
60 i+=1
61 except:
62 a=0
63 special_gene(target_file,database,gene_list)