Add files via upload

access_data_rest.py

@@ -0,0 +1,49 @@
+import json
+import requests
+import progressbar
+from create_synteny_matrix import update
+
+def update_rest(data):
+    gids={}
+    with open("processed/not_found.json","r") as file:
+        gids=dict(json.load(file))
+
+    gids=list(gids.keys())
+
+    geneseq={}
+
+    server = "https://rest.ensembl.org"
+    ext = "/sequence/id?type=cds"
+    headers={ "Content-Type" : "application/json", "Accept" : "application/json"}
+
+    for i in progressbar.progressbar(range(0,len(gids)-50,50)):
+        ids=dict(ids=list(gids[i:i+50]))
+        while(1):
+            try:
+                r = requests.post(server+ext, headers=headers, data=str(json.dumps(ids)))
+                if not r.ok:
+                    r.raise_for_status()
+                gs=r.json()
+                tgs={}
+                for g in gs:
+                    tgs[g["query"]]=g["seq"]
+                geneseq.update(tgs)
+                break
+            except Exception as e:
+                print("Error:",e)
+                continue
+
+    data.update(geneseq)
+    for genes in gids:
+        try:
+            _=data[genes]
+        except:
+            print(genes)
+            update(data,genes)
+
+    with open("processed/gene_sequences.json","w") as file:
+        json.dump(data,file)
+
+    print("Gene Sequences Updated Successfully")
+    return data
+            

create_synteny_matrix.py

@@ -12,17 +12,31 @@ def update(gene_seq,gene):
     while(1):
         try:
             server = "https://rest.ensembl.org"
-            ext = "/sequence/id/"+str(gene)+"?"
+            ext = "/sequence/id/"+str(gene)+"?type=cds;multiple_sequences=1"
 
-            r = requests.get(server+ext, headers={ "Content-Type" : "text/plain"})
+            r = requests.get(server+ext, headers={ "Content-Type" : "application/json"})
 
             if not r.ok:
                 r.raise_for_status()
                 sys.exit()
-
-            gene_seq[gene]=str(r.text)
+            r=r.json()
+            if len(r)==1:
+                r=dict(r[0])
+                gene_seq[gene]=str(r["seq"])
                 return
-        except:
+            else:
+                maxi=0
+                maxlen=0
+                for i in range(len(r)):
+                    m=r[i]
+                    m=dict(m)
+                    if len(m["seq"])>maxlen:
+                        maxi=i
+                r=dict(r[maxi])
+                gene_seq[gene]=str(r["seq"])
+                return
+        except Exception as e:
+            print("\nError:",e)
             continue
 
 def create_synteny_matrix_mul(gene_seq,g1,g2,n):
@@ -32,7 +46,7 @@ def create_synteny_matrix_mul(gene_seq,g1,g2,n):
         try:
             temp=gene_seq[gene]
         except:
-            #print("Updating gene sequences for gene:",gene)
+            print("Updating gene sequences for gene:",gene)
             update(gene_seq,gene)
     for gene in g2:
         if gene=="NULL_GENE":
@@ -40,7 +54,7 @@ def create_synteny_matrix_mul(gene_seq,g1,g2,n):
         try:
             temp=gene_seq[gene]
         except:
-            #print("Updating gene sequences for gene:",gene)
+            print("Updating gene sequences for gene:",gene)
             update(gene_seq,gene)
     #print(n)
     sm=np.zeros((n,n,2))
@@ -61,6 +75,7 @@ def create_synteny_matrix_mul(gene_seq,g1,g2,n):
                 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)
+
             except:
                 return np.zeros((n,n,2)),np.zeros((n,n,2))
     return sm,sml
@@ -70,12 +85,12 @@ def synteny_matrix(gene_seq,hdf,lsy,n,enable_break):
     sl=[]
     t=0
     ind=[]
-    start=time.time()
     for index,row in progressbar.progressbar(hdf.iterrows()):
         g1=str(row["gene_stable_id"])
         g2=str(row["homology_gene_stable_id"])
         x=[]
         y=[]
