diff --git a/create_synteny_matrix.py b/create_synteny_matrix.py index 5695dc7..0fe8370 100644 --- a/create_synteny_matrix.py +++ b/create_synteny_matrix.py @@ -9,7 +9,8 @@ from skbio.alignment import local_pairwise_align_ssw from skbio import DNA,TabularMSA,RNA def update(gene_seq,gene): - while(1): + t=0 + while(t!=100): try: server = "https://rest.ensembl.org" ext = "/sequence/id/"+str(gene)+"?type=cds;multiple_sequences=1" @@ -36,8 +37,10 @@ def update(gene_seq,gene): gene_seq[gene]=str(r["seq"]) return except Exception as e: + t+=1 print("\nError:",e) continue + gene_seq[gene]="" def create_synteny_matrix_mul(gene_seq,g1,g2,n): for gene in g1: @@ -46,7 +49,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": @@ -54,7 +57,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)) @@ -75,7 +78,6 @@ 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 diff --git a/get_data.py b/get_data.py index 40dd959..7127274 100644 --- a/get_data.py +++ b/get_data.py @@ -36,21 +36,21 @@ def get_data_genome(arg,dir): def get_data_homology(arg,dir): a_h=[] d_h={} - if arg[2]=="-l": + if arg[0]=="-l": if not os.path.exists(dir): os.mkdir(dir) - download_data(arg[3],dir) - elif arg[2]=="-f": - get_data_file(arg[3],dir) - elif arg[2]=="-d": - if arg[4]=="-r": + download_data(arg[1],dir) + elif arg[0]=="-f": + get_data_file(arg[1],dir) + elif arg[0]=="-d": + if arg[2]=="-r": c=0 else: return a_h,d_h elif arg[2]=="-nd": return a_h,d_h - if arg[4]=="-r": + if arg[2]=="-r": a_h,d_h=read_data_homology(dir) assert(len(a_h)==len(d_h)) diff --git a/prepare_negative_dataset.py b/prepare_negative_dataset.py index 151ebec..2cf718e 100644 --- a/prepare_negative_dataset.py +++ b/prepare_negative_dataset.py @@ -5,6 +5,7 @@ import os import gc import progressbar import random +import pickle from read_data import read_data_homology,read_data_genome arg=sys.argv @@ -23,9 +24,12 @@ for df in progressbar.progressbar(a_h): for h in hgids: gmap[h]=1 -a=[] -d={} -a,d=read_data_genome("data",a,d) +data={} +with open("genome_maps","rb") as file: + data=pickle.load(file) +a=data["a"] +d=data["d"] +data={} assert(len(a)==len(d)) indexes_gd=[list(df.index.values) for df in a] d=list(d.keys()) diff --git a/prepare_other_features_negative.py b/prepare_other_features_negative.py new file mode 100644 index 0000000..6da9040 --- /dev/null +++ b/prepare_other_features_negative.py @@ -0,0 +1,55 @@ +import pandas as pd +import numpy as np +import json +import gc +import pickle +import progressbar +import sys +from tree_data import create_tree_data + + + +smg_name="negative_dataset_synteny_matrices_global" +sml_name="negative_dataset_synteny_matrices_local" +smi_name="negative_dataset_indexes" +dir_name="processed/synteny_matrices/" +label=2 + +rows=[] +df=pd.read_hdf("negative_dataset.h5",key="ndf") +try: + smg=np.load(dir_name+smg_name+".npy") + sml=np.load(dir_name+sml_name+".npy") + indexes=np.load(dir_name+smi_name+".npy") +except: + print("Incomplete data") + sys.exit(1) +df=df.loc[indexes] +for _,row in progressbar.progressbar(df.iterrows()): + row["homology_species"]=row["homology_species"].lower() +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 + 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("negative_dataset","wb") as file: + pickle.dump(rows,file) + +print("Data Saved Successfully:)") \ No newline at end of file diff --git a/prepare_synteny_matrix_negative.py b/prepare_synteny_matrix_negative.py index 909f58a..b448cd2 100644 --- a/prepare_synteny_matrix_negative.py +++ b/prepare_synteny_matrix_negative.py @@ -13,7 +13,6 @@ if not os.path.isdir("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)) @@ -22,21 +21,16 @@ 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)) +gene_sequences=read_gene_sequences(a_h,lsy,"geneseq","gene_sequences") 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) +synteny_matrices_global,synteny_matrices_local,indexes=synteny_matrix(gene_sequences,df,lsy,n,0) +np.save(ndir+"negative_dataset"+"_"+nf1,synteny_matrices_global) +np.save(ndir+"negative_dataset"+"_"+nf2,synteny_matrices_local) +np.save(ndir+"negative_dataset"+"_"+nf3,indexes) print("Synteny Matrices Created Successfully :)") diff --git a/read_data.py b/read_data.py index f599dd7..d97c238 100644 --- a/read_data.py +++ b/read_data.py @@ -35,8 +35,9 @@ def read_data_genome(dir_name,a,dict_ind_genome): data_gene[y]=data_gene[y].apply(clear_data) except: continue - #print(data_gene[0:10]) - data_gene=data_gene[data_gene['gene_biotype']=='protein_coding'] + #print(data_gene[0:10]) + data_gene=data_gene[(data_gene['gene_biotype']=='protein_coding') | (data_gene['gene_source']=='protein_coding')] + #print(data_gene[data_gene["gene_id"]=="ENSPMGG00000022088"]) a.append(data_gene) n=lf[x].split(".")[0] dict_ind_genome[n]=len(a)-1 diff --git a/update_neighbor_genes_ndf.py b/update_neighbor_genes_ndf.py index d3fd0c7..17f93f4 100644 --- a/update_neighbor_genes_ndf.py +++ b/update_neighbor_genes_ndf.py @@ -8,23 +8,6 @@ import pickle 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)!=3: - print("No. of arguments more or less. Please check") - sys.exit(1) - data={} with open("genome_maps","rb") as file: data=pickle.load(file) @@ -45,6 +28,6 @@ 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) +lsy=create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,cmap,cimap,save_after,0) print(len(lsy)) print("Neighbor Genes Updated Successfully") \ No newline at end of file