From 3679a68b13473142c48a24192a49422627f021ab Mon Sep 17 00:00:00 2001 From: HarshitGupta11 <50410275+HarshitGupta11@users.noreply.github.com> Date: Mon, 17 Jun 2019 15:47:47 +0530 Subject: [PATCH] Add files via upload --- access_data_rest.py | 49 +++++++++++++ create_synteny_matrix.py | 38 ++++++---- open_synteny_matrices.py | 1 + open_synteny_matrices_negative.py | 112 +++++++++++++++++++++++++++++ prepare_negative_dataset.py | 63 +++++++++------- prepare_other_factors.py | 55 ++++++++++++++ prepare_synteny_matrix.py | 8 ++- prepare_synteny_matrix_negative.py | 42 +++++++++++ process_data.py | 2 +- read_get_gene_seq.py | 10 ++- read_negative_dataset.py | 3 +- tree_data.py | 15 +++- update_neighbor_genes_ndf.py | 42 +++++++++++ 13 files changed, 397 insertions(+), 43 deletions(-) create mode 100644 access_data_rest.py create mode 100644 open_synteny_matrices_negative.py create mode 100644 prepare_other_factors.py create mode 100644 prepare_synteny_matrix_negative.py create mode 100644 update_neighbor_genes_ndf.py diff --git a/access_data_rest.py b/access_data_rest.py new file mode 100644 index 0000000..368e336 --- /dev/null +++ b/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 + diff --git a/create_synteny_matrix.py b/create_synteny_matrix.py index ca7b2ee..5695dc7 100644 --- a/create_synteny_matrix.py +++ b/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) - return - except: + r=r.json() + if len(r)==1: + r=dict(r[0]) + gene_seq[gene]=str(r["seq"]) + return + 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] @@ -92,7 +107,7 @@ def synteny_matrix(gene_seq,hdf,lsy,n,enable_break): 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) @@ -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) diff --git a/open_synteny_matrices.py b/open_synteny_matrices.py index ee0dd85..5eff2ba 100644 --- a/open_synteny_matrices.py +++ b/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.") diff --git a/open_synteny_matrices_negative.py b/open_synteny_matrices_negative.py new file mode 100644 index 0000000..a37d8a1 --- /dev/null +++ b/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") + diff --git a/prepare_negative_dataset.py b/prepare_negative_dataset.py index 50daa50..151ebec 100644 --- a/prepare_negative_dataset.py +++ b/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 diff --git a/prepare_other_factors.py b/prepare_other_factors.py new file mode 100644 index 0000000..927faaa --- /dev/null +++ b/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:)") \ No newline at end of file diff --git a/prepare_synteny_matrix.py b/prepare_synteny_matrix.py index 109bd23..6034189 100644 --- a/prepare_synteny_matrix.py +++ b/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 :)") diff --git a/prepare_synteny_matrix_negative.py b/prepare_synteny_matrix_negative.py new file mode 100644 index 0000000..909f58a --- /dev/null +++ b/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 :)") diff --git a/process_data.py b/process_data.py index b09de42..b4da37a 100644 --- a/process_data.py +++ b/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] diff --git a/read_get_gene_seq.py b/read_get_gene_seq.py index 490e0af..cf31437 100644 --- a/read_get_gene_seq.py +++ b/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 diff --git a/read_negative_dataset.py b/read_negative_dataset.py index c8ba596..df1068f 100644 --- a/read_negative_dataset.py +++ b/read_negative_dataset.py @@ -3,4 +3,5 @@ import pandas as pd df=pd.read_hdf("negative_dataset.h5",key="ndf") print(df.info()) -print(df.loc[33333:33433]) \ No newline at end of file +print(df.loc[33333:33433]) +print(df[df["species"]=="mola_mola"]) diff --git a/tree_data.py b/tree_data.py index 8a19d55..035140f 100644 --- a/tree_data.py +++ b/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) diff --git a/update_neighbor_genes_ndf.py b/update_neighbor_genes_ndf.py new file mode 100644 index 0000000..fd30397 --- /dev/null +++ b/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") \ No newline at end of file