annotate predict/predict_RNAs.py @ 39:c5f80e82efb7 draft

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