Mercurial > repos > tyty > structurefold

comparison predict/predict_RNAs.py @ 64:a1ce42d5258d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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