Mercurial > repos > tyty > structurefold
comparison predict/predict_RNAs.py @ 67:96a827962750 draft
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 |