diff --git a/model.py b/model.py new file mode 100644 index 0000000..08c0f8d --- /dev/null +++ b/model.py @@ -0,0 +1,132 @@ +dim=10 +n=3 +fl=2*n+1 +maxbl=29 + +import tensorflow as tf + +def create_model(): + tf.reset_default_graph() + g=tf.Graph() + with g.as_default(): + synmg=tf.placeholder(dtype=tf.float64,shape=(None,2*n+1,2*n+1,2),name="Synteny_matrix_placeholder_Global") + synml=tf.placeholder(dtype=tf.float64,shape=(None,2*n+1,2*n+1,2),name="Synteny_matrix_placeholder_Local") + pfam=tf.placeholder(dtype=tf.float64,shape=(None,2*n+1,2*n+1,1),name="Pfam_matrix_placeholder") + bls=tf.placeholder(dtype=tf.float64,shape=(None,maxbl),name="Species_Branch_Length_Placeholder") + blhs=tf.placeholder(dtype=tf.float64,shape=(None,maxbl),name="Homology_Species_Branch_Length_Placeholder") + #gl=tf.placeholder(dtype=tf.float64,shape=(None,1),name="Mean_gene_length") + dps=tf.placeholder(dtype=tf.float64,shape=(None,1),name="mca_species_distance") + dphs=tf.placeholder(dtype=tf.float64,shape=(None,1),name="mca_homology_species_distance") + dis=tf.placeholder(dtype=tf.float64,shape=(None,1),name="total_distance") + lr=tf.placeholder(dtype=tf.float64,shape=(),name="learning_rate") + y=tf.placeholder(dtype=tf.int32,shape=(None),name="labels") + + lrs=tf.summary.scalar("Learning_Rate",lr) + + x=tf.concat([dps,dps-dphs,dis],1,name="Create_train_vector") + + print(synmg,"\n",synml,"\n",bls,"\n",blhs,"\n",dps,"\n",dphs,"\n",dis,"\n",x) + + reg_l2=tf.contrib.layers.l2_regularizer(0.001) + reg_l1 = tf.contrib.layers.l1_regularizer(scale=0.005, scope=None) + + def get_variable_by_shape(shape,name): + f=tf.get_variable(name,shape=shape,initializer=tf.glorot_uniform_initializer(),dtype=tf.float64,regularizer=reg_l2) + return f + + def create_synteny_aligner(name,synm): + with tf.variable_scope(name+"Synteny_Aligner",reuse=tf.AUTO_REUSE): + + fconv=get_variable_by_shape((2,2,2,dim),"fconv") + conv=tf.nn.conv2d(synm,fconv,(1,1,1,1),padding="VALID",name="Conv_aligner") + + fconv_1=get_variable_by_shape((2,2,dim,dim*2),"fconv_1") + conv_1=tf.nn.conv2d(conv,fconv_1,(1,1,1,1),padding="VALID",name="Conv_aligner_1") + + fxconv=get_variable_by_shape((fl,2,dim*2),"fxconv") + x_conv=tf.reshape(synm,(-1,fl*fl,2)) + x_conv=tf.nn.conv1d(x_conv,fxconv,stride=fl,padding="SAME",name="row_aligner") + + y_conv=tf.reshape(tf.transpose(synm,(0,2,1,3)),(-1,fl*fl,2)) + fyconv=get_variable_by_shape((fl,2,dim*2),"fyconv") + y_conv=tf.nn.conv1d(y_conv,fyconv,stride=fl,padding="SAME",name="column_aligner") + + + wconv=get_variable_by_shape((fl,fl,2,dim*2*10),"wconv") + w_conv=tf.nn.conv2d(synm,wconv,(1,1,1,1),padding="VALID",name="Global_Aligner_1") + + conv_1=tf.reshape(conv_1,(-1,25,dim*2)) + x_conv=tf.reshape(x_conv,(-1,fl,dim*2)) + y_conv=tf.reshape(y_conv,(-1,fl,dim*2)) + w_conv=tf.reshape(w_conv,(-1,10,dim*2)) + conv_final=tf.concat([conv_1,x_conv,y_conv,w_conv],1,name="Concatenate_All_Alignments") + return conv_final + + with tf.variable_scope("Pfam",reuse=tf.AUTO_REUSE): + wconv_pfam=get_variable_by_shape((fl,fl,1,dim*2*10),"wconv_pfam") + w_conv_pfam=tf.nn.conv2d(pfam,wconv_pfam,(1,1,1,1),padding="VALID",name="Global_Aligner_pfam") + w_conv_pfam=tf.reshape(w_conv_pfam,(-1,10,dim*2)) + + conv_final_g=create_synteny_aligner("Global_",synmg) + conv_final_l=create_synteny_aligner("Local_",synml) + final=tf.concat([conv_final_g,conv_final_l,w_conv_pfam],1) + #final=conv_final_l + + with tf.variable_scope("Combine_Renormalize",reuse=tf.AUTO_REUSE): + bl=tf.concat([bls,blhs],1) + #bl=bls-blhs + theta_bl=get_variable_by_shape((maxbl*2,1),"theta_bl") + theta_bl=tf.matmul(bl,theta_bl) + x=tf.concat([x,theta_bl],1) + theta=get_variable_by_shape((4,108),"theta") + bias=get_variable_by_shape((1,108),"b") + theta_2=tf.matmul(x,theta)+bias + theta_2=tf.reshape(theta_2,(-1,108,1)) + theta_2=tf.tile(theta_2,[1,1,dim*2]) + final=final*theta_2 + print(final) + + flat=tf.layers.flatten(final) + + zero=tf.constant(0.0,dtype=tf.float64) + diff=dps-dphs + print(diff) + diff_2=tf.cast(tf.equal(diff,zero),tf.float64) + print(diff_2) + diff=tf.tile(diff_2,[1,dim*10]) + + flat=tf.concat([flat,diff],1) + print(flat) + #dense=tf.layers.dense(flat,2048,kernel_regularizer=reg_l2,bias_regularizer=reg_l2) + #dense_2=tf.layers.dense(dense,1024,kernel_regularizer=reg_l2,bias_regularizer=reg_l2) + dense_3=tf.layers.dense(flat,512,kernel_regularizer=reg_l2,bias_regularizer=reg_l2) + + logits_pred=tf.layers.dense(dense_3,3,name="Predictions") + print(logits_pred) + entropy=tf.nn.sparse_softmax_cross_entropy_with_logits(logits=logits_pred,labels=y) + print(entropy) + + #weights = tf.trainable_variables() # all vars of your graph + #regl1 = tf.contrib.layers.apply_regularization(reg_l1, weights) + reg_losses = tf.get_collection(tf.GraphKeys.REGULARIZATION_LOSSES) + reg_constant =0.00000001 + loss=tf.reduce_mean(entropy)+reg_constant * sum(reg_losses) + #loss=tf.reduce_mean(entropy) + optimizer=tf.train.RMSPropOptimizer(lr) + #optimizer=tf.train.AdamOptimizer() + losses=tf.summary.scalar("Loss",loss) + t_op=optimizer.minimize(loss) + acc=tf.math.in_top_k(tf.cast(logits_pred,tf.float32),y,1) + accuracy=tf.reduce_mean(tf.cast(acc,tf.float32)) + accs=tf.summary.scalar("Accuracy",accuracy) + summary=tf.summary.merge_all() + init=tf.global_variables_initializer() + saver=tf.train.Saver() + + for node in (synmg,synml,pfam,bls,blhs,dps,dphs,dis,lr,y): + g.add_to_collection("input_nodes",node) + + for node in (loss,t_op,accuracy,init,summary): + g.add_to_collection("output_nodes",node) + + return g,saver diff --git a/prepare_synteny_matrix.py b/prepare_synteny_matrix.py index 9de052d..45fc93a 100644 --- a/prepare_synteny_matrix.py +++ b/prepare_synteny_matrix.py @@ -1,69 +1,90 @@ -import pandas as pd import numpy as np -import os -import sys -import progressbar +import pandas as pd import json -import gc -from read_data import read_data_homology +import os +import sys +import pickle +from select_data import read_db_homology +from threads import Procerssrunner from read_get_gene_seq import read_gene_sequences -from create_synteny_matrix import synteny_matrix -from access_data_rest import update_rest -from select import select,create_map_reverse +from access_data_rest import update_rest,update_rest_protein +from process_negative import write_fasta -if not os.path.isdir("processed/synteny_matrices"): - os.mkdir("processed/synteny_matrices") +def read_data_synteny(nop,name): + smg=[] + sml=[] + indexes=[] + for i in range(nop): + try: + with open("temp_"+name+"/thread_"+str(i+1)+"_smg.temp","rb") as file: + smg=smg+pickle.load(file) + with open("temp_"+name+"/thread_"+str(i+1)+"_sml.temp","rb") as file: + sml=sml+pickle.load(file) + with open("temp_"+name+"/thread_"+str(i+1)+"_indexes.temp","rb") as file: + indexes=indexes+pickle.load(file) + except Exception as e: + print("Problem with thread",i+1,"detected for",name,e) + continue + print(len(indexes)) + return smg,sml,indexes -arg=sys.argv -arg=arg[1:] +def load_neighbor_genes(): + with open("processed/neighbor_genes.json","r") as file: + lsy=dict(json.load(file)) + print(len(lsy)) + print("Neighbor Genes Loaded") + return lsy -nos=int(arg[0]) +def read_data_homology(dirname): + lf=os.listdir(dirname) + if len(lf)==0: + print("No Files in the Directory!!!!!!!") + sys.exit(1) + a_h=[] + d_h=[] + for x in lf: + df,n=read_db_homology(dirname,x) + n=n.split()[0] + try: + indexes=np.load("processed/"+n+"_selected_indexes.npy") + except: + print("Incomplete data for:",n) + df=df.loc[indexes] + a_h.append(df) + d_h.append(n) + return a_h,d_h -with open("dist_matrix","r") as file: - matrix=file.readlines() -matrix=[x.split("\t") for x in matrix] -matrix=[[float(y) for y in x] for x in matrix] -matrix=np.array(matrix) -with open("sp_names","r") as file: - dname=file.readlines() -dname=[x.split("\n")[0] for x in dname] -spnmap,nspmap=create_map_reverse(dname) - -lsy={} -a_h,d_h=read_data_homology("data_homology") -print(d_h) -d_h=list(d_h.keys()) -d_h=[x.split()[0] for x in d_h] -print("Homology Data Read") -for i in range(len(a_h)): - a_h[i]=select(a_h[i],nos,matrix,spnmap,nspmap,d_h[i]) - assert(len(a_h[i])==nos) -gc.collect() -print("Data Selected") - -with open("processed/neighbor_genes.json","r") as file: - lsy=dict(json.load(file)) -print(len(lsy)) -print("Neighbor Genes Loaded") - -gene_sequences={} -gene_sequences=read_gene_sequences(a_h,lsy,"geneseq","gene_sequences") -print("Gene Sequences Loaded") - -print("Going to update not found sequences:") -gene_sequences=update_rest(gene_sequences) - -n=3 -ndir="processed/synteny_matrices/" -nf1="synteny_matrices_global" -nf2="synteny_matrices_local" -nf3="indexes" -for i in range(len(a_h)): - df=a_h[i] - print(len(df)) - synteny_matrices_global,synteny_matrices_local,indexes=synteny_matrix(gene_sequences,df,lsy,n,0) - np.save(ndir+str(d_h[i])+"_"+nf1,synteny_matrices_global) - np.save(ndir+str(d_h[i])+"_"+nf2,synteny_matrices_local) - np.save(ndir+str(d_h[i])+"_"+nf3,indexes) - print(len(indexes)) -print("Synteny Matrices Created Successfully :)") +def main(): + arg=sys.argv + nop=int(arg[-1]) + n=3 + a_h,d_h=read_data_homology("data_homology") + print("Data Read") + lsy=load_neighbor_genes() + gene_sequences=read_gene_sequences(a_h,lsy,"geneseq","gene_seq_positive") + gene_sequences=update_rest(gene_sequences,"gene_seq_positive") + print("Gene Sequences Loaded.") + if not os.path.isdir("processed/synteny_matrices"): + os.mkdir("processed/synteny_matrices") + ndir="processed/synteny_matrices/" + nf1="synteny_matrices_global" + nf2="synteny_matrices_local" + nf3="indexes" + for i in range(len(a_h)): + df=a_h[i] + part=len(df)//nop + pr=Procerssrunner() + pr.start_processes(nop,df,gene_sequences,lsy,part,n,d_h[i]) + smg,sml,indexes=read_data_synteny(nop,d_h[i]) + print(len(indexes)) + np.save(ndir+str(d_h[i])+"_"+nf1,smg) + np.save(ndir+str(d_h[i])+"_"+nf2,sml) + np.save(ndir+str(d_h[i])+"_"+nf3,indexes) + a_h[i]=df.loc[indexes] + print("Synteny Matrices Created Successfully :)") + protein_sequences=read_gene_sequences(a_h,lsy,"pro_seq","pro_seq_positive") + protein_sequences=update_rest_protein(protein_sequences,"pro_seq_positive") + write_fasta(protein_sequences,"protein_seq_positive") + print("Protein Sequences Loaded.") +if __name__=="__main__": + main() diff --git a/train.py b/train.py new file mode 100644 index 0000000..a1dfb9a --- /dev/null +++ b/train.py @@ -0,0 +1,192 @@ +import pickle +import numpy as np +import sys +import tensorflow as tf +from model import create_model + +def create_branch_length_padding(bl): + maxlen=0 + for x in bl: + if len(x)>maxlen: + maxlen=len(x) + + for x in range(len(bl)): + temp=bl[x] + for i in range(len(temp),maxlen): + temp=np.append(temp,[0]) + bl[x]=temp + return bl + +def train(train_synteny_matrices_global,train_synteny_matrices_local,train_pfam_matrices,train_branch_length_species,train_branch_length_homology_species,train_dist_p_s,train_dist_p_hs,train_distance,train_labels,v,num_epochs,learning_rate,decay,size_train,batch_size,model_name): + print("Going to train model {} for:\n Batch Size:{} \n Learning Rate:{} \n Decay:{}\n On {} Samples".format(v,batch_size,learning_rate,decay,len(train_synteny_matrices_global))) + graph,saver=create_model() + synmg,synml,pfam,bls,blhs,dps,dphs,dis,lr,y=graph.get_collection("input_nodes") + loss,t_op,accuracy,init,summary=graph.get_collection("output_nodes") + with tf.Session(graph=graph) as sess: + writer = tf.summary.FileWriter('./'+model_name+'_v'+str(v), sess.graph) + sess.run(init) + learn=learning_rate + for j in range(num_epochs): + for i in range(size_train//batch_size): + feed_dict={ + synmg:train_synteny_matrices_global[i*batch_size:(i+1)*batch_size], + synml:train_synteny_matrices_local[i*batch_size:(i+1)*batch_size], + pfam:train_pfam_matrices[i*batch_size:(i+1)*batch_size].reshape((batch_size,7,7,1)), + bls:train_branch_length_species[i*batch_size:(i+1)*batch_size], + blhs:train_branch_length_homology_species[i*batch_size:(i+1)*batch_size], + dps:train_dist_p_s[i*batch_size:(i+1)*batch_size].reshape((batch_size,1)), + dphs:train_dist_p_hs[i*batch_size:(i+1)*batch_size].reshape((batch_size,1)), + dis:train_distance[i*batch_size:(i+1)*batch_size].reshape((batch_size,1)), + lr:learn, + y:train_labels[i*batch_size:(i+1)*batch_size] + } + sess.run(t_op,feed_dict=feed_dict) + + test_dict={ + synmg:train_synteny_matrices_global[size_train:], + synml:train_synteny_matrices_local[size_train:], + pfam:train_pfam_matrices[size_train:].reshape((len(train_synteny_matrices_global)-size_train,7,7,1)), + bls:train_branch_length_species[size_train:], + blhs:train_branch_length_homology_species[size_train:], + dps:train_dist_p_s[size_train:].reshape((len(train_synteny_matrices_global)-size_train,1)), + dphs:train_dist_p_hs[size_train:].reshape((len(train_synteny_matrices_global)-size_train,1)), + dis:train_distance[size_train:].reshape((len(train_synteny_matrices_global)-size_train,1)), + y:train_labels[size_train:], + lr:learn + } + accuracy_test,loss_test,summary_write=sess.run([accuracy,loss,summary],feed_dict=test_dict) + writer.add_summary(summary_write,i+1) + print("Epoch:{} Test Accuracy:{} Test Loss:{}".format(j+1,accuracy_test*100,loss_test)) + learn*=decay + saver.save(sess,model_name+"_v"+str(v)+"/model.ckpt") + writer.close() + +def read_positive(len_p,bls,blhs,dis,dps,dphs,sml,smg,pfam,label): + rowsh=[] + with open("dataset","rb") as file: + rowsh=pickle.load(file) + shi=np.random.permutation(len(rowsh)) + rows_shuffled=[] + for i in range(len(shi)): + rows_shuffled.append(rowsh[shi[i]]) + rowsh=rows_shuffled + spco={} + spcp={} + for row in rowsh: + if row["species"] not in spco: + spco[row["species"]]=0 + spcp[row["species"]]=0 + + maxcount_o=int((len_p)*0.3/14) + maxcount_p=int((len_p)*0.7/14) + for row in rowsh: + if row["label"]==2: + continue + if row["label"]==1 and spco[row["species"]]>maxcount_o: + continue + if row["label"]==0 and spcp[row["species"]]>maxcount_p: + continue + bls.append(np.array(row["bls"])) + blhs.append(np.array(row["blhs"])) + dis.append(row["dis"]) + dps.append(row["dps"]) + dphs.append(row["dphs"]) + sml.append(row["local_alignment_matrix"]) + smg.append(row["global_alignment_matrix"]) + pfam.append(row["pfam_matrix"]) + label.append(row["label"]) + if row["label"]==1: + spco[row["species"]]+=1 + if row["label"]==0: + spcp[row["species"]]+=1 + + +def read_negative(len_n,bls,blhs,dis,dps,dphs,sml,smg,pfam,label): + rows=[] + with open("dataset","rb") as file: + rows=pickle.load(file) + rows=[row for row in rows if row["label"]==2] + shi=np.random.permutation(len(rows)) + rows_shuffled=[] + for i in range(len(shi)): + rows_shuffled.append(rows[shi[i]]) + rows=rows_shuffled + rows=rows[:len_n] + for row in rows: + bls.append(np.array(row["bls"])) + blhs.append(np.array(row["blhs"])) + dis.append(row["dis"]) + dps.append(row["dps"]) + dphs.append(row["dphs"]) + sml.append(row["local_alignment_matrix"]) + smg.append(row["global_alignment_matrix"]) + label.append(row["label"]) + pfam.append(row["pfam_matrix"]) + +def train_models(model_name,start,end,num_epochs,learn_rate,decay,size_train,batch_size): + k=1 + for i in range(start//10,end//10+1): + bls=[] + blhs=[] + dis=[] + dps=[] + dphs=[] + sml=[] + smg=[] + pfam=[] + label=[] + portion=float(i/10) + len_n=int(size_train*portion) + len_p=int(size_train*(1-portion)) + read_positive(len_p,bls,blhs,dis,dps,dphs,sml,smg,pfam,label) + read_negative(len_n,bls,blhs,dis,dps,dphs,sml,smg,pfam,label) + bls=create_branch_length_padding(bls) + blhs=create_branch_length_padding(blhs) + bls=np.array(bls) + print(bls.shape) + blhs=np.array(blhs) + print(blhs.shape) + dis=np.array(dis) + print(dis.shape) + dps=np.array(dps) + print(dps.shape) + dphs=np.array(dphs) + print(dphs.shape) + sml=np.array(sml) + print(sml.shape) + smg=np.array(smg) + print(smg.shape) + pfam=np.array(pfam) + print(pfam.shape) + label=np.array(label) + print(label.shape) + shi=np.random.permutation(len(label)) + labels=label[shi] + bls=bls[shi] + blhs=blhs[shi] + dis=dis[shi] + dps=dps[shi] + dphs=dphs[shi] + sml=sml[shi] + smg=smg[shi] + pfam=pfam[shi] + train(smg,sml,pfam,bls,blhs,dps,dphs,dis,labels,k,num_epochs,learn_rate,decay,int(0.9*size_train),batch_size,model_name) + k+=1 + +def main(): + arg=sys.argv + model_name=arg[-8] + start_p=int(arg[-7]) + end_p=int(arg[-6]) + num_epochs=int(arg[-5]) + learn_rate=float(arg[-4]) + decay=float(arg[-3]) + size_train=float(arg[-2]) + batch_size=int(arg[-1]) + train_models(model_name,start_p,end_p,num_epochs,learn_rate,decay,size_train,batch_size) + + +if __name__=="__main__": + main() + +