|
4
|
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
|
|
|
10
|
|
|
11 id_file = sys.argv[1]
|
|
|
12 seq_file = sys.argv[2]
|
|
|
13 output_file = sys.argv[4]
|
|
|
14
|
|
|
15
|
|
|
16 flag = 0
|
|
|
17 if sys.argv[3]!='None': #input reactivity file if provided
|
|
|
18 react_file = sys.argv[3]
|
|
|
19 react = parse_dist(react_file)
|
|
|
20 react = react[1]
|
|
|
21 flag = 1
|
|
|
22
|
|
|
23 ospath = os.path.realpath(sys.argv[0])
|
|
|
24 ost = ospath.split('/')
|
|
|
25 syspath = ""
|
|
|
26 for i in range(len(ost)-1):
|
|
|
27 syspath = syspath+ost[i].strip()
|
|
|
28 syspath = syspath+'/'
|
|
|
29
|
|
|
30 ids = read_t_file(id_file)
|
|
|
31 sequences = SeqIO.parse(seq_file, 'fasta')
|
|
|
32
|
|
|
33
|
|
|
34 seqs = {}
|
|
|
35 for seq in sequences:
|
|
|
36 seqs[seq.id] = seq.seq.tostring()
|
|
|
37
|
|
|
38 if len(ids)>10: #setup a limit of the number of sequence to be predicted
|
|
|
39 print("Number of sequences exceeds limitation!")
|
|
|
40 sys.exit(0)
|
|
|
41
|
|
|
42
|
|
|
43 #predict RNA structures
|
|
|
44 os.system("mkdir "+syspath+"output_f")
|
|
|
45 for i in range(len(ids)):
|
|
|
46 id_s = ids[i][0]
|
|
|
47 print(id_s)
|
|
|
48 #Put RNA sequence and reactivities into files
|
|
|
49 if id_s in seqs:
|
|
|
50 f = file(syspath+"temp.txt", 'w')
|
|
|
51 f.write('>'+id_s)
|
|
|
52 f.write('\n')
|
|
|
53 f.write(seqs[id_s])
|
|
|
54 f.close()
|
|
|
55 if flag == 0:
|
|
|
56 os.system("Fold "+syspath+"temp.txt"+" "+syspath+"output_f/"+id_s+".ct")
|
|
|
57 if flag == 1:
|
|
|
58 if id_s in react:
|
|
|
59 f = file(syspath+"constraint.txt",'w')
|
|
|
60 make_plot(react[id_s],id_s,(syspath+"output_f/")) #make a plot of the distribution of the reactivites of the input RNA
|
|
|
61 #h = file(syspath+"output_f/transcript_reactivities.txt", 'w')
|
|
|
62 #h.write(id_s)
|
|
|
63 #h.write('\n')
|
|
|
64 for j in range(0, (len(react[id_s]))):
|
|
|
65 if react[id_s][j]!='NA':
|
|
|
66 f.write(str(j+1))
|
|
|
67 f.write('\t')
|
|
|
68 f.write(str(react[id_s][j]))
|
|
|
69 f.write('\n')
|
|
|
70 #h.write(str(react[id_s][j])) #Output the reactivities
|
|
|
71 #h.write('\t')
|
|
|
72 f.close()
|
|
|
73 #h.write('\n')
|
|
|
74 #h.write('\n')
|
|
|
75 os.system("Fold "+syspath+"temp.txt"+" -sh"+" "+syspath+"constraint.txt"+" "+syspath+"output_f/"+id_s+".ct")
|
|
|
76 else:
|
|
|
77 print(id_s+" not in the data of react!")
|
|
|
78 os.system("draw "+syspath+"output_f/"+id_s+".ct "+syspath+"output_f/"+id_s+".ps")
|
|
|
79 else:
|
|
|
80 print(id_s+" not in the data of sequences!")
|
|
|
81
|
|
|
82 #Remove the unnecessary files
|
|
|
83 os.system("tar -zcvPf "+output_file+" "+syspath+"output_f/"+"*.*")
|
|
|
84 os.system("rm -f "+syspath+"temp.txt")
|
|
|
85 os.system("rm -r "+syspath+"output_f")
|
|
|
86 if flag == 1:
|
|
|
87 os.system("rm -f "+syspath+"constraint.txt")
|
|
|
88 # h.close()
|
|
|
89
|
|
|
90
|
|
|
91
|
|
|
92
|
|
|
93
|
