Mercurial > repos > tyty > structurefold
view 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 source
#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()