Mercurial > repos > tyty > structurefold
comparison predict/predict_RNAs.py @ 64:a1ce42d5258d draft
Uploaded
author | tyty |
---|---|
date | Tue, 18 Nov 2014 15:54:31 -0500 |
parents | |
children | 96a827962750 |
comparison
equal
deleted
inserted
replaced
63:c1f1b552c1b8 | 64:a1ce42d5258d |
---|---|
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 |