diff --git a/create_synteny_matrix.py b/create_synteny_matrix.py new file mode 100644 index 0000000..eff7d6d --- /dev/null +++ b/create_synteny_matrix.py @@ -0,0 +1,99 @@ +import numpy as np +import requests +import edlib as ed +import pandas as pd +import time +import sys +from skbio.alignment import local_pairwise_align_ssw +from skbio import DNA,TabularMSA,RNA + +def update(gene_seq,gene): + server = "https://rest.ensembl.org" + ext = "/sequence/id/"+str(gene)+"?" + + r = requests.get(server+ext, headers={ "Content-Type" : "text/plain"}) + + if not r.ok: + r.raise_for_status() + sys.exit() + + gene_seq[gene]=str(r.text) + +def create_synteny_matrix_mul(gene_seq,g1,g2,n): + for gene in g1: + if gene=="NULL_GENE": + continue + try: + temp=gene_seq[gene] + except: + print("Updating gene sequences for gene:",gene) + update(gene_seq,gene) + for gene in g2: + if gene=="NULL_GENE": + continue + try: + temp=gene_seq[gene] + except: + print("Updating gene sequences for gene:",gene) + update(gene_seq,gene) + #print(n) + sm=np.zeros((n,n,2)) + sml=np.zeros((n,n,2)) + for i in range(n): + if g1[i]=="NULL_GENE": + continue + for j in range(n): + if g2[j]=="NULL_GENE": + continue + norm_len=(len(gene_seq[g1[i]])+len(gene_seq[g2[j]])) + result = ed.align(gene_seq[g1[i]],gene_seq[g2[j]], mode="NW", task="distance") + sm[i][j][0]=result["editDistance"]/(norm_len) + result = ed.align(gene_seq[g1[i]],gene_seq[g2[j]][::-1], mode="NW", task="distance") + sm[i][j][1]=result["editDistance"]/(norm_len) + _,result,_=local_pairwise_align_ssw(DNA(gene_seq[g1[i]]),DNA(gene_seq[g2[j]])) + sml[i][j][0]=result/(norm_len) + _,result,_=local_pairwise_align_ssw(DNA(gene_seq[g1[i]]),DNA(gene_seq[g2[j]][::-1])) + sml[i][j][1]=result/(norm_len) + return sm,sml + +def synteny_matrix(gene_seq,hdf,lsy,n): + sg=[] + sl=[] + t=0 + ind=[] + start=time.time() + for index,row in hdf.iterrows(): + g1=str(row["gene_stable_id"]) + g2=str(row["homology_gene_stable_id"]) + x=[] + y=[] + try: + temp=lsy[g1] + temp=lsy[g2] + except: + continue + for i in range(len(lsy[g1]['b'])-1,-1,-1): + x.append(lsy[g1]['b'][i]) + x.append(g1) + for k in lsy[g1]['f']: + x.append(k) + + for i in range(len(lsy[g2]['b'])-1,-1,-1): + y.append(lsy[g2]['b'][i]) + y.append(g2) + for k in lsy[g2]['f']: + y.append(k) + + assert(len(x)==len(y)) + assert(len(x)==(2*n+1)) + smgtemp,smltemp=create_synteny_matrix_mul(gene_seq,x,y,2*n+1) + sg.append(smgtemp) + sl.append(smltemp) + ind.append(index) + t+=1 + if t==5: + break + end=time.time() + print("Time Taken:",end-start) + print("Average Time:",(end-start)/len(sg)) + return np.array(sg),np.array(sl),np.array(ind) diff --git a/create_train_data.py b/create_train_data.py new file mode 100644 index 0000000..33a96e5 --- /dev/null +++ b/create_train_data.py @@ -0,0 +1,104 @@ +import pandas as pd +import numpy as np +def create_branch_length_padding(bl): + maxlen=0 + for x in bl: + if len(x)>maxlen: + maxlen=len(x) + + for x in bl: + for i in range(len(x),maxlen): + x.append(0) + +def train_data(indexes,synteny_matrices,df,branch_length_species,branch_length_homology_species,distance,dist_p_s,dist_p_hs,gene_sequences): + + homology_type_counts=dict(df.homology_type.value_counts()) + homology_species_counts=dict(df.homology_species.value_counts()) + + #limits for each species and homology type in the trainig data so that the dataset is balanced + max_species_count=10000 + max_homology_type_count=25000 + for species in homology_species_counts: + homology_species_counts[species]=0 + + for homology_type in homology_type_counts: + homology_type_counts[homology_type]=0 + + labels=dict(ortholog_one2one=0,other_paralog=1,ortholog_one2many=2,ortholog_many2many=3,within_species_paralog=4) + + train_data_dataframe=pd.DataFrame() + train_branch_length_species=[] + train_branch_length_homology_species=[] + train_distance=[] + train_dist_p_s=[] + train_dist_p_hs=[] + train_labels=[] + train_indexes=[] + train_mean_gene_length=[] + for i in range(len(indexes)): + row=df.loc[indexes[i]] + if homology_species_counts[row["homology_species"]]>=max_species_count and row["homology_type"]!="within_species_paralog": + continue + if homology_type_counts[row["homology_type"]]>=max_homology_type_count: + continue + + train_branch_length_species.append(branch_length_species[i]) + train_branch_length_homology_species.append(branch_length_homology_species[i]) + train_distance.append(distance[i]) + train_dist_p_s.append(dist_p_s[i]) + train_dist_p_hs.append(dist_p_hs[i]) + train_data_dataframe=train_data_dataframe.append(row) + train_indexes.append(i) + train_labels.append(labels[row["homology_type"]]) + homology_species_counts[row["homology_species"]]+=1 + homology_type_counts[row["homology_type"]]+=1 + train_mean_gene_length.append((len(gene_sequences[row["gene_stable_id"]])+len(gene_sequences[row["homology_gene_stable_id"]]))/2) + + train_synteny_matrices=synteny_matrices[train_indexes] + + create_branch_length_padding(train_branch_length_species) + train_branch_length_species=np.array(train_branch_length_species) + create_branch_length_padding(train_branch_length_homology_species) + train_branch_length_homology_species=np.array(train_branch_length_homology_species) + + #renormalize the train_mean_gene_length by (x-mean)/std + train_mean_gene_length=(train_mean_gene_length-np.mean(train_mean_gene_length))/np.std(train_mean_gene_length) + + #create a random array of permutations to shuffle the indices + shi=np.random.permutation(len(train_labels)) + + train_branch_length_species=train_branch_length_species[shi] + print(train_branch_length_species.shape) + + train_branch_length_homology_species=train_branch_length_homology_species[shi] + print(train_branch_length_homology_species.shape) + + train_dist_p_s=np.array(train_dist_p_s) + train_dist_p_s=train_dist_p_s[shi] + print(train_dist_p_s.shape) + + train_dist_p_hs=np.array(train_dist_p_hs) + train_dist_p_hs=train_dist_p_hs[shi] + print(train_dist_p_hs.shape) + + train_synteny_matrices=train_synteny_matrices[shi] + print(train_synteny_matrices.shape) + + train_indexes=np.array(train_indexes) + train_indexes=train_indexes[shi] + print(train_indexes.shape) + + train_labels=np.array(train_labels) + train_labels=train_labels[shi] + print(train_labels.shape) + + train_mean_gene_length=np.array(train_mean_gene_length) + train_mean_gene_length=train_mean_gene_length[shi] + print(train_mean_gene_length.shape) + + train_distance=np.array(train_distance) + train_distance=train_distance[shi] + train_distance=(train_distance-np.mean(train_distance))/np.std(train_distance) + train_distance.shape + + return train_synteny_matrices,train_branch_length_species,train_branch_length_homology_species,train_mean_gene_length,train_dist_p_s,train_dist_p_hs,train_distance,train_labels diff --git a/ftpg.py b/ftpg.py new file mode 100644 index 0000000..3139d00 --- /dev/null +++ b/ftpg.py @@ -0,0 +1,47 @@ +from ftplib import FTP +from req_data import get_data_file + + +#This wil download all the fasta files for the coding sequences. To change the directory, change the argument in the get_data_file argument. +host ="ftp.ensembl.org" +user = "anonymous" +password = "" + +print("Connecting to {}".format(host)) +ftp = FTP(host) +ftp.login(user, password) +print("Connected to {}".format(host)) +l=ftp.nlst("/pub/release-96/fasta") +lt=[] +for x in l: + y=ftp.nlst(x+"/cds") + for z in y: + if z.endswith(".cds.all.fa.gz"): + lt.append(z) + +base_link="ftp://ftp.ensembl.org" +with open("seq_link.txt","w") as file: + for x in lt: + file.write(base_link+x) + file.write("\n") + +l=ftp.nlst("/pub/release-96/gtf") +lt=[] +for x in l: + y=ftp.nlst(x) + for z in y: + if z.endswith(".96.gtf.gz"): + lt.append(z) + +base_link="ftp://ftp.ensembl.org" +with open("gtf_link.txt","w") as file: + for x in lt: + file.write(base_link+x) + file.write("\n") + +ch=input("Do you want to download the data?[y/n]") +if ch=='y': + print("Downloading Data.................") + get_data_file("gtf_link.txt","data") + get_data_file("seq_link.txt","geneseq") + print("Download Complete.................") diff --git a/get_data.py b/get_data.py index 2691d81..a2ee986 100644 --- a/get_data.py +++ b/get_data.py @@ -4,7 +4,11 @@ from req_data import get_data_file,download_data from read_data import read_data_genome,read_data_homology from process_data import list_dict_genomes -def get_data_genome(arg,dir,a,d,ld,ldg): +def get_data_genome(arg,dir): + a=[] + d={} + ld=[] + ldg=[] if arg[0]=='-d': if arg[4]=="-r": c=0 @@ -28,7 +32,9 @@ def get_data_genome(arg,dir,a,d,ld,ldg): return ld,ldg,a,d -def get_data_homology(arg,dir,a_h,d_h): +def get_data_homology(arg,dir): + a_h=[] + d_h={} if arg[2]=="-l": if not os.path.exists(dir): os.mkdir(dir) @@ -44,7 +50,7 @@ def get_data_homology(arg,dir,a_h,d_h): return a_h,d_h if arg[4]=="-r": - a_h,d_h=read_data_homology(dir,a_h,d_h) + a_h,d_h=read_data_homology(dir) assert(len(a_h)==len(d_h)) return a_h,d_h diff --git a/main.py b/main.py index e2770d5..47d56e2 100644 --- a/main.py +++ b/main.py @@ -1,7 +1,14 @@ import sys - +import numpy as np +import pandas as pd +import json from get_data import get_data_homology,get_data_genome from process_data import create_data_homology_ls +from read_get_gene_seq import read_gene_sequences +from create_synteny_matrix import synteny_matrix +from tree_data import create_tree_data +from create_train_data import train_data +from train import train arg=sys.argv arg=arg[1:] @@ -10,33 +17,37 @@ if len(arg)!=5: print("No. of arguments more or less. Please check") sys.exit(1) -a=[] -d={} -a_h=[] -d_h={} -ld=[] -ldg=[] dir_g="data" +ld,ldg,a,d=get_data_genome(arg,dir_g) -ld,ldg,a,d=get_data_genome(arg,dir_g,a,d,ld,ldg) -#print(a[0][0:10],"\n",d,"\n",ld[0][0:10],"\n") dir_hom="data_homology" +a_h,d_h=get_data_homology(arg,dir_hom) -a_h,d_h=get_data_homology(arg,dir_hom,a_h,d_h) -#print(a_h[0][0:10],"\n",d_h) +if arg[-1]=="-d": + sys.exit(1) + +print("Data Read") n=2 #no. of numbers neighbors +save_after=50000 #to save data after n steps +lsy=create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,save_after) +print(len(lsy)) -lsy,lcmap=create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg) -""" -print(lsy,"\n",lcmap) -lt=ldg[0] -for x in lsy: +print("Neighbor Genes Found") - xr=lsy[x]['f'] - xl=lsy[x]['b'] - for g in range(len(xl)-1,-1,-1): - print(a[0].iloc[lt[xl[g]],[3,4]]) - print("------------------\n",a[0].iloc[lt[x],[3,4]],"\n-------------------") - for g in xr: - print(a[0].iloc[lt[g],[3,4]])""" +gene_sequences=read_gene_sequences(a_h[0],lsy,"geneseq","gene_sequences") + +synteny_matrices,indexes=synteny_matrix(gene_sequences,a_h[0],lsy,n) +print("Synteny Matrices are created successfully\n",len(indexes),"\n",len(synteny_matrices)) +np.save("synteny_matrices",synteny_matrices) +np.save("indexes",indexes) + +with open("gene_seq_updated.json","w") as file: + json.dump(gene_sequences,file) + +df=a_h[0].loc[indexes] +branch_length_species,branch_length_homology_species,distance,dist_p_s,dist_p_hs=create_tree_data("species_tree.tree",df) + +train_synteny_matrices,train_branch_length_species,train_branch_length_homology_species,train_mean_gene_length,train_dist_p_s,train_dist_p_hs,train_distance,train_labels=train_data(indexes,synteny_matrices,df,branch_length_species,branch_length_homology_species,distance,dist_p_s,dist_p_hs,gene_sequences) + +train(train_synteny_matrices,train_branch_length_species,train_branch_length_homology_species,train_mean_gene_length,train_dist_p_s,train_dist_p_hs,train_distance,train_labels) diff --git a/model.py b/model.py new file mode 100644 index 0000000..e547bf4 --- /dev/null +++ b/model.py @@ -0,0 +1,117 @@ +dim=30 + +import tensorflow as tf + +def create_model(): + tf.reset_default_graph() + g=tf.Graph() + with g.as_default(): + synm=tf.placeholder(dtype=tf.float64,shape=(None,5,5,2),name="Synteny_matrix_placeholder") + bls=tf.placeholder(dtype=tf.float64,shape=(None,21),name="Species_Branch_Length_Placeholder") + blhs=tf.placeholder(dtype=tf.float64,shape=(None,21),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_with_default(0.1,(),"learning_rate") + y=tf.placeholder(dtype=tf.int32,shape=(None),name="labels") + + #reshape the variables + #gl=tf.expand_dims(gl,-1) + #dps=tf.expand_dims(dps,-1) + #dphs=tf.expand_dims(dphs,-1) + #dis=tf.expand_dims(dis,-1) + + x=tf.concat([gl,dps,dphs,dis],1,name="Create_train_vector") + + print(synm,"\n",bls,"\n",blhs,"\n",gl,"\n",dps,"\n",dphs,"\n",dis,"\n",x) + + reg=tf.contrib.layers.l2_regularizer(0.01) + + def get_variable_by_shape(shape,name): + f=tf.get_variable(name,shape=shape,initializer=tf.glorot_uniform_initializer(),dtype=tf.float64,regularizer=reg) + return f + + with tf.variable_scope("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") + + print(conv) + print(conv_1) + #fxconv=get_variable_by_shape((1,5,2,dim),"fxconv") + fxconv=get_variable_by_shape((5,2,dim*2),"fxconv") + x_conv=tf.reshape(synm,(-1,25,2)) + x_conv=tf.nn.conv1d(x_conv,fxconv,stride=5,padding="SAME",name="row_aligner") + print(x_conv) + #fyconv=get_variable_by_shape((5,1,2,dim),"fyconv") + #y_conv=tf.nn.conv2d(synm,fyconv,(1,1,1,1),padding="SAME",name="column_aligner") + y_conv=tf.reshape(tf.transpose(synm,(0,2,1,3)),(-1,25,2)) + print(y_conv) + fyconv=get_variable_by_shape((5,2,dim*2),"fyconv") + y_conv=tf.nn.conv1d(y_conv,fyconv,stride=5,padding="SAME",name="column_aligner") + + print(y_conv) + conv_1=tf.reshape(conv_1,(-1,9,dim*2)) + print(conv) + x_conv=tf.reshape(x_conv,(-1,5,dim*2)) + print(x_conv) + y_conv=tf.reshape(y_conv,(-1,5,dim*2)) + print(y_conv) + conv_final=tf.concat([conv_1,x_conv,y_conv],1,name="Concatenate_All_Alignments") + print(conv_final) + + W=get_variable_by_shape((5,5,2),"W") + print(W) + W=W*synm + print(W) + #W=tf.reduce_sum(W,[1,2,3]) + #W=tf.reshape(W,(-1,1)) + #print(W) + + + with tf.variable_scope("Combine_Renormalize",reuse=tf.AUTO_REUSE): + bl=tf.concat([bls,blhs],1) + theta_bl=get_variable_by_shape((42,1),"theta_bl") + theta_bl=tf.matmul(bl,theta_bl) + x=tf.concat([x,theta_bl],1) + theta=get_variable_by_shape((5,19),"theta") + bias=get_variable_by_shape((1,19),"b") + theta_2=tf.matmul(x,theta)+bias + print(theta_2) + theta_2=tf.reshape(theta_2,(-1,19,1)) + theta_2=tf.tile(theta_2,[1,1,dim*2]) + print(theta_2) + conv_final=conv_final*theta_2 + print(conv_final) + + flat=tf.layers.flatten(conv_final) + flat_w=tf.layers.flatten(W) + flat=tf.concat([flat,flat_w],1) + print(flat) + dense=tf.layers.dense(flat,2048,kernel_regularizer=reg,bias_regularizer=reg) + dense_2=tf.layers.dense(dense,1024,kernel_regularizer=reg,bias_regularizer=reg) + dense_3=tf.layers.dense(dense_2,512,kernel_regularizer=reg,bias_regularizer=reg) + logits_pred=tf.layers.dense(dense_3,5,name="Predictions") + entropy=tf.nn.sparse_softmax_cross_entropy_with_logits(logits=logits_pred,labels=y) + print(entropy) + reg_losses = tf.get_collection(tf.GraphKeys.REGULARIZATION_LOSSES) + reg_constant = 0.0001 + loss=tf.reduce_mean(entropy)+reg_constant * sum(reg_losses) + optimizer=tf.train.RMSPropOptimizer(lr) + 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)) + init=tf.global_variables_initializer() + saver=tf.train.Saver() + + for node in (synm,bls,blhs,gl,dps,dphs,dis,lr,y): + g.add_to_collection("input_nodes",node) + + for node in (loss,t_op,accuracy,init): + g.add_to_collection("output_nodes",node) + + return g,saver diff --git a/prepare_data.py b/prepare_data.py new file mode 100644 index 0000000..3b991cd --- /dev/null +++ b/prepare_data.py @@ -0,0 +1,87 @@ +import sys +import numpy as np +import pandas as pd +import json +import os +from get_data import get_data_homology,get_data_genome +from process_data import create_data_homology_ls +from read_get_gene_seq import read_gene_sequences +from create_synteny_matrix import synteny_matrix +from tree_data import create_tree_data +from prepare_train_data import train_data + +if not os.path.exists("processed"): + os.mkdir("processed") + +arg=sys.argv +arg=arg[1:] + +enable_break=0 + +if arg[-1]=="-test": + enable_break=1 + +arg=arg[:-1] + +if len(arg)!=5: + print("No. of arguments more or less. Please check") + sys.exit(1) + +dir_g="data" +ld,ldg,a,d=get_data_genome(arg,dir_g) + +dir_hom="data_homology" +a_h,d_h=get_data_homology(arg,dir_hom) + +if arg[-1]=="-d": + sys.exit(1) + +print("Data Read") + +n=3 #no. of numbers neighbors +save_after=0 #to save data after n steps + +if enable_break==1: + save_after=3 +else: + save_after=1000000 + +lsy=create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,save_after,enable_break) +print(len(lsy)) + +print("Neighbor Genes Found") + +if enable_break==1: + gene_sequences=read_gene_sequences(a_h[0][0:10],lsy,"geneseq","gene_sequences") +else: + gene_sequences=read_gene_sequences(a_h[0],lsy,"geneseq","gene_sequences") + +if enable_break==1: + synteny_matrices_global,synteny_matrices_local,indexes=synteny_matrix(gene_sequences,a_h[0][0:10],lsy,n) +else: + synteny_matrices_global,synteny_matrices_local,indexes=synteny_matrix(gene_sequences,a_h[0],lsy,n) + +print("Synteny Matrices are created successfully\n",len(indexes),"\n",len(synteny_matrices_global)) +np.save("processed/synteny_matrices_global",synteny_matrices_global) +np.save("processed/synteny_matrices_local",synteny_matrices_local) +np.save("processed/indexes",indexes) + +with open("processed/gene_seq_updated.json","w") as file: + json.dump(gene_sequences,file) + +df=a_h[0].loc[indexes] +branch_length_species,branch_length_homology_species,distance,dist_p_s,dist_p_hs=create_tree_data("species_tree.tree",df) + +train_synteny_matrices_global,train_synteny_matrices_local,train_branch_length_species,train_branch_length_homology_species,train_mean_gene_length,train_dist_p_s,train_dist_p_hs,train_distance,train_labels=train_data(indexes,synteny_matrices_global,synteny_matrices_local,df,branch_length_species,branch_length_homology_species,distance,dist_p_s,dist_p_hs,gene_sequences) + +np.save("processed/train_synteny_matrices_global",train_synteny_matrices_global) +np.save("processed/train_synteny_matrices_local",train_synteny_matrices_local) +np.save("processed/train_branch_length_species",train_branch_length_species) +np.save("processed/train_branch_length_homology_species",train_branch_length_homology_species) +np.save("processed/train_mean_gene_length",train_mean_gene_length) +np.save("processed/train_dist_p_s",train_dist_p_s) +np.save("processed/train_dist_p_hs",train_dist_p_hs) +np.save("processed/train_distance",train_distance) +np.save("processed/train_labels",train_labels) + +print("Data Saved Successfully to processed :)") diff --git a/prepare_train_data.py b/prepare_train_data.py new file mode 100644 index 0000000..953522d --- /dev/null +++ b/prepare_train_data.py @@ -0,0 +1,98 @@ +import pandas as pd +import numpy as np +def create_branch_length_padding(bl): + maxlen=0 + for x in bl: + if len(x)>maxlen: + maxlen=len(x) + + for x in bl: + for i in range(len(x),maxlen): + x.append(0) + +def train_data(indexes,synteny_matrices_global,synteny_matrices_local,df,branch_length_species,branch_length_homology_species,distance,dist_p_s,dist_p_hs,gene_sequences): + """ + homology_type_counts=dict(df.homology_type.value_counts()) + homology_species_counts=dict(df.homology_species.value_counts()) + + #limits for each species and homology type in the trainig data so that the dataset is balanced + max_species_count=10000 + max_homology_type_count=25000 + for species in homology_species_counts: + homology_species_counts[species]=0 + + for homology_type in homology_type_counts: + homology_type_counts[homology_type]=0 + """ + labels=dict(ortholog_one2one=0,other_paralog=1,ortholog_one2many=2,ortholog_many2many=3,within_species_paralog=4) + + train_data_dataframe=pd.DataFrame() + train_labels=[] + train_indexes=[] + train_mean_gene_length=[] + for i in range(len(indexes)): + row=df.loc[indexes[i]] + if row["homology_type"]=="within_species_paralog": + continue + train_data_dataframe=train_data_dataframe.append(row) + train_indexes.append(i) + train_labels.append(labels[row["homology_type"]]) + train_mean_gene_length.append((len(gene_sequences[row["gene_stable_id"]])+len(gene_sequences[row["homology_gene_stable_id"]]))/2) + + train_synteny_matrices_global=synteny_matrices_global[train_indexes] + train_synteny_matrices_local=synteny_matrices_local[train_indexes] + train_distance=distance[train_indexes] + train_dist_p_s=dist_p_s[train_indexes] + train_dist_p_hs=dist_p_hs[train_indexes] + train_branch_length_species=branch_length_species[train_indexes] + train_branch_length_homology_species=branch_length_homology_species[train_indexes] + + create_branch_length_padding(train_branch_length_species) + train_branch_length_species=np.array(train_branch_length_species) + create_branch_length_padding(train_branch_length_homology_species) + train_branch_length_homology_species=np.array(train_branch_length_homology_species) + + #renormalize the train_mean_gene_length by (x-mean)/std + train_mean_gene_length=(train_mean_gene_length-np.mean(train_mean_gene_length))/np.std(train_mean_gene_length) + + #create a random array of permutations to shuffle the indices + shi=np.random.permutation(len(train_labels)) + + train_branch_length_species=train_branch_length_species[shi] + print(train_branch_length_species.shape) + + train_branch_length_homology_species=train_branch_length_homology_species[shi] + print(train_branch_length_homology_species.shape) + + train_dist_p_s=np.array(train_dist_p_s) + train_dist_p_s=train_dist_p_s[shi] + print(train_dist_p_s.shape) + + train_dist_p_hs=np.array(train_dist_p_hs) + train_dist_p_hs=train_dist_p_hs[shi] + print(train_dist_p_hs.shape) + + train_synteny_matrices_global=train_synteny_matrices_global[shi] + print(train_synteny_matrices_global.shape) + + train_synteny_matrices_local=train_synteny_matrices_local[shi] + print(train_synteny_matrices_local.shape) + + train_indexes=np.array(train_indexes) + train_indexes=train_indexes[shi] + print(train_indexes.shape) + + train_labels=np.array(train_labels) + train_labels=train_labels[shi] + print(train_labels.shape) + + train_mean_gene_length=np.array(train_mean_gene_length) + train_mean_gene_length=train_mean_gene_length[shi] + print(train_mean_gene_length.shape) + + train_distance=np.array(train_distance) + train_distance=train_distance[shi] + train_distance=(train_distance-np.mean(train_distance))/np.std(train_distance) + train_distance.shape + + return train_synteny_matrices_global,train_synteny_matrices_local,train_branch_length_species,train_branch_length_homology_species,train_mean_gene_length,train_dist_p_s,train_dist_p_hs,train_distance,train_labels diff --git a/process_data.py b/process_data.py index a9dd9b0..33f373e 100644 --- a/process_data.py +++ b/process_data.py @@ -1,42 +1,74 @@ import pandas import gc import numpy as np +from save_data import save_data_json +from save_data import write_dict_json + +#process_data.py + + +def create_map_list(l): #this function maps the indexes to values + t={} + for i in range(len(l)): + t[l[i]]=i + + return t def get_nearest_neighbors(g,gs,n,a,d,ld,ldg): - ne=[] - nr=[] - gi=d[gs.capitalize()] - sldf=a[gi] - sld=ld[gi] - sldg=ldg[gi] - if g not in sldg: - return ne - i=sldg[g] + ne=[] #list to store the backward genes + nr=[] #list to store the forward genes + gi=d[gs.capitalize()] #get the address of the corresponding species to which the gene belongs whose neighbor has to be found + sldf=a[gi]#select the dataframe + try: + sld=ld[gi]#see if the corresponding gene map exists + except: + #print("Length of Dataframes:{} \t Length of Loaded Genes:{} \t Length of Loaded Genomes Dictionaries:{}".format(len(a),len(ld),len(ldg))) + return ne,nr + sldg=ldg[gi]#select the corresponding map + if g not in sldg:#if the gene is not present in the dataframe return empty lists + return ne,nr + i=sldg[g]#find the index of the gnes #get the -n neighbors - start=int(sldf.iloc[i,[3]]) + start=int(sldf.iloc[i,[3]])#get the start location of the gene + flag=0 for j in range(n): - #select the column + if flag==1: + ne.append("NULL_GENE") + continue itemp=0 + #select the column end=sldf.iloc[:,4] end=np.array(end) assert(len(end)==len(sld)) - end=end-start - end_s=np.argsort(end) - for k in end_s: - if end[k]<0 and end[k+1]>=0: + end=end-start #subtract start from it so as to get relative position + end_s=np.argsort(end)#sort them by the order of distance + if end[end_s[0]]>=0:#if all the genes end ahead of the one in considertion + flag=1#increment the pointer + ne.append("NULL_GENE")#append the NULL_GENE value + continue + for k in end_s:#iterate through the sorted array + if end[k]<0 and end[k+1]>=0:#find the first value that is negative and the next one is positive to get the nearest gene itemp=k break - ne.append(sld[itemp]) - start=int(sldf.iloc[itemp,[3]]) + ne.append(sld[itemp])#push the gene in the array + start=int(sldf.iloc[itemp,[3]])#make "start" the start location of the current gene #print(start) #get the +n neighbors + flag=0 end=int(sldf.iloc[i,[4]]) for j in range(n): + if flag==1: + nr.append("NULL_GENE") + continue itemp=0 start=sldf.iloc[:,3] start=np.array(start) start=start-end start_s=np.argsort(start) + if start[start_s[-1]]<0: + flag=1 + nr.append("NULL_GENE") + continue for k in start_s: if start[k]>0: itemp=k @@ -46,38 +78,62 @@ def get_nearest_neighbors(g,gs,n,a,d,ld,ldg): return ne,nr -ls=[] -ld=[] def list_dict_genomes(a,n): + lst=[] + ldt=[] for x in a: - ldg={} + ldgt={} uc=list(x["gene_id"]) for i in range(len(uc)): - ldg[uc[i]]=i - ls.append(uc) - ld.append(ldg) - return ls,ld + ldgt[uc[i]]=i + lst.append(uc) + ldt.append(ldgt) + return lst,ldt -def create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg): + +def create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,save_after,enable_break): lsy={} #dictionary which stores +/- n genes of the given gene by id. Each key is a gene id which corresponds to the one in center. - lcmap={} #dictionary which stores the gene pairs already considered t=0 + c=0 + lsytemp={} + name="neighbor_genes" for df in a_h: for index,row in df.iterrows(): x=row["gene_stable_id"] y=row["homology_gene_stable_id"] xs=row["species"] ys=row["homology_species"] - if x+y in lcmap or y+x in lcmap: - continue - if x not in lsy: - xl,xr=get_nearest_neighbors(x,xs,n,a,d,ld,ldg) - lsy[x]=dict(b=xl,f=xr) - if y not in lsy: - yarr=[] - yl,yr=get_nearest_neighbors(y,ys,n,a,d,ld,ldg) - lsy[y]=dict(b=yl,f=yr) - lcmap[x+y]=1 - - - return lsy,lcmap + try: + z=lsy[x] + except: + try: + t2=d[xs.capitalize()]#see if the species exist in genomic maps + xl,xr=get_nearest_neighbors(x,xs,n,a,d,ld,ldg) + if len(xl)!=0:#check if neighboring genes were successfully found + lsy[x]=dict(b=xl,f=xr) + lsytemp[x]=dict(b=xl,f=xr) + except: + continue + try: + z=lsy[y] + except: + try: + t2=d[ys.capitalize()] + yl,yr=get_nearest_neighbors(y,ys,n,a,d,ld,ldg) + if len(yl)!=0: + lsy[y]=dict(b=yl,f=yr) + lsytemp[y]=dict(b=yl,f=yr) + except: + continue + t+=1 + if t>=save_after: + t=0 + c+=1 + write_dict_json(name+str(c),"processed",lsytemp) + lsytemp={} + if enable_break==1: + break + c+=1 + write_dict_json(name+str(c),"processed",lsytemp) + write_dict_json(name,"processed",lsy) + return lsy diff --git a/read_data.py b/read_data.py index 2211bc8..ba43f96 100644 --- a/read_data.py +++ b/read_data.py @@ -18,9 +18,9 @@ def read_data_genome(dir_name,a,dict_ind_genome): if len(lf)==0: print("No files in the data directory!!!!!!") sys.exit(1) - for x in lf: - - data_gene=pd.read_csv(dir_name+"/"+x,compression='gzip',sep='\t',comment='#',header=None) + print("Going to read data:") + for x in range(len(lf)): + data_gene=pd.read_csv(dir_name+"/"+lf[x],compression='gzip',sep='\t',comment='#',header=None) #print(data_gene.head) data_gene=data_gene[data_gene[2]=="gene"] data_gene=data_gene.sort_values(3) @@ -29,16 +29,21 @@ def read_data_genome(dir_name,a,dict_ind_genome): data_gene[["gene_id","gene_version","gene_name","gene_source","gene_biotype"]]=tmp data_gene=data_gene.drop(8,axis=1) #print(data_gene[0:10]) - for y in ["gene_version","gene_name","gene_source","gene_biotype","gene_id"]: - data_gene[y]=data_gene[y].apply(clear_data) - #print(data_gene[0:10]) + try: + for y in ["gene_version","gene_name","gene_source","gene_biotype","gene_id"]: + data_gene[y]=data_gene[y].apply(clear_data) + except: + continue + #print(data_gene[0:10]) data_gene=data_gene[data_gene['gene_biotype']=='protein_coding'] a.append(data_gene) - n=x.split(".")[0] + n=lf[x].split(".")[0] dict_ind_genome[n]=len(a)-1 return a,dict_ind_genome -def read_data_homology(dir,a_h,d_h): +def read_data_homology(dir): + a_h=[] + d_h={} lf=os.listdir(dir) if len(lf)==0: print("No Files in the Directory!!!!!!!") diff --git a/read_get_gene_seq.py b/read_get_gene_seq.py new file mode 100644 index 0000000..c1221e8 --- /dev/null +++ b/read_get_gene_seq.py @@ -0,0 +1,106 @@ +import json +from Bio import SeqIO +import pandas as pd +import os +import gzip + +def read_from_multiple_lsy(lsyfl): + lsy={} + for f in lsyfl: + d={} + with open(f,"r") as file: + d=dict(json.load(file)) + for t in d: + lsy[t]=d[t] + return lsy + +#this function updates the given dictionary with the given keys and values list +def create_dict(keys,values,dictionary): + for i in range(len(keys)): + if keys[i] not in dictionary: + dictionary[keys[i]]=values[i] + + return dictionary + +#this function maps all the genes to their respective species. +#(Function: when finding the species of any gene we do not need to search the entire dataframe) +def group_seq_by_species(df): + g_to_sp={} + sph=list(df.homology_species) + ghsp=list(df.homology_gene_stable_id) + sp=list(df.species) + gsp=list(df.gene_stable_id) + create_dict(gsp,sp,g_to_sp) + create_dict(ghsp,sph,g_to_sp) + return g_to_sp + +#this function returns the gene-id and gene-biotype from the description in the fasta file record. +def description_cleaner(description): + description=description.split() + t="" + gbt="" + for x in description: + try: + x=x.split(":") + if x[0]=="gene": + t=x[1].split(".")[0] + if x[0]=="gene_biotype": + gbt=x[1] + except: + return "aa","aa" + return t,gbt + +def read_gene_seq(dirname,s,genes_by_species): + lof=os.listdir(dirname)#list all the files in the sequences directory + ftr=[] + for f in lof: + if f.split(".")[0] in s:#check whether the species is present in the species to read list. Will skip those species which are not present in the dataframe + ftr.append(f) + data={} + for f in ftr: + species=f.split(".")[0].lower() + with gzip.open(dirname+"/"+f,"rt") as file: + record=SeqIO.parse(file,"fasta") + for r in record: + gid,gbt=description_cleaner(r.description) + if str(gid) not in data and str(gid) in genes_by_species[species] and gbt=="protein_coding": + data[gid]=str(r.seq) + return data + +def read_gene_sequences(df,lsy,data_dir,fname): + + """The basic idea here is to create a list/dictionary of all the genes by their species. + Once the mapping is done, all the respective fasta sequence files are read by Species + and the CDNA sequences for each gene in the species record are read and stored. + Thus we don't have to read the same file multiple times.""" + + grouped_genes=group_seq_by_species(df) + gene_by_species_dict={} + for i in df.homology_species.unique(): + gene_by_species_dict[i]=[] + for x in lsy: + species=grouped_genes[x]#get the species + if x not in gene_by_species_dict[species]:#check if the gene already exists in the species dict or not. + gene_by_species_dict[species].append(x) + xl=lsy[x]['b'] + xr=lsy[x]['f'] + for gxl in xl: + if gxl=="NULL_GENE": + break + if gxl not in gene_by_species_dict[species]: + gene_by_species_dict[species].append(gxl) + for gxr in xr: + if gxr=="NULL_GENE": + break + if gxr not in gene_by_species_dict[species]: + gene_by_species_dict[species].append(gxr) + + s=[x for x in gene_by_species_dict if len(gene_by_species_dict[x])!=0]#select those species only whose gene sequences we have to read. + s=[x.capitalize() for x in s] + + data=read_gene_seq(data_dir,s,gene_by_species_dict) + + with open("processed/"+fname+".json","w") as file:#save the data + json.dump(data,file) + + return data diff --git a/req_data.py b/req_data.py index 82016c5..a58b4cb 100644 --- a/req_data.py +++ b/req_data.py @@ -3,6 +3,11 @@ import pandas import sys import urllib.request as urllib import pandas as pd +import json +import requests +import time + +from process_data import create_map_list lf=[] @@ -24,3 +29,4 @@ def get_data_file(file,dir): os.mkdir(dir) for x in lf: download_data(x,dir) + diff --git a/save_data.py b/save_data.py new file mode 100644 index 0000000..dad130f --- /dev/null +++ b/save_data.py @@ -0,0 +1,27 @@ +import os +import pickle +import json +import sys + +def write_dict_json(name,dir,d): + if not os.path.exists(dir): + os.mkdir(dir) + + path=os.path.join(dir,name+".json") + with open(path,'w') as file: + json.dump(d,file) + +def write_file_multiple_json(name,dir,l): + if not os.path.exists(dir): + os.mkdir(dir) + + path=os.path.join(dir,name) + with open(path,"w") as file: + for x in l: + out=json.dumps(x) + file.write(out) + file.write('\n') + +def save_data_json(name,data): + with open(name+".json","w")as file: + json.dump(data,file) diff --git a/species_tree.tree b/species_tree.tree new file mode 100644 index 0000000..17e7df8 --- /dev/null +++ b/species_tree.tree @@ -0,0 +1 @@ +(saccharomyces_cerevisiae:0.12053,(((ciona_intestinalis:0.151113,ciona_savignyi:0.154677)1:0.0327515,((eptatretus_burgeri:0.128378,petromyzon_marinus:0.169182)1:0.00095149,((lepisosteus_oculatus:0.156711,((scleropages_formosus:0.116055,paramormyrops_kingsleyae:0.119385)1:0.0152653,((((astyanax_mexicanus:0.0831834,pygocentrus_nattereri:0.0898266)1:0.0199602,ictalurus_punctatus:0.10061)1:0.00963079,danio_rerio:0.110582)1:0.00197338,(esox_lucius:0.130765,(gadus_morhua:0.117116,(periophthalmus_magnuspinnatus:0.114713,(hippocampus_comes:0.114444,(((gasterosteus_aculeatus:0.103752,labrus_bergylta:0.0866677)1:0.00446595,((tetraodon_nigroviridis:0.0674859,takifugu_rubripes:0.0744041)1:0.0351015,mola_mola:0.0978235)1:0.00229189)1:0.00589275,((((cynoglossus_semilaevis:0.0954686,scophthalmus_maximus:0.100441)1:0.0001,(seriola_lalandi_dorsalis:0.0127928,seriola_dumerili:0.0127672)1:0.0693065)1:0.0001,((mastacembelus_armatus:0.0748029,monopterus_albus:0.0860571)1:0.00345198,anabas_testudineus:0.078993)1:0.0110165)1:0.0121107,(((((amphiprion_ocellaris:0.00675482,amphiprion_percula:0.00552518)1:0.0268836,acanthochromis_polyacanthus:0.0400714)1:0.0212642,stegastes_partitus:0.0517763)1:0.0269893,(amphilophus_citrinellus:0.0553949,(oreochromis_niloticus:0.0230128,(neolamprologus_brichardi:0.0155877,(((astatotilapia_calliptera:0.00237541,maylandia_zebra:0.00222459)1:0.000515714,pundamilia_nyererei:0.00887929)1:0.000673283,haplochromis_burtoni:0.00746981)1:0.00820467)1:0.0105078)1:0.0381336)1:0.0211695)1:0.0181227,((((oryzias_latipes_hni:0.0165617,oryzias_latipes:0.0147983)1:0.00102726,oryzias_latipes_hsok:0.0171027)1:0.0446085,oryzias_melastigma:0.0589057)1:0.0512701,(kryptolebias_marmoratus:0.0930933,((fundulus_heteroclitus:0.0803113,cyprinodon_variegatus:0.0838287)1:0.0001,(((xiphophorus_couchianus:0.0110953,xiphophorus_maculatus:0.00742472)1:0.0241526,gambusia_affinis:0.0308374)1:0.0098865,(poecilia_reticulata:0.0286208,(poecilia_mexicana:0.00793647,(poecilia_formosa:0.00521391,poecilia_latipinna:0.00440609)1:0.00198853)1:0.0163947)1:0.010153)1:0.043138)1:0.0137214)1:0.0121596)1:0.0001)1:0.0001)1:0.00168516)1:0.0151377)1:0.0017499)1:0.00383953)1:0.00032757)1:0.0116671)1:0.0001)1:0.00645315)1:0.0001,(latimeria_chalumnae:0.119274,(xenopus_tropicalis:0.14571,((ornithorhynchus_anatinus:0.127358,((((phascolarctos_cinereus:0.0589639,notamacropus_eugenii:0.0630161)1:0.0126856,sarcophilus_harrisii:0.0666494)1:0.00245254,monodelphis_domestica:0.0752439)1:0.0417722,((((loxodonta_africana:0.0624292,procavia_capensis:0.0761308)1:0.0167753,echinops_telfairi:0.0966547)1:0.019774,(dasypus_novemcinctus:0.0776558,choloepus_hoffmanni:0.0782942)1:0.0226792)1:0.00493446,(((erinaceus_europaeus:0.101031,sorex_araneus:0.107679)1:0.00110986,((((panthera_pardus:0.00386325,panthera_tigris_altaica:0.00452675)1:0.00752221,felis_catus:0.0116228)1:0.0534858,((((ursus_americanus:0.00367107,ursus_maritimus:0.00321893)1:0.0138769,ailuropoda_melanoleuca:0.0173581)1:0.0376917,(neovison_vison:0.0160509,mustela_putorius_furo:0.0160591)1:0.0386994)1:0.00641809,((canis_lupus_dingo:0.000945102,canis_familiaris:0.0013649)1:0.0104717,vulpes_vulpes:0.0119883)1:0.0472989)1:0.00602788)1:0.0163724,((pteropus_vampyrus:0.0800095,myotis_lucifugus:0.0826505)1:0.00681959,((((((capra_hircus:0.0112773,ovis_aries:0.0110627)1:0.0167979,((bos_mutus:0.00522287,bos_taurus:0.00464713)1:0.0001,bison_bison_bison:0.00477814)1:0.0226132)1:0.0437826,tursiops_truncatus:0.066084)1:0.0114156,sus_scrofa:0.075862)1:0.0001,vicugna_pacos:0.0815824)1:0.00472105,(equus_caballus:0.00639876,equus_asinus_asinus:0.00664124)1:0.0763951)1:0.00509347)1:0.0001)1:0.0158911)1:0.0001,((((((prolemur_simus:0.0371602,propithecus_coquereli:0.0381898)1:0.00515634,microcebus_murinus:0.0414087)1:0.0276342,otolemur_garnettii:0.0746654)1:0.00700933,(carlito_syrichta:0.0797062,((((saimiri_boliviensis_boliviensis:0.0262284,cebus_capucinus:0.0261016)1:0.00267284,callithrix_jacchus:0.0301072)1:0.0001,aotus_nancymaae:0.0251968)1:0.0236875,((((((pan_paniscus:0.00332066,pan_troglodytes:0.00221934)1:0.0043011,homo_sapiens:0.0066089)1:0.00183326,gorilla_gorilla:0.00859971)1:0.00841205,pongo_abelii:0.0171753)1:0.00280069,nomascus_leucogenys:0.0195845)1:0.0110611,(((piliocolobus_tephrosceles:0.00971396,colobus_angolensis_palliatus:0.011096)1:0.0016942,(rhinopithecus_roxellana:0.00209338,rhinopithecus_bieti:0.00301662)1:0.0096733)1:0.00485545,(chlorocebus_sabaeus:0.0116617,(((mandrillus_leucophaeus:0.00675832,cercocebus_atys:0.00569168)1:0.000532647,(theropithecus_gelada:0.00408503,papio_anubis:0.00409497)1:0.00243485)1:0.00115762,((macaca_mulatta:0.00237678,macaca_fascicularis:0.00210322)1:0.000897997,macaca_nemestrina:0.004232)1:0.0043418)1:0.00406851)1:0.00498204)1:0.013507)1:0.0178328)1:0.0275728)1:6.4199e-05)1:0.0159811,tupaia_belangeri:0.0928607)1:0.0060941,((oryctolagus_cuniculus:0.0779437,ochotona_princeps:0.0908063)1:0.0251776,(((((((cavia_aperea:0.0125214,cavia_porcellus:0.0081786)1:0.0662348,chinchilla_lanigera:0.0704902)1:0.00415561,octodon_degus:0.0759278)1:0.00136511,((heterocephalus_glaber_female:0.000514822,heterocephalus_glaber_male:0.000755178)1:0.0562411,fukomys_damarensis:0.0576989)1:0.0193318)1:0.0206755,(((urocitellus_parryii:0.012115,spermophilus_dauricus:0.013075)1:0.000664133,ictidomys_tridecemlineatus:0.0138159)1:0.00572007,marmota_marmota_marmota:0.0183052)1:0.0711537)1:0.00255843,(dipodomys_ordii:0.0916093,castor_canadensis:0.0925507)1:0.00517574)1:0.00430836,(jaculus_jaculus:0.0947278,(nannospalax_galili:0.0880758,((((((cricetulus_griseus_chok1gshd:0.000317157,cricetulus_griseus_crigri:0.000692843)1:0.00075912,cricetulus_griseus_picr:0.00068588)1:0.0478454,mesocricetus_auratus:0.0525732)1:0.016122,peromyscus_maniculatus_bairdii:0.0636996)1:0.0034916,microtus_ochrogaster:0.0703622)1:0.010456,(meriones_unguiculatus:0.0762326,(rattus_norvegicus:0.0630119,(mus_pahari:0.0365573,(mus_caroli:0.0204404,((mus_spretus:0.00897036,mus_spicilegus:0.00853964)1:0.00109415,(mus_musculus_casteij:0.00513862,(mus_musculus_pwkphj:0.00538657,(mus_musculus_wsbeij:0.00195001,((((mus_musculus_c57bl6nj:0.000735102,mus_musculus:0.0001)1:0.000893176,mus_musculus_nzohlltj:0.00151682)1:0.0001,(((((mus_musculus_c3hhej:0.000389873,mus_musculus_cbaj:0.000530127)1:0.000227207,mus_musculus_dba2j:0.00115279)1:0.000313341,(mus_musculus_aj:0.000333528,mus_musculus_balbcj:0.000946472)1:6.09235e-05)1:0.000331345,mus_musculus_akrj:0.000692839)1:0.0001,(mus_musculus_fvbnj:0.000902944,mus_musculus_nodshiltj:0.00149706)1:0.000237201)1:0.000361644)1:0.0001,(mus_musculus_lpj:0.000605635,mus_musculus_129s1svimj:0.000464365)1:0.000977946)1:0.000417414)1:0.00280773)1:0.00039707)1:0.00550964)1:0.0108989)1:0.0168799)1:0.0230902)1:0.0163714)1:0.000606242)1:0.0121791)1:0.0098349)1:0.00734652)1:0.00362496)1:0.0001)1:0.0001)1:0.00193823)1:0.0127074)1:0.00755806)1:0.00172511,(((((pogona_vitticeps:0.104873,anolis_carolinensis:0.105197)1:0.0134109,notechis_scutatus:0.111449)1:0.0001,salvator_merianae:0.115498)1:0.0123771,sphenodon_punctatus:0.127209)1:0.00101294,((((chelonoidis_abingdonii:0.0226001,gopherus_agassizii:0.0228099)1:0.0125923,chrysemys_picta_bellii:0.0340277)1:0.0413409,pelodiscus_sinensis:0.078325)1:0.0404938,(crocodylus_porosus:0.112492,(((((apteryx_owenii:0.00132365,apteryx_haastii:0.00212635)1:0.00320653,apteryx_rowi:0.00451347)1:0.0332048,dromaius_novaehollandiae:0.0390731)1:0.0306607,nothoprocta_perdicaria:0.0673109)1:0.016828,(((((gallus_gallus:0.0370661,meleagris_gallopavo:0.0429639)1:0.00643237,coturnix_japonica:0.0436126)1:0.00137178,numida_meleagris:0.0488591)1:0.0296977,(anser_brachyrhynchus:0.031382,anas_platyrhynchos_platyrhynchos:0.028548)1:0.0385031)1:0.00748951,(((manacus_vitellinus:0.013022,lepidothrix_coronata:0.011958)1:0.0535377,(((cyanistes_caeruleus:0.0148277,parus_major:0.0153823)1:0.0308174,ficedula_albicollis:0.0455876)1:0.000531936,(((taeniopygia_guttata:0.0208879,lonchura_striata_domestica:0.0200021)1:0.0209237,serinus_canaria:0.0413113)1:0.0001,(zonotrichia_albicollis:0.0116596,junco_hyemalis:0.0129704)1:0.0283468)1:0.00859797)1:0.0174028)1:0.0105652,((calidris_pugnax:0.0169552,calidris_pygmaea:0.0181848)1:0.0470247,melopsittacus_undulatus:0.0749353)1:0.00851032)1:0.0112749)1:0.00147191)1:0.0214409)1:0.0020761)1:0.0172119)1:0.00693331)1:0.0127467)1:0.0001)1:0.0180963)1:0.00379002)1:0.0416798)1:0.00981196,(drosophila_melanogaster:0.155163,caenorhabditis_elegans:0.188157)1:0.0001)1:0.12053); \ No newline at end of file diff --git a/test_main.py b/test_main.py new file mode 100644 index 0000000..dddb487 --- /dev/null +++ b/test_main.py @@ -0,0 +1,54 @@ +import sys +import numpy as np +import pandas as pd +import json +from get_data import get_data_homology,get_data_genome +from process_data import create_data_homology_ls +from read_get_gene_seq import read_gene_sequences +from create_synteny_matrix import synteny_matrix +from tree_data import create_tree_data +from create_train_data import train_data +from train import train + +arg=sys.argv +arg=arg[1:] + +if len(arg)!=5: + print("No. of arguments more or less. Please check") + sys.exit(1) + +dir_g="data" +ld,ldg,a,d=get_data_genome(arg,dir_g) + +dir_hom="data_homology" +a_h,d_h=get_data_homology(arg,dir_hom) + +if arg[-1]=="-d": + sys.exit(1) + +print("Data Read") + +n=2 #no. of numbers neighbors +save_after=3 #to save data after n steps +lsy=create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,save_after) +print(len(lsy)) + +print("Neighbor Genes Found") + +gene_sequences=read_gene_sequences(a_h[0],lsy,"geneseq","gene_sequences") + +synteny_matrices,indexes=synteny_matrix(gene_sequences,a_h[0][0:1000],lsy,n) +print("Synteny Matrices are created successfully\n",len(indexes),"\n",len(synteny_matrices)) +np.save("synteny_matrices",synteny_matrices) +np.save("indexes",indexes) + +with open("gene_seq_updated.json","w") as file: + json.dump(gene_sequences,file) + +df=a_h[0].loc[indexes] +branch_length_species,branch_length_homology_species,distance,dist_p_s,dist_p_hs=create_tree_data("species_tree.tree",df) + +train_synteny_matrices,train_branch_length_species,train_branch_length_homology_species,train_mean_gene_length,train_dist_p_s,train_dist_p_hs,train_distance,train_labels=train_data(indexes,synteny_matrices,df,branch_length_species,branch_length_homology_species,distance,dist_p_s,dist_p_hs,gene_sequences) + +train(train_synteny_matrices,train_branch_length_species,train_branch_length_homology_species,train_mean_gene_length, + train_dist_p_s,train_dist_p_hs,train_distance,train_labels) diff --git a/train.py b/train.py new file mode 100644 index 0000000..0621ce5 --- /dev/null +++ b/train.py @@ -0,0 +1,67 @@ +import tensorflow as tf +import numpy as np + +from model import create_model + +def train(train_synteny_matrices,train_branch_length_species,train_branch_length_homology_species,train_mean_gene_length,train_dist_p_s,train_dist_p_hs,train_distance,train_labels): + graph,saver=create_model() + synm,bls,blhs,gl,dps,dphs,dis,lr,y=graph.get_collection("input_nodes") + loss,t_op,accuracy,init=graph.get_collection("output_nodes") + with tf.Session(graph=graph) as sess: + sess.run(init) + batch_size=64 + num_epochs=30 + learn=0.001 + for j in range(num_epochs): + for i in range(50000//batch_size): + feed_dict={ + synm:train_synteny_matrices[i*batch_size:(i+1)*batch_size], + 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], + gl:train_mean_gene_length[i*batch_size:(i+1)*batch_size].reshape((batch_size,1)), + 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) + + feed_dict={ + synm:train_synteny_matrices[i*batch_size:(i+1)*batch_size].transpose((0,2,1,3)), + blhs:train_branch_length_species[i*batch_size:(i+1)*batch_size], + bls:train_branch_length_homology_species[i*batch_size:(i+1)*batch_size], + gl:train_mean_gene_length[i*batch_size:(i+1)*batch_size].reshape((batch_size,1)), + dphs:train_dist_p_s[i*batch_size:(i+1)*batch_size].reshape((batch_size,1)), + dps: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={ + synm:train_synteny_matrices[-15000:], + bls:train_branch_length_species[-15000:], + blhs:train_branch_length_homology_species[-15000:], + gl:train_mean_gene_length[-15000:].reshape((15000,1)), + dps:train_dist_p_s[-15000:].reshape((15000,1)), + dphs:train_dist_p_hs[-15000:].reshape((15000,1)), + dis:train_distance[-15000:].reshape((15000,1)), + y:train_labels[-15000:] + } + train_dict={ + synm:train_synteny_matrices[:15000], + bls:train_branch_length_species[:15000], + blhs:train_branch_length_homology_species[:15000], + gl:train_mean_gene_length[:15000].reshape((15000,1)), + dps:train_dist_p_s[:15000].reshape((15000,1)), + dphs:train_dist_p_hs[:15000].reshape((15000,1)), + dis:train_distance[:15000].reshape((15000,1)), + y:train_labels[:15000] + } + accuracy_train,loss_train=sess.run([accuracy,loss],feed_dict=train_dict) + accuracy_test,loss_test=sess.run([accuracy,loss],feed_dict=test_dict) + print("Epoch:{} Train Accuracy:{} Train Loss:{} Test Accuracy:{} Test Loss:{}".format(j+1,accuracy_train*100,loss_train,accuracy_test*100,loss_test)) + learn*=0.97 + saver.save(sess,"saved_models/model.ckpt") diff --git a/tree_data.py b/tree_data.py new file mode 100644 index 0000000..8a19d55 --- /dev/null +++ b/tree_data.py @@ -0,0 +1,37 @@ +from ete3 import Tree +import numpy as np + +def create_tree_data(treename,df): + t=Tree(treename) + branch_lengths_s=[] + branch_lengths_hs=[] + dist=[] + ns=[] + nhs=[] + for index,row in df.iterrows(): + d=0 + x=row["species"] + y=row["homology_species"] + bl=[] + c=0 + mca=t.get_common_ancestor(x,y) + node=t&x + while node.up!=mca: + d+=node.dist + bl.append(node.dist) + node=node.up + c+=1 + ns.append(c) + c=0 + branch_lengths_s.append(bl) + bl=[] + node=t&y + while node.up!=mca: + d+=node.dist + bl.append(node.dist) + node=node.up + c+=1 + nhs.append(c) + branch_lengths_hs.append(bl) + dist.append(d) + return np.array(branch_lengths_s),np.array(branch_lengths_hs),np.array(dist),np.array(ns),np.array(nhs)