+        t+=1
         try:
             temp=lsy[g1]
             temp=lsy[g2]
@@ -101,10 +116,9 @@ def synteny_matrix(gene_seq,hdf,lsy,n,enable_break):
         sg.append(smgtemp)
         sl.append(smltemp)
         ind.append(index)
-        t+=1
         if t==5 and enable_break==1:
             break
-    end=time.time()
     #print("Time Taken:",end-start)
     #print("Average Time:",(end-start)/len(sg))
+    print(t)
     return np.array(sg),np.array(sl),np.array(ind)

open_synteny_matrices.py

@@ -28,6 +28,7 @@ while(1):
         synteny_matrices_global=np.load(ndir+str(d_h[ch-1])+nf1+".npy")
         synteny_matrices_local=np.load(ndir+str(d_h[ch-1])+nf2+".npy")
         indexes=np.load(ndir+str(d_h[ch-1])+nf3+".npy")
+        print(len(indexes))
         break
     except:
         print("Choice invalid or incomplete files!!!!!. Try Another Index.")

open_synteny_matrices_negative.py

@@ -0,0 +1,112 @@
+import numpy as np 
+import pandas as pd
+import os
+import json
+import gc
+from read_data import read_data_homology
+import matplotlib.pyplot as plt
+import seaborn as sns
+from matplotlib.backends.backend_pdf import PdfPages
+
+a_h=[]
+d_h=[]
+df=pd.read_hdf("negative_dataset.h5",key="ndf")
+
+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))
+print("Neighbor Genes Loaded")
+
+synteny_matrices_global=np.load(ndir+"negative_dataset"+nf1+str(6)+".npy")
+synteny_matrices_local=np.load(ndir+"negative_dataset"+nf2+str(6)+".npy")
+indexes=np.load(ndir+"negative_dataset"+nf3+str(6)+".npy")
+df=df.loc[indexes]
+
+print(indexes[0:1000])
+ng=["Levenshtein Distance","Levenshtein Distance Reverse"]
+nl=["Local Alignment Score","Local Alignment Score Reverse"]
+
+inddict={}
+for i in range(len(indexes)):
+    inddict[indexes[i]]=i
+
+while(1):
+    try:
+        i=int(input("Enter the index:"))
+    except:
+        break
+    if i in inddict:
+        font = {'family': 'sans-serif',
+        'color':  'darkturquoise',
+        'weight': 'heavy',
+        'size': 20,
+        }
+        pdf=PdfPages(str(i)+".pdf")
+        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)
+        text="Species:"+df.loc[i].species
+        text+="\n"+"Gene Stable Id:"+df.loc[i].gene_stable_id
+        text+="\n"+"Homology Species:"+df.loc[i].homology_species
+        text+="\n"+"Homology Gene Stable Id:"+df.loc[i].homology_gene_stable_id
+        text+="\n"+"Homology Type:"+"Non Homology"
+        fp=plt.figure(figsize=(10,10))
+        fp.text(0.5,0.5,text,ha="center",fontdict=font)
+        pdf.savefig()
+        plt.close()
+        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)            
+            fig,ax=plt.subplots(figsize=(8,8))
+            ax.set_xlabel(str(df.loc[i].homology_species),fontsize=10)
+            ax.set_ylabel(str(df.loc[i].species),fontsize=10)
+            hmap=sns.heatmap(matrix,xticklabels=y, yticklabels=x,annot=True,ax=ax,linewidths=.5,cmap="YlGnBu",annot_kws={"size": 10})
+            hmap.figure.subplots_adjust(left=0.33,bottom=0.33,right=0.79,top=0.79)
+            ax.set_title(ng[m])
+            #plt.text(1,0.5,text,size=10)
+            pdf.savefig()
+            plt.show()
+
+        print("Local Alignment Matrix")
+        for m in range(sg.shape[-1]):
+            matrix=sl[:,:,m]
+            print(matrix)
+            fig,ax=plt.subplots(figsize=(8,8))
+            ax.set_xlabel(df.loc[i].homology_species)
+            ax.set_ylabel(df.loc[i].species)
+            hmap=sns.heatmap(matrix,xticklabels=y, yticklabels=x,annot=True,ax=ax,linewidths=.5,cmap="YlGnBu",annot_kws={"size": 10})
+            hmap.figure.subplots_adjust(left=0.27,bottom=0.29,right=0.92)
+            ax.set_title(nl[m])
+            pdf.savefig()
+            plt.show()
+        pdf.close()
+    else:
+        print("Index not found")
+

