Added New Features

New Data Feature and Command Line Variable For Testing

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)

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

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.................")

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

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)
+if arg[-1]=="-d":
-#print(a_h[0][0:10],"\n",d_h)
+    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("Neighbor Genes Found")
-"""
-print(lsy,"\n",lcmap)
-lt=ldg[0]
-for x in lsy:
 
-    xr=lsy[x]['f']
+gene_sequences=read_gene_sequences(a_h[0],lsy,"geneseq","gene_sequences")
-    xl=lsy[x]['b']
+
-    for g in range(len(xl)-1,-1,-1):
+synteny_matrices,indexes=synteny_matrix(gene_sequences,a_h[0],lsy,n)
-        print(a[0].iloc[lt[xl[g]],[3,4]])
+print("Synteny Matrices are created successfully\n",len(indexes),"\n",len(synteny_matrices))
-    print("------------------\n",a[0].iloc[lt[x],[3,4]],"\n-------------------")
+np.save("synteny_matrices",synteny_matrices)
-    for g in xr:
+np.save("indexes",indexes)
-        print(a[0].iloc[lt[g],[3,4]])"""
+
+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)

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

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 :)")

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

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=[]
+    ne=[] #list to store the backward genes
-    nr=[]
+    nr=[] #list to store the forward genes
-    gi=d[gs.capitalize()]
+    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]
+    sldf=a[gi]#select the dataframe
-    sld=ld[gi]
+    try:
-    sldg=ldg[gi]
+        sld=ld[gi]#see if the corresponding gene map exists
-    if g not in sldg:
+    except:
-        return ne
+        #print("Length of Dataframes:{} \t Length of Loaded Genes:{} \t Length of Loaded Genomes Dictionaries:{}".format(len(a),len(ld),len(ldg)))
-    i=sldg[g]
+        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=end-start #subtract start from it so as to get relative position
-        end_s=np.argsort(end)
+        end_s=np.argsort(end)#sort them by the order of distance
-        for k in end_s:
+        if end[end_s[0]]>=0:#if all the genes end ahead of the one in considertion
-            if end[k]<0 and end[k+1]>=0:
+            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])
+        ne.append(sld[itemp])#push the gene in the array
-        start=int(sldf.iloc[itemp,[3]])
+        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
+            ldgt[uc[i]]=i
-        ls.append(uc)
+        lst.append(uc)
-        ld.append(ldg)
+        ldt.append(ldgt)
-    return ls,ld
+    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:
+            try:
-                continue
+                z=lsy[x]
-            if x not in lsy:
+            except:
-                xl,xr=get_nearest_neighbors(x,xs,n,a,d,ld,ldg)
+                try:
-                lsy[x]=dict(b=xl,f=xr)
+                    t2=d[xs.capitalize()]#see if the species exist in genomic maps
-            if y not in lsy:
+                    xl,xr=get_nearest_neighbors(x,xs,n,a,d,ld,ldg)
-                yarr=[]
+                    if len(xl)!=0:#check if neighboring genes were successfully found
-                yl,yr=get_nearest_neighbors(y,ys,n,a,d,ld,ldg)
+                        lsy[x]=dict(b=xl,f=xr)
-                lsy[y]=dict(b=yl,f=yr)
+                        lsytemp[x]=dict(b=xl,f=xr)
-            lcmap[x+y]=1
+                except:
-
+                    continue
-
+            try:
-    return lsy,lcmap
+                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

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:
+    print("Going to read data:")
-
+    for x in range(len(lf)):
-        data_gene=pd.read_csv(dir_name+"/"+x,compression='gzip',sep='\t',comment='#',header=None)
+        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"]:
+        try:
-            data_gene[y]=data_gene[y].apply(clear_data)
+            for y in ["gene_version","gene_name","gene_source","gene_biotype","gene_id"]:
-        #print(data_gene[0:10])
+                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!!!!!!!")

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

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)
+

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)

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)

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")

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)