Add files via upload

create_synteny_matrix.py

@@ -45,7 +45,7 @@ def create_synteny_matrix_mul(gene_seq,g1,g2,n):
         for j in range(n):
             if g2[j]=="NULL_GENE":
                 continue
-            norm_len=(len(gene_seq[g1[i]])+len(gene_seq[g2[j]]))
+            norm_len=max(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")
@@ -94,6 +94,6 @@ def synteny_matrix(gene_seq,hdf,lsy,n,enable_break):
         if t==5 and enable_break==1:
             break
     end=time.time()
-    print("Time Taken:",end-start)
-    print("Average Time:",(end-start)/len(sg))
+    #print("Time Taken:",end-start)
+    #print("Average Time:",(end-start)/len(sg))
     return np.array(sg),np.array(sl),np.array(ind)

open_synteny_matrices.py

@@ -0,0 +1,75 @@
+import numpy as np 
+import os
+import json
+from read_data import read_data_homology
+import matplotlib.pyplot as plt
+import seaborn as sns
+
+a_h,_=read_data_homology("data_homology")
+
+ndir="processed/synteny_matrices/"
+nf1="synteny_matrices_global_"
+nf2="synteny_matrices_local_"
+nf3="indexes_"
+
+lsy={}
+with open("processed/neighbor_genes.json","r") as file:
+    lsy=dict(json.load(file))
+
+synteny_matrices_global=np.load(ndir+nf1+"11"+".npy")
+synteny_matrices_local=np.load(ndir+nf2+"11"+".npy")
+indexes=np.load(ndir+nf3+"11"+".npy")
+
+df=a_h[0].loc[indexes]
+
+inddict={}
+for i in range(len(indexes)):
+    inddict[indexes[i]]=i
+
+while(1):
+    i=int(input("Enter the index"))
+    if i in inddict:
+        print("Species",df.loc[i].species)
+        print("Homology Species",df.loc[i].homology_species)
+        print("Gene Stable Id:",df.loc[i].gene_stable_id)
+        print("Homology Gene Stable Id:",df.loc[i].homology_gene_stable_id)
+        print("Global aligned matrix:")
+
+        g1=df.loc[i].gene_stable_id
+        g2=df.loc[i].homology_gene_stable_id
+        x=[]
+        y=[]
+        for n in range(len(lsy[g1]['b'])-1,-1,-1):
+            x.append(lsy[g1]['b'][n])
+        x.append(g1)
+        for k in lsy[g1]['f']:
+            x.append(k)
+
+        for n in range(len(lsy[g2]['b'])-1,-1,-1):
+            y.append(lsy[g2]['b'][n])
+        y.append(g2)
+        for k in lsy[g2]['f']:
+            y.append(k)
+
+        loc=inddict[i]
+        sg=synteny_matrices_global[loc]
+        sl=synteny_matrices_local[loc]
+        for m in range(sg.shape[-1]):
+            matrix=sg[:,:,m]
+            print(matrix)
+            hmap=sns.heatmap(matrix,xticklabels=y, yticklabels=x,annot=True)
+            plt.show()
+
+        print("Local Alignment Matrix")
+        for m in range(sg.shape[-1]):
+            matrix=sl[:,:,m]
+            print(matrix)
+            hmap=sns.heatmap(matrix,xticklabels=y, yticklabels=x,annot=True)
+            plt.show()
+
+            
+
+    else:
+        print("Index not found")
+
+

prepare_synteny_matrix.py

@@ -0,0 +1,57 @@
+import pandas as pd 
+import numpy as np 
+import os 
+import sys
+import progressbar
+import json
+from read_data import read_data_homology
+from read_get_gene_seq import read_gene_sequences
+from create_synteny_matrix import synteny_matrix
+
+if not os.path.exists("processed/synteny_matrcies"):
+    os.mkdir("processed/synteny_matrices")
+
+arg=sys.argv
+arg=arg[1:]
+
+enable_break=0
+
+if arg[-1]=="-test":
+    enable_break=1
+
+lsy={}
+a_h,d_h=read_data_homology("data_homology")
+print("Homology Data Read")
+with open("processed/neighbor_genes.json","r") as file:
+    lsy=dict(json.load(file))
+print(len(lsy))
+print("Neighbor Genes Loaded")
+if enable_break==1:
+    gene_sequences=read_gene_sequences(a_h,lsy,"geneseq","gene_sequences")
+else:
+    gene_sequences=read_gene_sequences(a_h[0],lsy,"geneseq","gene_sequences")
+print("Gene Sequences Loaded")
+
+
+if enable_break==1:
+    save_after=10
+else:
+    save_after=250000
+n=3
+c=0
+j=0
+ndir="processed/synteny_matrices/"
+nf1="synteny_matrices_global"
+nf2="synteny_matrices_local"
+nf3="indexes"
+for df in a_h:
+    j+=1
+    for i in progressbar.progressbar(range(0,len(df)//save_after)):
+        synteny_matrices_global,synteny_matrices_local,indexes=synteny_matrix(gene_sequences,df[i*save_after:(i+1)*save_after],lsy,n,enable_break)
+        np.save(ndir+nf1+"_"+str(j)+str(c+1),synteny_matrices_global)
+        np.save(ndir+nf2+"_"+str(j)+str(c+1),synteny_matrices_local)
+        np.save(ndir+nf3+"_"+str(j)+str(c+1),indexes)
+        c+=1
+        if enable_break==1:
+            break
+print("Synteny Matrices Created Successfully :)")

process_data.py

@@ -1,6 +1,8 @@
 import pandas
 import gc
 import numpy as np
+import json
+import progressbar
 from save_data import save_data_json
 from save_data import write_dict_json
 
@@ -14,6 +16,18 @@ def create_map_list(l): #this function maps the indexes to values
 
     return t
 
+def list_dict_genomes(a,n):
+    lst=[]
+    ldt=[]
+    for x in a:
+        ldgt={}
+        uc=list(x["gene_id"])
+        for i in range(len(uc)):
+            ldgt[uc[i]]=i
+        lst.append(uc)
+        ldt.append(ldgt)
+    return lst,ldt
+
 def get_nearest_neighbors(g,gs,n,a,d,ld,ldg):
     ne=[] #list to store the backward genes
     nr=[] #list to store the forward genes
@@ -27,7 +41,7 @@ def get_nearest_neighbors(g,gs,n,a,d,ld,ldg):
     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
+    i=sldg[g]#find the index of the genes
     #get the -n neighbors
     start=int(sldf.iloc[i,[3]])#get the start location of the gene
     flag=0
@@ -42,7 +56,7 @@ def get_nearest_neighbors(g,gs,n,a,d,ld,ldg):
         assert(len(end)==len(sld))
         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
+        if end[end_s[0]]>=0:#if all the genes end ahead of the one in consideration
             flag=1#increment the pointer
             ne.append("NULL_GENE")#append the NULL_GENE value
             continue
@@ -78,27 +92,18 @@ def get_nearest_neighbors(g,gs,n,a,d,ld,ldg):
 
     return ne,nr
 
-def list_dict_genomes(a,n):
-    lst=[]
-    ldt=[]
-    for x in a:
-        ldgt={}
-        uc=list(x["gene_id"])
-        for i in range(len(uc)):
-            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,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.
     t=0
-    c=0
+    with open("processed/neighbor_genes.json","r") as file:
+        lsy=dict(json.load(file))
+    print("Existing neighbor genes read!!")
+    print(len(lsy))
+    c=4
     lsytemp={}
     name="neighbor_genes"
     for df in a_h:
-        for index,row in df.iterrows():
+        for _,row in progressbar.progressbar(df.iterrows()):
             x=row["gene_stable_id"]
             y=row["homology_gene_stable_id"]
             xs=row["species"]

read_get_gene_seq.py

@@ -24,8 +24,7 @@ def create_dict(keys,values,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={}
+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)
@@ -67,19 +66,24 @@ def read_gene_seq(dirname,s,genes_by_species):
                     data[gid]=str(r.seq)
     return data
 
-def read_gene_sequences(df,lsy,data_dir,fname):
+def read_gene_sequences(hdf,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)
+    grouped_genes={}
     gene_by_species_dict={}
-    for i in df.homology_species.unique():
-        gene_by_species_dict[i]=[]
+    for df in hdf:
+        grouped_genes=group_seq_by_species(df,grouped_genes)        
+        for i in df.homology_species.unique():
+            if i not in gene_by_species_dict:
+                gene_by_species_dict[i]=[]
     for x in lsy:
-        species=grouped_genes[x]#get the species
+        try:
+            species=grouped_genes[x]#get the species
+        except:
+            continue
         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']
@@ -95,11 +99,10 @@ def read_gene_sequences(df,lsy,data_dir,fname):
             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]
     
-    print(len(s))
-
     data=read_gene_seq(data_dir,s,gene_by_species_dict)
 
     with open("processed/"+fname+".json","w") as file:#save the data