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