diff --git a/create_synteny_matrix.py b/create_synteny_matrix.py index cbc0f11..2d6b356 100644 --- a/create_synteny_matrix.py +++ b/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) diff --git a/open_synteny_matrices.py b/open_synteny_matrices.py new file mode 100644 index 0000000..c784f94 --- /dev/null +++ b/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") + + diff --git a/prepare_synteny_matrix.py b/prepare_synteny_matrix.py new file mode 100644 index 0000000..7bf417c --- /dev/null +++ b/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 :)") diff --git a/process_data.py b/process_data.py index 33f373e..7d1c77a 100644 --- a/process_data.py +++ b/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"] diff --git a/read_get_gene_seq.py b/read_get_gene_seq.py index 71c9062..05e22b1 100644 --- a/read_get_gene_seq.py +++ b/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'] @@ -94,12 +98,11 @@ def read_gene_sequences(df,lsy,data_dir,fname): 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] - - 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