Mercurial > repos > tyty > structurefold
comparison predict/predict_RNAs.py @ 59:afd114ef8857 draft
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 |