From d66638422f145fa234cbf4525b2566077bfa5e89 Mon Sep 17 00:00:00 2001 From: HarshitGupta11 <50410275+HarshitGupta11@users.noreply.github.com> Date: Wed, 12 Jun 2019 13:55:02 +0530 Subject: [PATCH] Add files via upload --- create_synteny_matrix.py | 7 +++--- prepare_synteny_matrix.py | 47 ++++++++++++++++++--------------------- process_data.py | 19 +++++++++------- read_get_gene_seq.py | 14 +++++++----- update_neighbor_genes.py | 47 +++++++++++++++++++++++++++++++++++++++ 5 files changed, 93 insertions(+), 41 deletions(-) create mode 100644 update_neighbor_genes.py diff --git a/create_synteny_matrix.py b/create_synteny_matrix.py index a65ce95..ca7b2ee 100644 --- a/create_synteny_matrix.py +++ b/create_synteny_matrix.py @@ -4,6 +4,7 @@ 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 @@ -31,7 +32,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": @@ -39,7 +40,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)) @@ -70,7 +71,7 @@ def synteny_matrix(gene_seq,hdf,lsy,n,enable_break): t=0 ind=[] start=time.time() - for index,row in hdf.iterrows(): + for index,row in progressbar.progressbar(hdf.iterrows()): g1=str(row["gene_stable_id"]) g2=str(row["homology_gene_stable_id"]) x=[] diff --git a/prepare_synteny_matrix.py b/prepare_synteny_matrix.py index f100691..109bd23 100644 --- a/prepare_synteny_matrix.py +++ b/prepare_synteny_matrix.py @@ -8,50 +8,47 @@ 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"): +if not os.path.isdir("processed/synteny_matrices"): os.mkdir("processed/synteny_matrices") arg=sys.argv arg=arg[1:] -enable_break=0 - -if arg[-1]=="-test": - enable_break=1 +nos=int(arg[0]) lsy={} a_h,d_h=read_data_homology("data_homology") +print(d_h) +d_h=list(d_h.keys()) print("Homology Data Read") +for i in range(len(a_h)): + df=a_h[i] + random_indexes=np.random.permutation(len(df)) + random_indexes=random_indexes[:nos] + df=df.loc[random_indexes] + assert(len(df)==nos) + a_h[i]=df + 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,lsy,"geneseq","gene_sequences") + + +gene_sequences=read_gene_sequences(a_h,lsy,"geneseq","gene_sequences") print("Gene Sequences Loaded") -if enable_break==1: - save_after=10 -else: - save_after=500000 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 +for i in range(len(a_h)): + df=a_h[i] + 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("Synteny Matrices Created Successfully :)") diff --git a/process_data.py b/process_data.py index 7d1c77a..8a2cee1 100644 --- a/process_data.py +++ b/process_data.py @@ -2,6 +2,7 @@ import pandas import gc import numpy as np import json +import os import progressbar from save_data import save_data_json from save_data import write_dict_json @@ -92,13 +93,14 @@ def get_nearest_neighbors(g,gs,n,a,d,ld,ldg): return ne,nr -def create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,save_after,enable_break): +def create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,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 - with open("processed/neighbor_genes.json","r") as file: - lsy=dict(json.load(file)) - print("Existing neighbor genes read!!") - print(len(lsy)) + if os.path.exists("processed/neighbor_genes.json"): + 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" @@ -131,14 +133,15 @@ def create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,save_after,enable_break): except: continue t+=1 - if t>=save_after: + if t>=save_after and update==0: t=0 c+=1 write_dict_json(name+str(c),"processed",lsytemp) lsytemp={} if enable_break==1: break - c+=1 - write_dict_json(name+str(c),"processed",lsytemp) + if not update: + c+=1 + write_dict_json(name+str(c),"processed",lsytemp) write_dict_json(name,"processed",lsy) return lsy diff --git a/read_get_gene_seq.py b/read_get_gene_seq.py index 05e22b1..490e0af 100644 --- a/read_get_gene_seq.py +++ b/read_get_gene_seq.py @@ -1,8 +1,10 @@ import json from Bio import SeqIO +import numpy as np import pandas as pd import os import gzip +import progressbar def read_from_multiple_lsy(lsyfl): lsy={} @@ -56,7 +58,7 @@ def read_gene_seq(dirname,s,genes_by_species): if f.split(".")[0] in s:#check whether the species is present in the species to read list. Will skip those species which are not present in the dataframe ftr.append(f) data={} - for f in ftr: + for f in progressbar.progressbar(ftr): species=f.split(".")[0].lower() with gzip.open(dirname+"/"+f,"rt") as file: record=SeqIO.parse(file,"fasta") @@ -74,12 +76,14 @@ def read_gene_sequences(hdf,lsy,data_dir,fname): Thus we don't have to read the same file multiple times.""" grouped_genes={} gene_by_species_dict={} - for df in hdf: + for df in progressbar.progressbar(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: + gene_by_species_dict[i]=[] + for i in df.species.unique(): + gene_by_species_dict[i]=[] + + for x in progressbar.progressbar(lsy): try: species=grouped_genes[x]#get the species except: diff --git a/update_neighbor_genes.py b/update_neighbor_genes.py new file mode 100644 index 0000000..8e6fdfb --- /dev/null +++ b/update_neighbor_genes.py @@ -0,0 +1,47 @@ +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" +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 + +lsy=create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,save_after,enable_break,1) +print(len(lsy)) + +print("Neighbor Genes Updated Successfully") + + +