Mercurial > repos > tyty > structurefold

comparison predict/predict_RNAs.py @ 59:afd114ef8857 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author tyty
date Tue, 18 Nov 2014 01:01:52 -0500
parents
children
comparison
equal deleted inserted replaced
58:971ef91f6b68 59:afd114ef8857
1 #RNA structure prediction & Output and illustrate reactivities
2
3 import sys
4 from parse_dis_pac import *
5 from read_file import *
6 from Bio import SeqIO
7 import os
8 from rtts_plot import *
9 import random
10 import string
11
12
13 id_file = sys.argv[1]
14 seq_file = sys.argv[2]
15 output_file = sys.argv[4]
16
17
18 flag = 0
19 if sys.argv[3]!='None': #input reactivity file if provided
20 react_file = sys.argv[3]
21 react = parse_dist(react_file)
22 react = react[1]
23 flag = 1
24
25 ospath = os.path.realpath(sys.argv[0])
26 ost = ospath.split('/')
27 syspath = ""
28 for i in range(len(ost)-1):
29 syspath = syspath+ost[i].strip()
30 syspath = syspath+'/'
31
32 ids = read_t_file(id_file)
33 sequences = SeqIO.parse(seq_file, 'fasta')
34
35 rs = ''.join(random.sample(string.ascii_letters + string.digits, 8))
36
37
38 seqs = {}
39 for seq in sequences:
40 seqs[seq.id] = seq.seq.tostring()
41
42 if len(ids)>100: #setup a limit of the number of sequence to be predicted
43 print("Number of sequences exceeds limitation!")
44 sys.exit(0)
45
46
47 #predict RNA structures
48 os.system("mkdir "+syspath+"output_"+rs)
49 for i in range(len(ids)):
50 id_s = ids[i][0]
51 print(id_s)
52 #Put RNA sequence and reactivities into files
53 if id_s in seqs:
54 f = file(syspath+"temp.txt", 'w')
55 f.write('>'+id_s)
56 f.write('\n')
57 f.write(seqs[id_s])
58 f.close()
59 if flag == 0:
60 os.system("Fold "+syspath+"temp.txt"+" "+syspath+"output_"+rs+"/"+id_s+".ct")
61 if flag == 1:
62 if id_s in react:
63 f = file(syspath+"constraint.txt",'w')
64 make_plot(react[id_s],id_s,(syspath+"output_"+rs+"/")) #make a plot of the distribution of the reactivites of the input RNA
65 #h = file(syspath+"output_f/transcript_reactivities.txt", 'w')
66 #h.write(id_s)
67 #h.write('\n')
68 for j in range(0, (len(react[id_s]))):
69 if react[id_s][j]!='NA':
70 f.write(str(j+1))
71 f.write('\t')
72 f.write(str(react[id_s][j]))
73 f.write('\n')
74 #h.write(str(react[id_s][j])) #Output the reactivities
75 #h.write('\t')
76 f.close()
77 #h.write('\n')
78 #h.write('\n')
79 os.system("Fold "+syspath+"temp.txt"+" -sh"+" "+syspath+"constraint.txt"+" "+syspath+"output_"+rs+"/"+id_s+".ct")
80 else:
81 print(id_s+" not in the data of react!")
82 os.system("draw "+syspath+"output_"+rs+"/"+id_s+".ct "+syspath+"output_"+rs+"/"+id_s+".ps")
83 else:
84 print(id_s+" not in the data of sequences!")
85
86 #Remove the unnecessary files
87 os.system("tar -zcvPf "+output_file+" "+syspath+"output_"+rs+"/"+"*.* 2>"+syspath+"log.txt")
88 os.system("rm -f "+syspath+"temp.txt")
89 os.system("rm -r "+syspath+"output_"+rs)
90 if flag == 1:
91 os.system("rm -f "+syspath+"constraint.txt")
92 # h.close()
93
94
95
96
97