|
59
|
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
|