prepare_negative_dataset.py

@@ -15,13 +15,6 @@ seed=int(arg[1])
 random.seed(seed)
 
 a_h,d_h=read_data_homology("data_homology")
-a=[]
-d={}
-a,d=read_data_genome("data",a,d)
-assert(len(a)==len(d))
-indexes_gd=[list(df.index.values) for df in a]
-d=list(d.keys())
-
 gmap={}
 indexes_hd=[]
 for df in progressbar.progressbar(a_h):
@@ -30,35 +23,53 @@ for df in progressbar.progressbar(a_h):
     for h in hgids:
         gmap[h]=1
 
+a=[]
+d={}
+a,d=read_data_genome("data",a,d)
+assert(len(a)==len(d))
+indexes_gd=[list(df.index.values) for df in a]
+d=list(d.keys())
+ld=[]
+for i in progressbar.progressbar(range(len(a))):
+    df=a[i]
+    ldg=[]
+    for _,row in df.iterrows():
+        gid=row.gene_id
+        try:
+            temp=gmap[gid]
+        except:
+            ldg.append(gid)
+    ld.append(ldg)
+assert(len(ld)==len(a))
 
-
-col_names=["gid","species","hgid","h_species"]
+negativesamp={}
 nohd=len(a_h)
 nogd=len(a)
 rows=[]
 for i in progressbar.progressbar(range(nos)):
     while(1):
         try:
-            slh=random.randrange(nohd)
-            slg=random.randrange(nogd)
-            slgd=a[slg]
-            n1=d[slg]
-            indexes=indexes_gd[slg]
-            ind=random.randrange(len(indexes))
+            slh=random.randrange(nohd)#sample a homology database
+            slg=random.randrange(nogd)#sample a gene annotation file
+            slhd=a_h[slh]#select the homology database
+            indexes=indexes_hd[slh]#select the respective indexes
+            ldg=ld[slg]#select the given gene_id annotations
+            ind1=random.randrange(len(ldg))#sample a gene
+            g1=ldg[ind1]#get the gene id
+            ind2=random.randrange(len(indexes))#sample a row
+            row=slhd.loc[indexes[ind2]]#get the row from the database
+            r={}
             try:
-                _=gmap[slgd.loc[indexes[ind]].gene_id]
+                _=negativesamp[row.gene_stable_id+g1]#check if they exist in the database
+                _=negativesamp[g1+row.gene_stable_id]#check if they exist in the database
                 continue
             except:
-                slhd=a_h[slh]
-                lid=indexes_hd[slh]
-                ind_1=random.randrange(len(slhd))
-                row=slhd.loc[lid[ind_1]]
-                r={}
-                r["gene_stable_id"]=row.gene_stable_id 
-                r["species"]=row.species
-                r["homology_gene_stable_id"]=slgd.loc[indexes[ind]].gene_id
-                r["homology_species"]=n1  
-                rows.append(r) 
+                r["gene_stable_id"]=row.gene_stable_id#add it to the row
+                r["species"]=row.species#add speccies to the row
+                r["homology_gene_stable_id"]=g1#add the gene_id to the row
+                r["homology_species"]=d[slg]#add the gene species to the row
+                rows.append(r)#add it to the rows dict
+                negativesamp[row.gene_stable_id+g1]=1#add to the map so duplicate samples are avoided
                 break 
         except:
             continue  

prepare_other_factors.py

