Mercurial > repos > tyty > structurefold

comparison predict/predict_RNAs.py @ 4:a292aaf51735 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author tyty
date Mon, 15 Sep 2014 14:52:43 -0400
parents
children
comparison
equal deleted inserted replaced
3:f4cc06e92530 4:a292aaf51735
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