annotate predict/predict_RNAs.py @ 59:afd114ef8857 draft

Uploaded
author tyty
date Tue, 18 Nov 2014 01:01:52 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
59
afd114ef8857 Uploaded
tyty
parents:
diff changeset
1 #RNA structure prediction & Output and illustrate reactivities
afd114ef8857 Uploaded
tyty
parents:
diff changeset
2
afd114ef8857 Uploaded
tyty
parents:
diff changeset
3 import sys
afd114ef8857 Uploaded
tyty
parents:
diff changeset
4 from parse_dis_pac import *
afd114ef8857 Uploaded
tyty
parents:
diff changeset
5 from read_file import *
afd114ef8857 Uploaded
tyty
parents:
diff changeset
6 from Bio import SeqIO
afd114ef8857 Uploaded
tyty
parents:
diff changeset
7 import os
afd114ef8857 Uploaded
tyty
parents:
diff changeset
8 from rtts_plot import *
afd114ef8857 Uploaded
tyty
parents:
diff changeset
9 import random
afd114ef8857 Uploaded
tyty
parents:
diff changeset
10 import string
afd114ef8857 Uploaded
tyty
parents:
diff changeset
11
afd114ef8857 Uploaded
tyty
parents:
diff changeset
12
afd114ef8857 Uploaded
tyty
parents:
diff changeset
13 id_file = sys.argv[1]
afd114ef8857 Uploaded
tyty
parents:
diff changeset
14 seq_file = sys.argv[2]
afd114ef8857 Uploaded
tyty
parents:
diff changeset
15 output_file = sys.argv[4]
afd114ef8857 Uploaded
tyty
parents:
diff changeset
16
afd114ef8857 Uploaded
tyty
parents:
diff changeset
17
afd114ef8857 Uploaded
tyty
parents:
diff changeset
18 flag = 0
afd114ef8857 Uploaded
tyty
parents:
diff changeset
19 if sys.argv[3]!='None': #input reactivity file if provided
afd114ef8857 Uploaded
tyty
parents:
diff changeset
20 react_file = sys.argv[3]
afd114ef8857 Uploaded
tyty
parents:
diff changeset
21 react = parse_dist(react_file)
afd114ef8857 Uploaded
tyty
parents:
diff changeset
22 react = react[1]
afd114ef8857 Uploaded
tyty
parents:
diff changeset
23 flag = 1
afd114ef8857 Uploaded
tyty
parents:
diff changeset
24
afd114ef8857 Uploaded
tyty
parents:
diff changeset
25 ospath = os.path.realpath(sys.argv[0])
afd114ef8857 Uploaded
tyty
parents:
diff changeset
26 ost = ospath.split('/')
afd114ef8857 Uploaded
tyty
parents:
diff changeset
27 syspath = ""
afd114ef8857 Uploaded
tyty
parents:
diff changeset
28 for i in range(len(ost)-1):
afd114ef8857 Uploaded
tyty
parents:
diff changeset
29 syspath = syspath+ost[i].strip()
afd114ef8857 Uploaded
tyty
parents:
diff changeset
30 syspath = syspath+'/'
afd114ef8857 Uploaded
tyty
parents:
diff changeset
31
afd114ef8857 Uploaded
tyty
parents:
diff changeset
32 ids = read_t_file(id_file)
afd114ef8857 Uploaded
tyty
parents:
diff changeset
33 sequences = SeqIO.parse(seq_file, 'fasta')
afd114ef8857 Uploaded
tyty
parents:
diff changeset
34
afd114ef8857 Uploaded
tyty
parents:
diff changeset
35 rs = ''.join(random.sample(string.ascii_letters + string.digits, 8))
afd114ef8857 Uploaded
tyty
parents:
diff changeset
36
afd114ef8857 Uploaded
tyty
parents:
diff changeset
37
afd114ef8857 Uploaded
tyty
parents:
diff changeset
38 seqs = {}
afd114ef8857 Uploaded
tyty
parents:
diff changeset
39 for seq in sequences:
afd114ef8857 Uploaded
tyty
parents:
diff changeset
40 seqs[seq.id] = seq.seq.tostring()
afd114ef8857 Uploaded
tyty
parents:
diff changeset
41
afd114ef8857 Uploaded
tyty
parents:
diff changeset
42 if len(ids)>100: #setup a limit of the number of sequence to be predicted
afd114ef8857 Uploaded
tyty
parents:
diff changeset
43 print("Number of sequences exceeds limitation!")
afd114ef8857 Uploaded
tyty
parents:
diff changeset
44 sys.exit(0)
afd114ef8857 Uploaded
tyty
parents:
diff changeset
45
afd114ef8857 Uploaded
tyty
parents:
diff changeset
46
afd114ef8857 Uploaded
tyty
parents:
diff changeset
47 #predict RNA structures
afd114ef8857 Uploaded
tyty
parents:
diff changeset
48 os.system("mkdir "+syspath+"output_"+rs)
afd114ef8857 Uploaded
tyty
parents:
diff changeset
49 for i in range(len(ids)):
afd114ef8857 Uploaded
tyty
parents:
diff changeset
50 id_s = ids[i][0]
afd114ef8857 Uploaded
tyty
parents:
diff changeset
51 print(id_s)
afd114ef8857 Uploaded
tyty
parents:
diff changeset
52 #Put RNA sequence and reactivities into files
afd114ef8857 Uploaded
tyty
parents:
diff changeset
53 if id_s in seqs:
afd114ef8857 Uploaded
tyty
parents:
diff changeset
54 f = file(syspath+"temp.txt", 'w')
afd114ef8857 Uploaded
tyty
parents:
diff changeset
55 f.write('>'+id_s)
afd114ef8857 Uploaded
tyty
parents:
diff changeset
56 f.write('\n')
afd114ef8857 Uploaded
tyty
parents:
diff changeset
57 f.write(seqs[id_s])
afd114ef8857 Uploaded
tyty
parents:
diff changeset
58 f.close()
afd114ef8857 Uploaded
tyty
parents:
diff changeset
59 if flag == 0:
afd114ef8857 Uploaded
tyty
parents:
diff changeset
60 os.system("Fold "+syspath+"temp.txt"+" "+syspath+"output_"+rs+"/"+id_s+".ct")
afd114ef8857 Uploaded
tyty
parents:
diff changeset
61 if flag == 1:
afd114ef8857 Uploaded
tyty
parents:
diff changeset
62 if id_s in react:
afd114ef8857 Uploaded
tyty
parents:
diff changeset
63 f = file(syspath+"constraint.txt",'w')
afd114ef8857 Uploaded
tyty
parents:
diff changeset
64 make_plot(react[id_s],id_s,(syspath+"output_"+rs+"/")) #make a plot of the distribution of the reactivites of the input RNA
afd114ef8857 Uploaded
tyty
parents:
diff changeset
65 #h = file(syspath+"output_f/transcript_reactivities.txt", 'w')
afd114ef8857 Uploaded
tyty
parents:
diff changeset
66 #h.write(id_s)
afd114ef8857 Uploaded
tyty
parents:
diff changeset
67 #h.write('\n')
afd114ef8857 Uploaded
tyty
parents:
diff changeset
68 for j in range(0, (len(react[id_s]))):
afd114ef8857 Uploaded
tyty
parents:
diff changeset
69 if react[id_s][j]!='NA':
afd114ef8857 Uploaded
tyty
parents:
diff changeset
70 f.write(str(j+1))
afd114ef8857 Uploaded
tyty
parents:
diff changeset
71 f.write('\t')
afd114ef8857 Uploaded
tyty
parents:
diff changeset
72 f.write(str(react[id_s][j]))
afd114ef8857 Uploaded
tyty
parents:
diff changeset
73 f.write('\n')
afd114ef8857 Uploaded
tyty
parents:
diff changeset
74 #h.write(str(react[id_s][j])) #Output the reactivities
afd114ef8857 Uploaded
tyty
parents:
diff changeset
75 #h.write('\t')
afd114ef8857 Uploaded
tyty
parents:
diff changeset
76 f.close()
afd114ef8857 Uploaded
tyty
parents:
diff changeset
77 #h.write('\n')
afd114ef8857 Uploaded
tyty
parents:
diff changeset
78 #h.write('\n')
afd114ef8857 Uploaded
tyty
parents:
diff changeset
79 os.system("Fold "+syspath+"temp.txt"+" -sh"+" "+syspath+"constraint.txt"+" "+syspath+"output_"+rs+"/"+id_s+".ct")
afd114ef8857 Uploaded
tyty
parents:
diff changeset
80 else:
afd114ef8857 Uploaded
tyty
parents:
diff changeset
81 print(id_s+" not in the data of react!")
afd114ef8857 Uploaded
tyty
parents:
diff changeset
82 os.system("draw "+syspath+"output_"+rs+"/"+id_s+".ct "+syspath+"output_"+rs+"/"+id_s+".ps")
afd114ef8857 Uploaded
tyty
parents:
diff changeset
83 else:
afd114ef8857 Uploaded
tyty
parents:
diff changeset
84 print(id_s+" not in the data of sequences!")
afd114ef8857 Uploaded
tyty
parents:
diff changeset
85
afd114ef8857 Uploaded
tyty
parents:
diff changeset
86 #Remove the unnecessary files
afd114ef8857 Uploaded
tyty
parents:
diff changeset
87 os.system("tar -zcvPf "+output_file+" "+syspath+"output_"+rs+"/"+"*.* 2>"+syspath+"log.txt")
afd114ef8857 Uploaded
tyty
parents:
diff changeset
88 os.system("rm -f "+syspath+"temp.txt")
afd114ef8857 Uploaded
tyty
parents:
diff changeset
89 os.system("rm -r "+syspath+"output_"+rs)
afd114ef8857 Uploaded
tyty
parents:
diff changeset
90 if flag == 1:
afd114ef8857 Uploaded
tyty
parents:
diff changeset
91 os.system("rm -f "+syspath+"constraint.txt")
afd114ef8857 Uploaded
tyty
parents:
diff changeset
92 # h.close()
afd114ef8857 Uploaded
tyty
parents:
diff changeset
93
afd114ef8857 Uploaded
tyty
parents:
diff changeset
94
afd114ef8857 Uploaded
tyty
parents:
diff changeset
95
afd114ef8857 Uploaded
tyty
parents:
diff changeset
96
afd114ef8857 Uploaded
tyty
parents:
diff changeset
97