@@ -0,0 +1,55 @@
+import pandas as pd
+import numpy as np 
+import json
+import gc
+import pickle
+from read_data import read_data_homology
+from tree_data import create_tree_data
+
+a_h,d_h=read_data_homology("data_homology")
+d_h=list(d_h.keys())
+
+smg_name="_synteny_matrices_global.npy"
+sml_name="_synteny_matrices_local.npy"
+smi_name="_indexes.npy"
+dir_name="processed/synteny_matrices/"
+label=dict(ortholog_one2one=0,other_paralog=1,ortholog_one2many=1,ortholog_many2many=1,within_species_paralog=0)
+
+rows=[]
+for i in range(len(a_h)):
+    df=a_h[i]
+    n=d_h[i]
+    try:
+        smg=np.load(dir_name+n+smg_name)
+        sml=np.load(dir_name+n+sml_name)
+        indexes=np.load(dir_name+n+smi_name)
+    except:
+        print("Incomplete data for:",n)
+        continue
+    df=df.loc[indexes]
+    branch_length_species,branch_length_homology_species,distance,dist_p_s,dist_p_hs=create_tree_data("species_tree.tree",df)
+    assert(len(branch_length_species)==len(df))
+    assert(len(sml)==len(distance))
+    for i in range(len(df)):
+        index=indexes[i]
+        row=df.loc[index]
+        r={}
+        r["species"]=row["species"]
+        r["homology_species"]=row["homology_species"]
+        r["gene_stable_id"]=row["gene_stable_id"]
+        r["homology_gene_stable_id"]=row["homology_gene_stable_id"]
+        r["label"]=label[row["homology_type"]]
+        r["global_alignment_matrix"]=smg[i]
+        r["local_alignment_matrix"]=sml[i]
+        r["index_homology_dataset"]=index
+        r["bls"]=branch_length_species[i]
+        r["blhs"]=branch_length_homology_species[i]
+        r["dis"]=distance[i]
+        r["dps"]=dist_p_s[i]
+        r["dphs"]=dist_p_hs[i]
+        rows.append(r)
+
+with open("dataset","wb") as file:
+    pickle.dump(rows,file)
+
+print("Data Saved Successfully:)")

prepare_synteny_matrix.py

@@ -7,6 +7,7 @@ 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
+from access_data_rest import update_rest
 
 if not os.path.isdir("processed/synteny_matrices"):
     os.mkdir("processed/synteny_matrices")
@@ -29,6 +30,8 @@ for i in range(len(a_h)):
     assert(len(df)==nos)
     a_h[i]=df
 
+print("Data Selected")
+
 with open("processed/neighbor_genes.json","r") as file:
     lsy=dict(json.load(file))
 print(len(lsy))
@@ -38,6 +41,8 @@ print("Neighbor Genes Loaded")
 gene_sequences=read_gene_sequences(a_h,lsy,"geneseq","gene_sequences")
 print("Gene Sequences Loaded")
 
+print("Going to update not found sequences:")
+gene_sequences=update_rest(gene_sequences)
 
 n=3
 ndir="processed/synteny_matrices/"
@@ -46,9 +51,10 @@ nf2="synteny_matrices_local"
 nf3="indexes"
 for i in range(len(a_h)):
     df=a_h[i]
+    print(len(df))
     synteny_matrices_global,synteny_matrices_local,indexes=synteny_matrix(gene_sequences,df,lsy,n,0)
     np.save(ndir+str(d_h[i])+"_"+nf1,synteny_matrices_global)
     np.save(ndir+str(d_h[i])+"_"+nf2,synteny_matrices_local)
     np.save(ndir+str(d_h[i])+"_"+nf3,indexes)
-        
+    print(len(indexes))    
 print("Synteny Matrices Created Successfully :)")

prepare_synteny_matrix_negative.py

