Mercurial > repos > tyty > structurefold

comparison predict/predict_RNAs.py @ 67:96a827962750 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded updated scripts, removed *.pyc and .DS_Store
author devteam
date Fri, 21 Nov 2014 11:28:59 -0500
parents a1ce42d5258d
children fb80870002a3
comparison
equal deleted inserted replaced
66:d2817a631a7b 67:96a827962750
1 #RNA structure prediction & Output and illustrate reactivities 1 #RNA structure prediction & Output and illustrate reactivities
2 2
3 import sys 3 import sys
4 import shlex
5 import subprocess
6 import tarfile
4 from parse_dis_pac import * 7 from parse_dis_pac import *
5 from read_file import * 8 from read_file import *
6 from Bio import SeqIO 9 from Bio import SeqIO
7 import os 10 import os
8 from rtts_plot import * 11 from rtts_plot import *
13 id_file = sys.argv[1] 16 id_file = sys.argv[1]
14 seq_file = sys.argv[2] 17 seq_file = sys.argv[2]
15 output_file = sys.argv[4] 18 output_file = sys.argv[4]
16 19
17 20
18 flag = 0 21 flag = False
19 if sys.argv[3]!='None': #input reactivity file if provided 22 if sys.argv[3]!='None': #input reactivity file if provided
20 react_file = sys.argv[3] 23 react_file = sys.argv[3]
21 react = parse_dist(react_file) 24 react = parse_dist(react_file)
22 react = react[1] 25 react = react[1]
23 flag = 1 26 flag = True
24 27
25 syspath = os.getcwd() 28 syspath = os.getcwd()
26 29
27 ids = read_t_file(id_file) 30 ids = read_t_file(id_file)
28 sequences = SeqIO.parse(seq_file, 'fasta') 31 sequences = SeqIO.parse(seq_file, 'fasta')
36 print("Number of sequences exceeds limitation!") 39 print("Number of sequences exceeds limitation!")
37 sys.exit(0) 40 sys.exit(0)
38 41
39 42
40 #predict RNA structures 43 #predict RNA structures
41 output_directory = os.path.join(syspath, "output_files/") 44 output_directory = os.path.join(syspath, "output_files")
42 os.makedirs(output_directory) 45 if not os.path.exists(output_directory):
46 os.makedirs(output_directory)
43 for i in range(len(ids)): 47 for i in range(len(ids)):
44 id_s = ids[i][0] 48 id_s = ids[i][0]
45 print(id_s) 49 print(id_s)
46 #Put RNA sequence and reactivities into files 50 #Put RNA sequence and reactivities into files
47 if id_s in seqs: 51 if id_s in seqs:
48 f = file(syspath+"temp.txt", 'w') 52 fh = file(os.path.join(syspath,"temp.txt"), 'w')
49 f.write('>'+id_s) 53 fh.write('>'+id_s)
50 f.write('\n') 54 fh.write('\n')
51 f.write(seqs[id_s]) 55 fh.write(seqs[id_s])
52 f.close() 56 fh.close()
53 if flag == 0: 57 if not flag:
54 os.system("Fold "+syspath+"temp.txt"+" "+output_directory+id_s+".ct") 58 command = shlex.split('Fold %s %s' % (os.path.join(syspath, temp.txt), os.path.join(output_directory, '%s.ct' % id_s)))
55 if flag == 1: 59 subprocess.call(command)
60 else:
56 if id_s in react: 61 if id_s in react:
57 f = file(syspath+"constraint.txt",'w') 62 fh = file(os.path.join(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 63 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]))): 64 for j in range(0, (len(react[id_s]))):
63 if react[id_s][j]!='NA': 65 if react[id_s][j]!='NA':
64 f.write(str(j+1)) 66 fh.write(str(j+1))
65 f.write('\t') 67 fh.write('\t')
66 f.write(str(react[id_s][j])) 68 fh.write(str(react[id_s][j]))
67 f.write('\n') 69 fh.write('\n')
68 #h.write(str(react[id_s][j])) #Output the reactivities 70 #h.write(str(react[id_s][j])) #Output the reactivities
69 #h.write('\t') 71 #h.write('\t')
70 f.close() 72 fh.close()
71 #h.write('\n') 73 #h.write('\n')
72 #h.write('\n') 74 #h.write('\n')
73 os.system("Fold "+syspath+"temp.txt"+" -sh"+" "+syspath+"constraint.txt"+" "+output_directory+id_s+".ct") 75 command = shlex.split("Fold %s -sh %s %s" % (os.path.join(syspath, "temp.txt"),
76 os.path.join(syspath, "constraint.txt"),
77 os.path.join(output_directory, "%s.ct" % id_s)))
78 subprocess.call(command)
74 else: 79 else:
75 print(id_s+" not in the data of react!") 80 print(id_s+" not in the data of react!")
76 os.system("draw "+output_directory+id_s+".ct "+output_directory+"/"+id_s+".ps") 81 command = shlex.split('draw %s.ct %s.ps' % (os.path.join(output_directory, id_s), os.path.join(output_directory, id_s)))
82 subprocess.call(command)
77 else: 83 else:
78 print(id_s+" not in the data of sequences!") 84 print(id_s+" not in the data of sequences!")
79 85
80 #Remove the unnecessary files 86 #Remove the unnecessary files
81 os.system("tar -zcvPf "+output_file+" "+output_directory+"/"+"*.* 2>"+output_directory+"log.txt") 87 tarball = tarfile.open(output_file, 'w:gz')
82 os.system("rm -f "+syspath+"temp.txt") 88 for filename in os.listdir(output_directory):
83 os.system("rm -r "+output_directory) 89 filepath = os.path.join(output_directory, filename)
84 if flag == 1: 90 print filepath
85 os.system("rm -f "+syspath+"constraint.txt") 91 tarball.add(filepath, arcname=filename)
86 # h.close() 92 print os.listdir(syspath)
87 93 print os.listdir(output_directory)
88 94 # tarball.add('%s.tif' % os.path.join(syspath, id_s), arcname='%s.tif' % id_s)
89 95 tarball.close()
90
91