|
0
|
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) |