@@ -0,0 +1,42 @@
+import pandas as pd 
+import numpy as np 
+import os 
+import sys
+import progressbar
+import json
+from read_get_gene_seq import read_gene_sequences
+from create_synteny_matrix import synteny_matrix
+
+if not os.path.isdir("processed/synteny_matrices"):
+    os.mkdir("processed/synteny_matrices")
+
+df=pd.read_hdf("negative_dataset.h5",key="ndf")
+for _,row in progressbar.progressbar(df.iterrows()):
+    row["homology_species"]=row["homology_species"].lower()
+print(df[0:10])
+
+with open("processed/neighbor_genes.json","r") as file:
+    lsy=dict(json.load(file))
+print(len(lsy))
+print("Neighbor Genes Loaded")
+
+a_h=[]
+a_h.append(df)
+#gene_sequences=read_gene_sequences(a_h,lsy,"geneseq","gene_sequences")
+with open("processed/gene_sequences.json","r") as file:
+    gene_sequences=dict(json.load(file))
+print("Gene Sequences Loaded")
+
+
+
+n=3
+ndir="processed/synteny_matrices/"
+nf1="synteny_matrices_global"
+nf2="synteny_matrices_local"
+nf3="indexes"
+for i in range(10):
+    synteny_matrices_global,synteny_matrices_local,indexes=synteny_matrix(gene_sequences,df[i*100000:(i+1)*100000],lsy,n,0)
+    np.save(ndir+"negative_dataset"+"_"+nf1+str(i),synteny_matrices_global)
+    np.save(ndir+"negative_dataset"+"_"+nf2+str(i),synteny_matrices_local)
+    np.save(ndir+"negative_dataset"+"_"+nf3+str(i),indexes)
+print("Synteny Matrices Created Successfully :)")

process_data.py

@@ -147,7 +147,7 @@ def create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,cmap,cimap,save_after,enable_br
                     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 Exception as e:
+                except:
                     continue
             try:
                 z=lsy[y]

read_get_gene_seq.py

@@ -108,8 +108,16 @@ def read_gene_sequences(hdf,lsy,data_dir,fname):
     s=[x.capitalize() for x in s]
     
     data=read_gene_seq(data_dir,s,gene_by_species_dict)
+    not_found={}
+    for species in gene_by_species_dict:
+        for gene in gene_by_species_dict[species]:
+            try:
+                _=data[gene]
+            except:
+                not_found[gene]=1
 
+    with open("processed/not_found.json","w") as file:
+        json.dump(not_found,file)
     with open("processed/"+fname+".json","w") as file:#save the data
         json.dump(data,file)
-
     return data

read_negative_dataset.py

@@ -4,3 +4,4 @@ df=pd.read_hdf("negative_dataset.h5",key="ndf")
 
 print(df.info())
 print(df.loc[33333:33433])
+print(df[df["species"]=="mola_mola"])

tree_data.py

@@ -1,5 +1,16 @@
 from ete3 import Tree
 import numpy as np
+import progressbar
+
+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 create_tree_data(treename,df):
     t=Tree(treename)
@@ -8,7 +19,7 @@ def create_tree_data(treename,df):
     dist=[]
     ns=[]
     nhs=[]
-    for index,row in df.iterrows():
+    for index,row in progressbar.progressbar(df.iterrows()):
         d=0
         x=row["species"]
         y=row["homology_species"]
@@ -34,4 +45,6 @@ def create_tree_data(treename,df):
         nhs.append(c)
         branch_lengths_hs.append(bl)
         dist.append(d)
+    create_branch_length_padding(branch_lengths_s)
+    create_branch_length_padding(branch_lengths_hs)
     return np.array(branch_lengths_s),np.array(branch_lengths_hs),np.array(dist),np.array(ns),np.array(nhs)

update_neighbor_genes_ndf.py

@@ -0,0 +1,42 @@
+import sys
+import numpy as np
+import pandas as pd
+import json
+import os
+import gc
+from get_data import get_data_homology,get_data_genome
+from process_data import create_data_homology_ls
+
+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"
+cmap,cimap,ld,ldg,a,d=get_data_genome(arg,dir_g)
+
+df=pd.read_hdf("negative_dataset.h5",key="ndf")
+a_h=[]
+d_h=[]
+a_h.append(df)
+d_h.append("negative_dataset")
+print("Data Read")
+
+n=3 #no. of numbers neighbors
+save_after=5 #to save data after n steps
+
+lsy=create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,cmap,cimap,save_after,enable_break,1)
+print(len(lsy))
+print("Neighbor Genes Updated Successfully")