Mercurial > repos > tyty > structurefold
diff predict/predict_RNAs.py @ 64:a1ce42d5258d draft
Uploaded
author | tyty |
---|---|
date | Tue, 18 Nov 2014 15:54:31 -0500 |
parents | |
children | 96a827962750 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/predict/predict_RNAs.py Tue Nov 18 15:54:31 2014 -0500 @@ -0,0 +1,91 @@ +#RNA structure prediction & Output and illustrate reactivities + +import sys +from parse_dis_pac import * +from read_file import * +from Bio import SeqIO +import os +from rtts_plot import * +import random +import string + + +id_file = sys.argv[1] +seq_file = sys.argv[2] +output_file = sys.argv[4] + + +flag = 0 +if sys.argv[3]!='None': #input reactivity file if provided + react_file = sys.argv[3] + react = parse_dist(react_file) + react = react[1] + flag = 1 + +syspath = os.getcwd() + +ids = read_t_file(id_file) +sequences = SeqIO.parse(seq_file, 'fasta') + + +seqs = {} +for seq in sequences: + seqs[seq.id] = seq.seq.tostring() + +if len(ids)>100: #setup a limit of the number of sequence to be predicted + print("Number of sequences exceeds limitation!") + sys.exit(0) + + +#predict RNA structures +output_directory = os.path.join(syspath, "output_files/") +os.makedirs(output_directory) +for i in range(len(ids)): + id_s = ids[i][0] + print(id_s) + #Put RNA sequence and reactivities into files + if id_s in seqs: + f = file(syspath+"temp.txt", 'w') + f.write('>'+id_s) + f.write('\n') + f.write(seqs[id_s]) + f.close() + if flag == 0: + os.system("Fold "+syspath+"temp.txt"+" "+output_directory+id_s+".ct") + if flag == 1: + if id_s in react: + f = file(syspath+"constraint.txt",'w') + make_plot(react[id_s],id_s,(output_directory)) #make a plot of the distribution of the reactivites of the input RNA + #h = file(syspath+"output_f/transcript_reactivities.txt", 'w') + #h.write(id_s) + #h.write('\n') + for j in range(0, (len(react[id_s]))): + if react[id_s][j]!='NA': + f.write(str(j+1)) + f.write('\t') + f.write(str(react[id_s][j])) + f.write('\n') + #h.write(str(react[id_s][j])) #Output the reactivities + #h.write('\t') + f.close() + #h.write('\n') + #h.write('\n') + os.system("Fold "+syspath+"temp.txt"+" -sh"+" "+syspath+"constraint.txt"+" "+output_directory+id_s+".ct") + else: + print(id_s+" not in the data of react!") + os.system("draw "+output_directory+id_s+".ct "+output_directory+"/"+id_s+".ps") + else: + print(id_s+" not in the data of sequences!") + +#Remove the unnecessary files +os.system("tar -zcvPf "+output_file+" "+output_directory+"/"+"*.* 2>"+output_directory+"log.txt") +os.system("rm -f "+syspath+"temp.txt") +os.system("rm -r "+output_directory) +if flag == 1: + os.system("rm -f "+syspath+"constraint.txt") + # h.close() + + + + +