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/process_data.py b/process_data.py index 8a2cee1..b09de42 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,18 +143,18 @@ 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) - except: + except Exception as e: continue try: z=lsy[y] 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/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")