Add files via upload

This commit is contained in:
HarshitGupta11 2019-06-13 15:01:55 +05:30 committed by GitHub
parent f3ca3a5870
commit 0d205b7743
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
4 changed files with 62 additions and 32 deletions

View file

@ -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=[]

View file

@ -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)

View file

@ -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]

View file

@ -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")