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/create_synteny_matrix_v2.py b/create_synteny_matrix_v2.py new file mode 100644 index 0000000..2efe511 --- /dev/null +++ b/create_synteny_matrix_v2.py @@ -0,0 +1,119 @@ +import numpy as np +import requests +import edlib as ed +import pandas as pd +import time +import sys +import progressbar +from skbio.alignment import local_pairwise_align_ssw +from skbio import DNA,TabularMSA,RNA + +def update(gene_seq,gene): + while(1): + try: + server = "https://rest.ensembl.org" + ext = "/sequence/id/"+str(gene)+"?type=cds;multiple_sequences=1" + + r = requests.get(server+ext, headers={ "Content-Type" : "application/json"}) + + if not r.ok: + r.raise_for_status() + sys.exit() + 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): + 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=max(len(gene_seq[g1[i]]),len(gene_seq[g2[j]])) + try: + 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) + except: + return np.zeros((n,n,2)),np.zeros((n,n,2)) + return sm,sml + +def synteny_matrix(gene_seq,hdf,lsy,n,enable_break): + sg=[] + sl=[] + t=0 + ind=[] + 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] + 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) + if np.all(smgtemp==0): + continue + sg.append(smgtemp) + sl.append(smltemp) + ind.append(index) + if t==5 and enable_break==1: + break + #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/get_data.py b/get_data.py index a2ee986..96634f4 100644 --- a/get_data.py +++ b/get_data.py @@ -2,7 +2,7 @@ import sys import os 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 +from process_data import list_dict_genomes,create_chromosome_maps def get_data_genome(arg,dir): a=[] @@ -22,15 +22,16 @@ def get_data_genome(arg,dir): if arg[4]=="-r": a,d=read_data_genome(dir,a,d) assert(len(a)==len(d)) - + print("Creating Maps:") ld,ldg=list_dict_genomes(a,d) + cmap,cimap=create_chromosome_maps(a,d) assert(len(ld)==len(ldg)) for i in range(len(ld)): assert(len(ld[i])==len(ldg[i])) - return ld,ldg,a,d + return cmap,cimap,ld,ldg,a,d def get_data_homology(arg,dir): a_h=[] diff --git a/main.py b/main.py deleted file mode 100644 index 47d56e2..0000000 --- a/main.py +++ /dev/null @@ -1,53 +0,0 @@ -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=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)) - -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],lsy,n) -print("Synteny Matrices are created successfully\n",len(indexes),"\n",len(synteny_matrices)) -np.save("synteny_matrices",synteny_matrices) -np.save("indexes",indexes) - -with open("gene_seq_updated.json","w") as file: - json.dump(gene_sequences,file) - -df=a_h[0].loc[indexes] -branch_length_species,branch_length_homology_species,distance,dist_p_s,dist_p_hs=create_tree_data("species_tree.tree",df) - -train_synteny_matrices,train_branch_length_species,train_branch_length_homology_species,train_mean_gene_length,train_dist_p_s,train_dist_p_hs,train_distance,train_labels=train_data(indexes,synteny_matrices,df,branch_length_species,branch_length_homology_species,distance,dist_p_s,dist_p_hs,gene_sequences) - -train(train_synteny_matrices,train_branch_length_species,train_branch_length_homology_species,train_mean_gene_length,train_dist_p_s,train_dist_p_hs,train_distance,train_labels) diff --git a/my.txt b/my.txt deleted file mode 100644 index 6263a18..0000000 --- a/my.txt +++ /dev/null @@ -1,3 +0,0 @@ -ftp://ftp.ensembl.org/pub/release-96/gtf/lepisosteus_oculatus/Lepisosteus_oculatus.LepOcu1.96.gtf.gz -ftp://ftp.ensembl.org/pub/release-96/gtf/mola_mola/Mola_mola.ASM169857v1.96.gtf.gz -ftp://ftp.ensembl.org/pub/release-96/gtf/homo_sapiens/Homo_sapiens.GRCh38.96.gtf.gz 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_data.py b/prepare_data.py deleted file mode 100644 index b5ea833..0000000 --- a/prepare_data.py +++ /dev/null @@ -1,87 +0,0 @@ -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,enable_break) -else: - synteny_matrices_global,synteny_matrices_local,indexes=synteny_matrix(gene_sequences,a_h[0],lsy,n,enable_break) - -print("Synteny Matrices are created successfully\n",len(indexes),"\n",len(synteny_matrices_global)) -np.save("processed/synteny_matrices_global",synteny_matrices_global) -np.save("processed/synteny_matrices_local",synteny_matrices_local) -np.save("processed/indexes",indexes) - -with open("processed/gene_seq_updated.json","w") as file: - json.dump(gene_sequences,file) - -df=a_h[0].loc[indexes] -branch_length_species,branch_length_homology_species,distance,dist_p_s,dist_p_hs=create_tree_data("species_tree.tree",df) - -train_synteny_matrices_global,train_synteny_matrices_local,train_branch_length_species,train_branch_length_homology_species,train_mean_gene_length,train_dist_p_s,train_dist_p_hs,train_distance,train_labels=train_data(indexes,synteny_matrices_global,synteny_matrices_local,df,branch_length_species,branch_length_homology_species,distance,dist_p_s,dist_p_hs,gene_sequences) - -np.save("processed/train_synteny_matrices_global",train_synteny_matrices_global) -np.save("processed/train_synteny_matrices_local",train_synteny_matrices_local) -np.save("processed/train_branch_length_species",train_branch_length_species) -np.save("processed/train_branch_length_homology_species",train_branch_length_homology_species) -np.save("processed/train_mean_gene_length",train_mean_gene_length) -np.save("processed/train_dist_p_s",train_dist_p_s) -np.save("processed/train_dist_p_hs",train_dist_p_hs) -np.save("processed/train_distance",train_distance) -np.save("processed/train_labels",train_labels) - -print("Data Saved Successfully to processed :)") diff --git a/prepare_negative_dataset.py b/prepare_negative_dataset.py new file mode 100644 index 0000000..151ebec --- /dev/null +++ b/prepare_negative_dataset.py @@ -0,0 +1,82 @@ +import sys +import numpy as np +import pandas as pd +import os +import gc +import progressbar +import random +from read_data import read_data_homology,read_data_genome + +arg=sys.argv +arg=arg[1:] + +nos=int(arg[0]) +seed=int(arg[1]) +random.seed(seed) + +a_h,d_h=read_data_homology("data_homology") +gmap={} +indexes_hd=[] +for df in progressbar.progressbar(a_h): + indexes_hd.append(list(df.index.values)) + hgids=df.homology_gene_stable_id.unique() + 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)) + +negativesamp={} +nohd=len(a_h) +nogd=len(a) +rows=[] +for i in progressbar.progressbar(range(nos)): + while(1): + try: + 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: + _=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: + 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 + +a_h=[] +gc.collect() + +ndf=pd.DataFrame(rows) +print(ndf[0:10]) +ndf.to_hdf("negative_dataset.h5",key="ndf",mode="w") 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 8a2cee1..b4da37a 100644 --- a/process_data.py +++ b/process_data.py @@ -16,6 +16,26 @@ def create_map_list(l): #this function maps the indexes to values t[l[i]]=i return t +def create_chromosome_maps(a,n): + cmap=[] + cimap=[] + for df in progressbar.progressbar(a): + chmap={} + chindmap={} + for index,row in df.iterrows(): + g=row.gene_id + try: + temp=chmap[g] + except: + chmap[g]=str(row.Chr) + if str(row.Chr) in chindmap: + chindmap[str(row.Chr)].append(index) + else: + chindmap[str(row.Chr)]=[] + chindmap[str(row.Chr)].append(index) + cmap.append(chmap) + cimap.append(chindmap) + return cmap,cimap def list_dict_genomes(a,n): lst=[] @@ -23,17 +43,20 @@ def list_dict_genomes(a,n): for x in a: ldgt={} uc=list(x["gene_id"]) - for i in range(len(uc)): - ldgt[uc[i]]=i + for i,r in x.iterrows(): + ldgt[r.gene_id]=i lst.append(uc) ldt.append(ldgt) return lst,ldt -def get_nearest_neighbors(g,gs,n,a,d,ld,ldg): +def get_nearest_neighbors(g,gs,n,a,d,ld,ldg,cmap,cimap): + #print("Finding Neighbor Genes") ne=[] #list to store the backward genes nr=[] #list to store the forward genes gi=d[gs.capitalize()] #get the address of the corresponding species to which the gene belongs whose neighbor has to be found sldf=a[gi]#select the dataframe + scmap=cmap[gi]#select the correct chromosome map + scimap=cimap[gi]#select the correct index maps try: sld=ld[gi]#see if the corresponding gene map exists except: @@ -42,9 +65,12 @@ 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 genes + i=sldg[g]#find the index of the gene + chromosome_id=scmap[g]#get the chromosome no from the database. + scimap=scimap[chromosome_id]#select the correct chromosomes indexes + sldf=sldf.loc[scimap]#select only the same chromosome genes. #get the -n neighbors - start=int(sldf.iloc[i,[3]])#get the start location of the gene + start=int(sldf.loc[i]['start'])#get the start location of the gene flag=0 for j in range(n): if flag==1: @@ -52,9 +78,9 @@ def get_nearest_neighbors(g,gs,n,a,d,ld,ldg): continue itemp=0 #select the column - end=sldf.iloc[:,4] + end=list(sldf.end) end=np.array(end) - assert(len(end)==len(sld)) + assert(len(end)==len(sldf)) 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 consideration @@ -65,18 +91,19 @@ def get_nearest_neighbors(g,gs,n,a,d,ld,ldg): 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])#push the gene in the array - start=int(sldf.iloc[itemp,[3]])#make "start" the start location of the current gene + itemp=scimap[itemp] + ne.append(sldf.loc[itemp].gene_id) + start=int(sldf.loc[itemp].start)#make "start" the start location of the current gene #print(start) #get the +n neighbors flag=0 - end=int(sldf.iloc[i,[4]]) + end=int(sldf.loc[i].end) for j in range(n): if flag==1: nr.append("NULL_GENE") continue itemp=0 - start=sldf.iloc[:,3] + start=list(sldf.start) start=np.array(start) start=start-end start_s=np.argsort(start) @@ -88,12 +115,12 @@ def get_nearest_neighbors(g,gs,n,a,d,ld,ldg): if start[k]>0: itemp=k break - nr.append(sld[itemp]) - end=int(sldf.iloc[itemp,[4]]) - + itemp=scimap[itemp] + nr.append(sldf.loc[itemp].gene_id) + end=int(sldf.loc[itemp].end) return ne,nr -def create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,save_after,enable_break,update): +def create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,cmap,cimap,save_after,enable_break,update): 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 if os.path.exists("processed/neighbor_genes.json"): @@ -101,10 +128,11 @@ def create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,save_after,enable_break,update) lsy=dict(json.load(file)) print("Existing neighbor genes read!!") print(len(lsy)) - c=4 + c=0 lsytemp={} name="neighbor_genes" for df in a_h: + #df=df[0:2] for _,row in progressbar.progressbar(df.iterrows()): x=row["gene_stable_id"] y=row["homology_gene_stable_id"] @@ -115,7 +143,7 @@ def create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,save_after,enable_break,update) except: try: t2=d[xs.capitalize()]#see if the species exist in genomic maps - xl,xr=get_nearest_neighbors(x,xs,n,a,d,ld,ldg) + xl,xr=get_nearest_neighbors(x,xs,n,a,d,ld,ldg,cmap,cimap) 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) @@ -126,7 +154,7 @@ def create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,save_after,enable_break,update) except: try: t2=d[ys.capitalize()] - yl,yr=get_nearest_neighbors(y,ys,n,a,d,ld,ldg) + yl,yr=get_nearest_neighbors(y,ys,n,a,d,ld,ldg,cmap,cimap) if len(yl)!=0: lsy[y]=dict(b=yl,f=yr) lsytemp[y]=dict(b=yl,f=yr) diff --git a/read_data.py b/read_data.py index ba43f96..f599dd7 100644 --- a/read_data.py +++ b/read_data.py @@ -2,6 +2,7 @@ import os import pandas as pd import gzip import sys +import progressbar def clear_data(x): x=x.split() @@ -18,23 +19,23 @@ def read_data_genome(dir_name,a,dict_ind_genome): if len(lf)==0: print("No files in the data directory!!!!!!") sys.exit(1) + colname=["Chr","source","feature","start","end","score","strand","frame","attribute"] print("Going to read data:") - for x in range(len(lf)): - data_gene=pd.read_csv(dir_name+"/"+lf[x],compression='gzip',sep='\t',comment='#',header=None) + for x in progressbar.progressbar(range(len(lf))): + data_gene=pd.read_csv(dir_name+"/"+lf[x],compression='gzip',sep='\t',comment='#',header=None,names=colname) #print(data_gene.head) - data_gene=data_gene[data_gene[2]=="gene"] - data_gene=data_gene.sort_values(3) - tmp=data_gene[8].str.split(";",expand=True) + data_gene=data_gene[data_gene["feature"]=="gene"] + tmp=data_gene["attribute"].str.split(";",expand=True) tmp=tmp.iloc[:,:5] data_gene[["gene_id","gene_version","gene_name","gene_source","gene_biotype"]]=tmp - data_gene=data_gene.drop(8,axis=1) + data_gene=data_gene.drop("attribute",axis=1) #print(data_gene[0:10]) try: for y in ["gene_version","gene_name","gene_source","gene_biotype","gene_id"]: data_gene[y]=data_gene[y].apply(clear_data) except: continue - #print(data_gene[0:10]) + #print(data_gene[0:10]) data_gene=data_gene[data_gene['gene_biotype']=='protein_coding'] a.append(data_gene) n=lf[x].split(".")[0] 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 new file mode 100644 index 0000000..df1068f --- /dev/null +++ b/read_negative_dataset.py @@ -0,0 +1,7 @@ +import pandas as pd + +df=pd.read_hdf("negative_dataset.h5",key="ndf") + +print(df.info()) +print(df.loc[33333:33433]) +print(df[df["species"]=="mola_mola"]) diff --git a/test_main.py b/test_main.py deleted file mode 100644 index dddb487..0000000 --- a/test_main.py +++ /dev/null @@ -1,54 +0,0 @@ -import sys -import numpy as np -import pandas as pd -import json -from get_data import get_data_homology,get_data_genome -from process_data import create_data_homology_ls -from read_get_gene_seq import read_gene_sequences -from create_synteny_matrix import synteny_matrix -from tree_data import create_tree_data -from create_train_data import train_data -from train import train - -arg=sys.argv -arg=arg[1:] - -if len(arg)!=5: - print("No. of arguments more or less. Please check") - sys.exit(1) - -dir_g="data" -ld,ldg,a,d=get_data_genome(arg,dir_g) - -dir_hom="data_homology" -a_h,d_h=get_data_homology(arg,dir_hom) - -if arg[-1]=="-d": - sys.exit(1) - -print("Data Read") - -n=2 #no. of numbers neighbors -save_after=3 #to save data after n steps -lsy=create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,save_after) -print(len(lsy)) - -print("Neighbor Genes Found") - -gene_sequences=read_gene_sequences(a_h[0],lsy,"geneseq","gene_sequences") - -synteny_matrices,indexes=synteny_matrix(gene_sequences,a_h[0][0:1000],lsy,n) -print("Synteny Matrices are created successfully\n",len(indexes),"\n",len(synteny_matrices)) -np.save("synteny_matrices",synteny_matrices) -np.save("indexes",indexes) - -with open("gene_seq_updated.json","w") as file: - json.dump(gene_sequences,file) - -df=a_h[0].loc[indexes] -branch_length_species,branch_length_homology_species,distance,dist_p_s,dist_p_hs=create_tree_data("species_tree.tree",df) - -train_synteny_matrices,train_branch_length_species,train_branch_length_homology_species,train_mean_gene_length,train_dist_p_s,train_dist_p_hs,train_distance,train_labels=train_data(indexes,synteny_matrices,df,branch_length_species,branch_length_homology_species,distance,dist_p_s,dist_p_hs,gene_sequences) - -train(train_synteny_matrices,train_branch_length_species,train_branch_length_homology_species,train_mean_gene_length, - train_dist_p_s,train_dist_p_hs,train_distance,train_labels) diff --git a/test_model.py b/test_model.py new file mode 100644 index 0000000..62ab0e6 --- /dev/null +++ b/test_model.py @@ -0,0 +1,181 @@ +import tensorflow as tf +import numpy as np +import pandas as pd +import requests +import sys +import pickle +from get_data import get_data_genome +from process_data import get_nearest_neighbors +from create_synteny_matrix_v2 import create_synteny_matrix_mul +from ete3 import Tree +import traceback + +""" +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)""" + +data={} +with open("genome_maps","rb") as file: + data=pickle.load(file) +cmap=data["cmap"] +cimap=data["cimap"] +ld=data["ld"] +ldg=data["ldg"] +a=data["a"] +d=data["d"] + +n=3 + +server = "https://rest.ensembl.org" +def get_gene_data(gid): + empty={} + ext = "/lookup/id/"+gid+"?" + try: + r = requests.get(server+ext, headers={ "Content-Type" : "application/json"}) + if not r.ok: + r.raise_for_status() + except: + return empty + decoded = r.json() + return dict(decoded) + +def create_branch_length_padding(bl): + maxlen=29 + for x in bl: + for i in range(len(x),maxlen): + x.append(0) + +def create_tree_data(x,y): + t=Tree("species_tree.tree") + branch_lengths_s=[] + branch_lengths_hs=[] + dist=[] + ns=[] + nhs=[] + d=0 + 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) + 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) + +def filter_gene_details(gid,n,a,d,ld,ldg,cmap,cimap): + gd=get_gene_data(gid) + if len(gd)==0 or gd["biotype"]!="protein_coding": + print("Gene Id is Incorrect!!!. Try again.") + return list() + sp=gd["species"] + print(sp) + try: + t2=d[sp.capitalize()] + except: + print("Gene Species does not exist in gene maps. Try another species.") + return list(),sp + nl,nr=get_nearest_neighbors(gid,sp,n,a,d,ld,ldg,cmap,cimap) + if len(nl)==0 or len(nr)==0: + print("Gene is not present in gtf file. Please Try another gene.") + return list(),sp + print("Found nearest neighbors") + x=[] + for i in range(len(nl)-1,-1,-1): + x.append(nl[i]) + x.append(gid) + for k in nr: + x.append(k) + assert(len(x)==2*n+1) + return x,sp + +def get_label(l): + if l==0: + return "Paralogous" + if l==1: + return "Orthologous" + return "Non Homologous" + +def softmax(arr): + arr=arr-np.amax(arr) + arr=np.exp(arr) + arr=arr/np.sum(arr) + return arr + +geneseq={} +try: + model=tf.train.import_meta_graph('saved_models/model.ckpt.meta') +except: + print("Something wrong with the model.") + sys.exit(1) +with tf.Session() as sess: + try: + model.restore(sess,"saved_models/model.ckpt") + graph = tf.get_default_graph() + synmgt,synmlt,blst,blhst,dpst,dphst,dist,lrt,yt=graph.get_collection("input_nodes") + predictions=graph.get_tensor_by_name("Predictions/BiasAdd:0") + print("Model Loaded Successfully :)") + except: + print(":(") + sys.exit() + while(1): + try: + ch=input("Do you want to enter a gene id [y/n]:") + if ch=='n': + break + g1=input("Enter the first gene id:") + g1,sp1=filter_gene_details(g1,3,a,d,ld,ldg,cmap,cimap) + if len(g1)==0: + continue + g2=input("Enter the second gene id:") + g2,sp2=filter_gene_details(g2,3,a,d,ld,ldg,cmap,cimap) + if len(g2)==0: + continue + smg,sml=create_synteny_matrix_mul(geneseq,g1,g2,2*n+1) + bls,blhs,dis,dps,dphs=create_tree_data(sp1,sp2) + fd={synmgt:smg.reshape((1,7,7,2)), + synmlt:sml.reshape((1,7,7,2)), + blst:bls.reshape((1,29)), + blhst:blhs.reshape((1,29)), + dpst:dps.reshape((1,1)), + dist:dis.reshape((1,1)), + dphst:dphs.reshape((1,1))} + preds=sess.run([predictions],feed_dict=fd) + label=np.argmax(preds) + print(get_label(label)) + print(softmax(preds)) + + except Exception as e: + print("Some Error Was There Try Again:(",e) + traceback.print_exc() + continue \ No newline at end of file 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.py b/update_neighbor_genes.py index 8e6fdfb..42bbf92 100644 --- a/update_neighbor_genes.py +++ b/update_neighbor_genes.py @@ -25,7 +25,7 @@ if len(arg)!=5: sys.exit(1) dir_g="data" -ld,ldg,a,d=get_data_genome(arg,dir_g) +cmap,cimap,ld,ldg,a,d=get_data_genome(arg,dir_g) dir_hom="data_homology" a_h,d_h=get_data_homology(arg,dir_hom) @@ -36,9 +36,9 @@ if arg[-1]=="-d": print("Data Read") n=3 #no. of numbers neighbors -save_after=0 #to save data after n steps +save_after=5 #to save data after n steps -lsy=create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,save_after,enable_break,1) +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") 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