priyatham/compara-deep-learning
1
0
Fork 0
mirror of https://github.com/Priyatham-sai-chand/compara-deep-learning.git synced 2026-10-05 08:11:34 -07:00
Code Issues Projects Releases Packages Wiki Activity Actions

fix flank8 errors

Browse source
This commit is contained in:
Priyatham Sai Chand 2021-04-08 21:11:16 +05:30
parent 07f9f23101
commit 3f885458a4
17 changed files with 791 additions and 819 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

165
access_data_rest.py
View file

@ -3,152 +3,163 @@ import requests
import progressbar import progressbar
import sys import sys
def update_protein(gene_seq,gene):
t=0 def update_protein(gene_seq, gene):
while(t!=2): t = 0
while(t != 2):
try: try:
server = "https://rest.ensembl.org" server = "https://rest.ensembl.org"
ext = "/sequence/id/"+str(gene)+"?type=protein;multiple_sequences=1" ext = "/sequence/id/" + \
str(gene)+"?type=protein;multiple_sequences=1"
r = requests.get(server+ext, headers={ "Content-Type" : "application/json"}) r = requests.get(
server+ext, headers={"Content-Type": "application/json"})
if not r.ok: if not r.ok:
r.raise_for_status() r.raise_for_status()
sys.exit() sys.exit()
r=r.json() r = r.json()
if len(r)==1: if len(r) == 1:
r=dict(r[0]) r = dict(r[0])
gene_seq[gene]=str(r["seq"]) gene_seq[gene] = str(r["seq"])
return return
else: else:
maxi=0 maxi = 0
maxlen=0 maxlen = 0
for i in range(len(r)): for i in range(len(r)):
m=r[i] m = r[i]
m=dict(m) m = dict(m)
if len(m["seq"])>maxlen: if len(m["seq"]) > maxlen:
maxi=i maxi = i
r=dict(r[maxi]) r = dict(r[maxi])
gene_seq[gene]=str(r["seq"]) gene_seq[gene] = str(r["seq"])
return return
except : except Exception:
t+=1 t += 1
#print("\nError:",e) # print("\nError:",e)
continue continue
gene_seq[gene]="" gene_seq[gene] = ""
def update_rest_protein(data,fname):
gids={}
with open("processed/not_found_"+fname+".json","r") as file:
gids=dict(json.load(file))
gids=list(gids.keys()) def update_rest_protein(data, fname):
gids = {}
with open("processed/not_found_"+fname+".json", "r") as file:
gids = dict(json.load(file))
geneseq={} gids = list(gids.keys())
geneseq = {}
server = "https://rest.ensembl.org" server = "https://rest.ensembl.org"
ext = "/sequence/id?type=protein" ext = "/sequence/id?type=protein"
headers={ "Content-Type" : "application/json", "Accept" : "application/json"} headers = {"Content-Type": "application/json",
"Accept": "application/json"}
for i in progressbar.progressbar(range(0,len(gids)-50,50)): for i in progressbar.progressbar(range(0, len(gids)-50, 50)):
ids=dict(ids=list(gids[i:i+50])) ids = dict(ids=list(gids[i:i+50]))
while(1): while(1):
try: try:
r = requests.post(server+ext, headers=headers, data=str(json.dumps(ids))) r = requests.post(server+ext, headers=headers,
data=str(json.dumps(ids)))
if not r.ok: if not r.ok:
r.raise_for_status() r.raise_for_status()
gs=r.json() gs = r.json()
tgs={} tgs = {}
for g in gs: for g in gs:
tgs[g["query"]]=g["seq"] tgs[g["query"]] = g["seq"]
geneseq.update(tgs) geneseq.update(tgs)
break break
except Exception as e: except Exception as e:
print("Error:",e) print("Error:", e)
continue continue
data.update(geneseq) data.update(geneseq)
for genes in gids: for genes in gids:
try: try:
_=data[genes] _ = data[genes]
except: except Exception:
print(genes) print(genes)
update_protein(data,genes) update_protein(data, genes)
print("Protein Sequences Updated Successfully") print("Protein Sequences Updated Successfully")
return data return data
def update(gene_seq,gene):
t=0 def update(gene_seq, gene):
while(t!=2): t = 0
while(t != 2):
try: try:
server = "https://rest.ensembl.org" server = "https://rest.ensembl.org"
ext = "/sequence/id/"+str(gene)+"?type=cds;multiple_sequences=1" ext = "/sequence/id/"+str(gene)+"?type=cds;multiple_sequences=1"
r = requests.get(server+ext, headers={ "Content-Type" : "application/json"}) r = requests.get(
server+ext, headers={"Content-Type": "application/json"})
if not r.ok: if not r.ok:
r.raise_for_status() r.raise_for_status()
sys.exit() sys.exit()
r=r.json() r = r.json()
if len(r)==1: if len(r) == 1:
r=dict(r[0]) r = dict(r[0])
gene_seq[gene]=str(r["seq"]) gene_seq[gene] = str(r["seq"])
return return
else: else:
maxi=0 maxi = 0
maxlen=0 maxlen = 0
for i in range(len(r)): for i in range(len(r)):
m=r[i] m = r[i]
m=dict(m) m = dict(m)
if len(m["seq"])>maxlen: if len(m["seq"]) > maxlen:
maxi=i maxi = i
r=dict(r[maxi]) r = dict(r[maxi])
gene_seq[gene]=str(r["seq"]) gene_seq[gene] = str(r["seq"])
return return
except Exception as e: except Exception:
t+=1 t += 1
#print("\nError:",e) # print("\nError:",e)
continue continue
gene_seq[gene]="" gene_seq[gene] = ""
def update_rest(data,fname):
gids={}
with open("processed/not_found_"+fname+".json","r") as file:
gids=dict(json.load(file))
gids=list(gids.keys()) def update_rest(data, fname):
gids = {}
with open("processed/not_found_"+fname+".json", "r") as file:
gids = dict(json.load(file))
geneseq={} gids = list(gids.keys())
geneseq = {}
server = "https://rest.ensembl.org" server = "https://rest.ensembl.org"
ext = "/sequence/id?type=cds" ext = "/sequence/id?type=cds"
headers={ "Content-Type" : "application/json", "Accept" : "application/json"} headers = {"Content-Type": "application/json",
"Accept": "application/json"}
for i in progressbar.progressbar(range(0,len(gids)-50,50)): for i in progressbar.progressbar(range(0, len(gids)-50, 50)):
ids=dict(ids=list(gids[i:i+50])) ids = dict(ids=list(gids[i:i+50]))
while(1): while(1):
try: try:
r = requests.post(server+ext, headers=headers, data=str(json.dumps(ids))) r = requests.post(server+ext, headers=headers,
data=str(json.dumps(ids)))
if not r.ok: if not r.ok:
r.raise_for_status() r.raise_for_status()
gs=r.json() gs = r.json()
tgs={} tgs = {}
for g in gs: for g in gs:
tgs[g["query"]]=g["seq"] tgs[g["query"]] = g["seq"]
geneseq.update(tgs) geneseq.update(tgs)
break break
except Exception as e: except Exception as e:
print("Error:",e) print("Error:", e)
continue continue
data.update(geneseq) data.update(geneseq)
for genes in gids: for genes in gids:
try: try:
_=data[genes] _ = data[genes]
except: except Exception:
print(genes) print(genes)
update(data,genes) update(data, genes)
print("Gene Sequences Updated Successfully") print("Gene Sequences Updated Successfully")
return data return data

2
finalize_dataset.py
View file

@ -60,7 +60,7 @@ def prepare_features(a_h, d_h, sptree, label):
branch_length_species, \ branch_length_species, \
branch_length_homology_species, \ branch_length_homology_species, \
distance, dist_p_s, dist_p_hs = create_tree_data( distance, dist_p_s, dist_p_hs = create_tree_data(
sptree, df) sptree, df)
assert(len(branch_length_species) == len(df)) assert(len(branch_length_species) == len(df))
assert(len(sml) == len(distance)) assert(len(sml) == len(distance))

69
ftpg.py
View file

@ -3,30 +3,31 @@ import progressbar
import os import os
import sys import sys
import urllib.request as urllib import urllib.request as urllib
import requests
def download_data(x,dir_name):
fname=x.split("/")[-1]
path=os.path.join(dir_name,fname)
urllib.urlretrieve(x,path)
def get_data_file(file,dir): def download_data(x, dir_name):
fname = x.split("/")[-1]
path = os.path.join(dir_name, fname)
urllib.urlretrieve(x, path)
def get_data_file(file, dir):
if not os.path.isfile(file): if not os.path.isfile(file):
print("The specified file does not exist!!!") print("The specified file does not exist!!!")
sys.exit(1) sys.exit(1)
with open(file,"r")as f: with open(file, "r")as f:
lf=f.read().splitlines() lf = f.read().splitlines()
if not os.path.exists(dir): if not os.path.exists(dir):
os.mkdir(dir) os.mkdir(dir)
for x in progressbar.progressbar(lf): for x in progressbar.progressbar(lf):
download_data(x,dir) download_data(x, dir)
# This wil download all the fasta files for the coding sequences.
#This wil download all the fasta files for the coding sequences. To change the directory, change the argument in the get_data_file argument. # To change the directory, change the argument in the get_data_file argument.
host ="ftp.ensembl.org" host = "ftp.ensembl.org"
user = "anonymous" user = "anonymous"
password = "" password = ""
@ -34,50 +35,50 @@ print("Connecting to {}".format(host))
ftp = FTP(host) ftp = FTP(host)
ftp.login(user, password) ftp.login(user, password)
print("Connected to {}".format(host)) print("Connected to {}".format(host))
base_link="ftp://ftp.ensembl.org" base_link = "ftp://ftp.ensembl.org"
#find sequences of all the cds files # find sequences of all the cds files
l=ftp.nlst("/pub/release-96/fasta") list_of_files = ftp.nlst("/pub/release-96/fasta")
lt=[] lt = []
for x in l: for x in list_of_files:
y=ftp.nlst(x+"/cds") y = ftp.nlst(x+"/cds")
for z in y: for z in y:
if z.endswith(".cds.all.fa.gz"): if z.endswith(".cds.all.fa.gz"):
lt.append(z) lt.append(z)
with open("seq_link.txt","w") as file: with open("seq_link.txt", "w") as file:
for x in lt: for x in lt:
file.write(base_link+x) file.write(base_link+x)
file.write("\n") file.write("\n")
#find all the files with protein sequences # find all the files with protein sequences
l=ftp.nlst("/pub/release-96/fasta") list_of_files = ftp.nlst("/pub/release-96/fasta")
lt=[] lt = []
for x in progressbar.progressbar(l): for x in progressbar.progressbar(list_of_files):
y=ftp.nlst(x+"/pep") y = ftp.nlst(x+"/pep")
for z in y: for z in y:
if z.endswith(".pep.all.fa.gz"): if z.endswith(".pep.all.fa.gz"):
lt.append(z) lt.append(z)
with open("protein_seq.txt","w") as file: with open("protein_seq.txt", "w") as file:
for x in lt: for x in lt:
file.write(base_link+x) file.write(base_link+x)
file.write("\n") file.write("\n")
#get link of all the gtf files # get link of all the gtf files
l=ftp.nlst("/pub/release-96/gtf") list_of_files = ftp.nlst("/pub/release-96/gtf")
lt=[] lt = []
for x in l: for x in list_of_files:
y=ftp.nlst(x) y = ftp.nlst(x)
for z in y: for z in y:
if z.endswith(".96.gtf.gz"): if z.endswith(".96.gtf.gz"):
lt.append(z) lt.append(z)
with open("gtf_link.txt","w") as file: with open("gtf_link.txt", "w") as file:
for x in lt: for x in lt:
file.write(base_link+x) file.write(base_link+x)
file.write("\n") file.write("\n")
print("Downloading Data") print("Downloading Data")
get_data_file("gtf_link.txt","data") get_data_file("gtf_link.txt", "data")
get_data_file("seq_link.txt","geneseq") get_data_file("seq_link.txt", "geneseq")
get_data_file("protein_seq.txt","pro_seq") get_data_file("protein_seq.txt", "pro_seq")
print("Download Complete.................") print("Download Complete.................")

235
model.py
View file

@ -1,132 +1,161 @@
dim=10
n=3
fl=2*n+1
maxbl=29
import tensorflow as tf import tensorflow as tf
dim = 10
n = 3
fl = 2*n+1
maxbl = 29
def create_model(): def create_model():
tf.reset_default_graph() tf.reset_default_graph()
g=tf.Graph() g = tf.Graph()
with g.as_default(): with g.as_default():
synmg=tf.placeholder(dtype=tf.float64,shape=(None,2*n+1,2*n+1,2),name="Synteny_matrix_placeholder_Global") synmg = tf.placeholder(dtype=tf.float64, shape=(
synml=tf.placeholder(dtype=tf.float64,shape=(None,2*n+1,2*n+1,2),name="Synteny_matrix_placeholder_Local") None, 2*n+1, 2*n+1, 2), name="Synteny_matrix_placeholder_Global")
pfam=tf.placeholder(dtype=tf.float64,shape=(None,2*n+1,2*n+1,1),name="Pfam_matrix_placeholder") synml = tf.placeholder(dtype=tf.float64, shape=(
bls=tf.placeholder(dtype=tf.float64,shape=(None,maxbl),name="Species_Branch_Length_Placeholder") None, 2*n+1, 2*n+1, 2), name="Synteny_matrix_placeholder_Local")
blhs=tf.placeholder(dtype=tf.float64,shape=(None,maxbl),name="Homology_Species_Branch_Length_Placeholder") pfam = tf.placeholder(dtype=tf.float64, shape=(
#gl=tf.placeholder(dtype=tf.float64,shape=(None,1),name="Mean_gene_length") None, 2*n+1, 2*n+1, 1), name="Pfam_matrix_placeholder")
dps=tf.placeholder(dtype=tf.float64,shape=(None,1),name="mca_species_distance") bls = tf.placeholder(dtype=tf.float64, shape=(
dphs=tf.placeholder(dtype=tf.float64,shape=(None,1),name="mca_homology_species_distance") None, maxbl), name="Species_Branch_Length_Placeholder")
dis=tf.placeholder(dtype=tf.float64,shape=(None,1),name="total_distance") blhs = tf.placeholder(dtype=tf.float64, shape=(
lr=tf.placeholder(dtype=tf.float64,shape=(),name="learning_rate") None, maxbl), name="Homology_Species_Branch_Length_Placeholder")
y=tf.placeholder(dtype=tf.int32,shape=(None),name="labels") # gl=tf.placeholder(dtype=tf.float64,shape=(None,1),name="Mean_gene_length")
dps = tf.placeholder(dtype=tf.float64, shape=(
lrs=tf.summary.scalar("Learning_Rate",lr) None, 1), name="mca_species_distance")
dphs = tf.placeholder(dtype=tf.float64, shape=(
x=tf.concat([dps,dps-dphs,dis],1,name="Create_train_vector") None, 1), name="mca_homology_species_distance")
dis = tf.placeholder(dtype=tf.float64, shape=(
None, 1), name="total_distance")
lr = tf.placeholder(dtype=tf.float64, shape=(), name="learning_rate")
y = tf.placeholder(dtype=tf.int32, shape=(None), name="labels")
print(synmg,"\n",synml,"\n",bls,"\n",blhs,"\n",dps,"\n",dphs,"\n",dis,"\n",x) # lrs = tf.summary.scalar("Learning_Rate", lr)
reg_l2=tf.contrib.layers.l2_regularizer(0.001) x = tf.concat([dps, dps-dphs, dis], 1, name="Create_train_vector")
reg_l1 = tf.contrib.layers.l1_regularizer(scale=0.005, scope=None)
def get_variable_by_shape(shape,name): print(synmg, "\n", synml, "\n", bls, "\n", blhs,
f=tf.get_variable(name,shape=shape,initializer=tf.glorot_uniform_initializer(),dtype=tf.float64,regularizer=reg_l2) "\n", dps, "\n", dphs, "\n", dis, "\n", x)
reg_l2 = tf.contrib.layers.l2_regularizer(0.001)
# reg_l1 = tf.contrib.layers.l1_regularizer(scale=0.005, scope=None)
def get_variable_by_shape(shape, name):
f = tf.get_variable(name, shape=shape,
initializer=tf.glorot_uniform_initializer(),
dtype=tf.float64,
regularizer=reg_l2)
return f return f
def create_synteny_aligner(name,synm): def create_synteny_aligner(name, synm):
with tf.variable_scope(name+"Synteny_Aligner",reuse=tf.AUTO_REUSE): with tf.variable_scope(
name+"Synteny_Aligner", reuse=tf.AUTO_REUSE):
fconv=get_variable_by_shape((2,2,2,dim),"fconv") fconv = get_variable_by_shape((2, 2, 2, dim), "fconv")
conv=tf.nn.conv2d(synm,fconv,(1,1,1,1),padding="VALID",name="Conv_aligner") conv = tf.nn.conv2d(synm, fconv, (1, 1, 1, 1),
padding="VALID", name="Conv_aligner")
fconv_1=get_variable_by_shape((2,2,dim,dim*2),"fconv_1") fconv_1 = get_variable_by_shape((2, 2, dim, dim*2), "fconv_1")
conv_1=tf.nn.conv2d(conv,fconv_1,(1,1,1,1),padding="VALID",name="Conv_aligner_1") conv_1 = tf.nn.conv2d(
conv, fconv_1, (1, 1, 1, 1), padding="VALID",
fxconv=get_variable_by_shape((fl,2,dim*2),"fxconv") name="Conv_aligner_1")
x_conv=tf.reshape(synm,(-1,fl*fl,2))
x_conv=tf.nn.conv1d(x_conv,fxconv,stride=fl,padding="SAME",name="row_aligner") fxconv = get_variable_by_shape((fl, 2, dim*2), "fxconv")
x_conv = tf.reshape(synm, (-1, fl*fl, 2))
y_conv=tf.reshape(tf.transpose(synm,(0,2,1,3)),(-1,fl*fl,2)) x_conv = tf.nn.conv1d(
fyconv=get_variable_by_shape((fl,2,dim*2),"fyconv") x_conv, fxconv, stride=fl, padding="SAME",
y_conv=tf.nn.conv1d(y_conv,fyconv,stride=fl,padding="SAME",name="column_aligner") name="row_aligner")
y_conv = tf.reshape(tf.transpose(
wconv=get_variable_by_shape((fl,fl,2,dim*2*10),"wconv") synm, (0, 2, 1, 3)), (-1, fl*fl, 2))
w_conv=tf.nn.conv2d(synm,wconv,(1,1,1,1),padding="VALID",name="Global_Aligner_1") fyconv = get_variable_by_shape((fl, 2, dim*2), "fyconv")
y_conv = tf.nn.conv1d(
conv_1=tf.reshape(conv_1,(-1,25,dim*2)) y_conv, fyconv, stride=fl, padding="SAME",
x_conv=tf.reshape(x_conv,(-1,fl,dim*2)) name="column_aligner")
y_conv=tf.reshape(y_conv,(-1,fl,dim*2))
w_conv=tf.reshape(w_conv,(-1,10,dim*2)) wconv = get_variable_by_shape((fl, fl, 2, dim*2*10), "wconv")
conv_final=tf.concat([conv_1,x_conv,y_conv,w_conv],1,name="Concatenate_All_Alignments") w_conv = tf.nn.conv2d(
synm, wconv, (1, 1, 1, 1), padding="VALID",
name="Global_Aligner_1")
conv_1 = tf.reshape(conv_1, (-1, 25, dim*2))
x_conv = tf.reshape(x_conv, (-1, fl, dim*2))
y_conv = tf.reshape(y_conv, (-1, fl, dim*2))
w_conv = tf.reshape(w_conv, (-1, 10, dim*2))
conv_final = tf.concat(
[conv_1, x_conv, y_conv, w_conv], 1,
name="Concatenate_All_Alignments")
return conv_final return conv_final
with tf.variable_scope("Pfam",reuse=tf.AUTO_REUSE): with tf.variable_scope("Pfam", reuse=tf.AUTO_REUSE):
wconv_pfam=get_variable_by_shape((fl,fl,1,dim*2*10),"wconv_pfam") wconv_pfam = get_variable_by_shape(
w_conv_pfam=tf.nn.conv2d(pfam,wconv_pfam,(1,1,1,1),padding="VALID",name="Global_Aligner_pfam") (fl, fl, 1, dim*2*10), "wconv_pfam")
w_conv_pfam=tf.reshape(w_conv_pfam,(-1,10,dim*2)) w_conv_pfam = tf.nn.conv2d(
pfam, wconv_pfam, (1, 1, 1, 1), padding="VALID",
name="Global_Aligner_pfam")
w_conv_pfam = tf.reshape(w_conv_pfam, (-1, 10, dim*2))
conv_final_g=create_synteny_aligner("Global_",synmg) conv_final_g = create_synteny_aligner("Global_", synmg)
conv_final_l=create_synteny_aligner("Local_",synml) conv_final_l = create_synteny_aligner("Local_", synml)
final=tf.concat([conv_final_g,conv_final_l,w_conv_pfam],1) final = tf.concat([conv_final_g, conv_final_l, w_conv_pfam], 1)
#final=conv_final_l # final=conv_final_l
with tf.variable_scope("Combine_Renormalize",reuse=tf.AUTO_REUSE): with tf.variable_scope("Combine_Renormalize", reuse=tf.AUTO_REUSE):
bl=tf.concat([bls,blhs],1) bl = tf.concat([bls, blhs], 1)
#bl=bls-blhs # bl=bls-blhs
theta_bl=get_variable_by_shape((maxbl*2,1),"theta_bl") theta_bl = get_variable_by_shape((maxbl*2, 1), "theta_bl")
theta_bl=tf.matmul(bl,theta_bl) theta_bl = tf.matmul(bl, theta_bl)
x=tf.concat([x,theta_bl],1) x = tf.concat([x, theta_bl], 1)
theta=get_variable_by_shape((4,108),"theta") theta = get_variable_by_shape((4, 108), "theta")
bias=get_variable_by_shape((1,108),"b") bias = get_variable_by_shape((1, 108), "b")
theta_2=tf.matmul(x,theta)+bias theta_2 = tf.matmul(x, theta)+bias
theta_2=tf.reshape(theta_2,(-1,108,1)) theta_2 = tf.reshape(theta_2, (-1, 108, 1))
theta_2=tf.tile(theta_2,[1,1,dim*2]) theta_2 = tf.tile(theta_2, [1, 1, dim*2])
final=final*theta_2 final = final*theta_2
print(final) print(final)
flat=tf.layers.flatten(final) flat = tf.layers.flatten(final)
zero=tf.constant(0.0,dtype=tf.float64) zero = tf.constant(0.0, dtype=tf.float64)
diff=dps-dphs diff = dps-dphs
print(diff) print(diff)
diff_2=tf.cast(tf.equal(diff,zero),tf.float64) diff_2 = tf.cast(tf.equal(diff, zero), tf.float64)
print(diff_2) print(diff_2)
diff=tf.tile(diff_2,[1,dim*10]) diff = tf.tile(diff_2, [1, dim*10])
flat=tf.concat([flat,diff],1) flat = tf.concat([flat, diff], 1)
print(flat) print(flat)
#dense=tf.layers.dense(flat,2048,kernel_regularizer=reg_l2,bias_regularizer=reg_l2) # dense=tf.layers.dense(flat,2048,kernel_regularizer=reg_l2,bias_regularizer=reg_l2)
#dense_2=tf.layers.dense(dense,1024,kernel_regularizer=reg_l2,bias_regularizer=reg_l2) # dense_2=tf.layers.dense(dense,1024,kernel_regularizer=reg_l2,bias_regularizer=reg_l2)
dense_3=tf.layers.dense(flat,512,kernel_regularizer=reg_l2,bias_regularizer=reg_l2) dense_3 = tf.layers.dense(
flat, 512, kernel_regularizer=reg_l2, bias_regularizer=reg_l2)
logits_pred=tf.layers.dense(dense_3,3,name="Predictions")
logits_pred = tf.layers.dense(dense_3, 3, name="Predictions")
print(logits_pred) print(logits_pred)
entropy=tf.nn.sparse_softmax_cross_entropy_with_logits(logits=logits_pred,labels=y) entropy = tf.nn.sparse_softmax_cross_entropy_with_logits(
logits=logits_pred, labels=y)
print(entropy) print(entropy)
#weights = tf.trainable_variables() # all vars of your graph # weights = tf.trainable_variables() # all vars of your graph
#regl1 = tf.contrib.layers.apply_regularization(reg_l1, weights) # regl1 = tf.contrib.layers.apply_regularization(reg_l1, weights)
reg_losses = tf.get_collection(tf.GraphKeys.REGULARIZATION_LOSSES) reg_losses = tf.get_collection(tf.GraphKeys.REGULARIZATION_LOSSES)
reg_constant =0.00000001 reg_constant = 0.00000001
loss=tf.reduce_mean(entropy)+reg_constant * sum(reg_losses) loss = tf.reduce_mean(entropy)+reg_constant * sum(reg_losses)
#loss=tf.reduce_mean(entropy) # loss=tf.reduce_mean(entropy)
optimizer=tf.train.RMSPropOptimizer(lr) optimizer = tf.train.RMSPropOptimizer(lr)
#optimizer=tf.train.AdamOptimizer() # optimizer=tf.train.AdamOptimizer()
losses=tf.summary.scalar("Loss",loss) # losses = tf.summary.scalar("Loss", loss)
t_op=optimizer.minimize(loss) t_op = optimizer.minimize(loss)
acc=tf.math.in_top_k(tf.cast(logits_pred,tf.float32),y,1) acc = tf.math.in_top_k(tf.cast(logits_pred, tf.float32), y, 1)
accuracy=tf.reduce_mean(tf.cast(acc,tf.float32)) accuracy = tf.reduce_mean(tf.cast(acc, tf.float32))
accs=tf.summary.scalar("Accuracy",accuracy) # accs = tf.summary.scalar("Accuracy", accuracy)
summary=tf.summary.merge_all() summary = tf.summary.merge_all()
init=tf.global_variables_initializer() init = tf.global_variables_initializer()
saver=tf.train.Saver() saver = tf.train.Saver()
for node in (synmg,synml,pfam,bls,blhs,dps,dphs,dis,lr,y): for node in (synmg, synml, pfam, bls, blhs, dps, dphs, dis, lr, y):
g.add_to_collection("input_nodes",node) g.add_to_collection("input_nodes", node)
for node in (loss,t_op,accuracy,init,summary): for node in (loss, t_op, accuracy, init, summary):
g.add_to_collection("output_nodes",node) g.add_to_collection("output_nodes", node)
return g,saver return g, saver

86
pfam_folder_pred.py
View file

@ -1,62 +1,66 @@
import json import pandas as pd
import gc
import pandas as pd
import numpy as np
import pickle
import sys import sys
import progressbar import progressbar
import os import os
from neighbor_genes import read_genome_maps from neighbor_genes import read_genome_maps
from process_data import create_data_homology_ls from process_data import create_data_homology_ls
from read_get_gene_seq import read_gene_sequences from read_get_gene_seq import read_gene_sequences
from access_data_rest import update_rest,update_rest_protein from access_data_rest import update_rest_protein
from prepare_synteny_matrix import write_fasta from prepare_synteny_matrix import write_fasta
from process_data import create_map_list
def read_database(fname,dirname):
df=pd.read_csv(dirname+"/"+fname,sep="\t",header=None) def read_database(fname, dirname):
label_dict=dict(ortholog_one2one=1, df = pd.read_csv(dirname+"/"+fname, sep="\t", header=None)
other_paralog=0, label_dict = dict(ortholog_one2one=1,
non_homolog=2, other_paralog=0,
ortholog_one2many=1, non_homolog=2,
ortholog_many2many=1, ortholog_one2many=1,
within_species_paralog=0, ortholog_many2many=1,
gene_split=4) within_species_paralog=0,
label=[] gene_split=4)
for _,row in df.iterrows(): label = []
for _, row in df.iterrows():
label.append(label_dict[row[7]]) label.append(label_dict[row[7]])
df=df.assign(label=label) df = df.assign(label=label)
df=df.drop(7,axis=1) df = df.drop(7, axis=1)
df=df.drop(0,axis=1) df = df.drop(0, axis=1)
df.columns=["gene_stable_id","species","homology_gene_stable_id","homology_species","goc","wga","label"] df.columns = ["gene_stable_id", "species", "homology_gene_stable_id",
"homology_species", "goc", "wga", "label"]
return df return df
def read_prediction_file_folder(dir_name): def read_prediction_file_folder(dir_name):
lf=os.listdir(dir_name) lf = os.listdir(dir_name)
a_h=[] a_h = []
d_h=[] d_h = []
for x in progressbar.progressbar(lf): for x in progressbar.progressbar(lf):
df=read_database(x,dir_name) df = read_database(x, dir_name)
a_h.append(df) a_h.append(df)
d_h.append(x.split(".")[0]) d_h.append(x.split(".")[0])
return a_h,d_h return a_h, d_h
def create_synteny_features(a_h,d_h,n,a,d,ld,ldg,cmap,cimap,name):
lsy=create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,cmap,cimap,0) def create_synteny_features(a_h, d_h, n, a, d, ld, ldg, cmap, cimap, name):
protein_sequences=read_gene_sequences(a_h,lsy,"pro_seq","prediction_"+name) lsy = create_data_homology_ls(a_h, d_h, n, a, d, ld, ldg, cmap, cimap, 0)
protein_sequences=update_rest_protein(protein_sequences,"prediction_"+name) protein_sequences = read_gene_sequences(
write_fasta(protein_sequences,"prediction_"+name) a_h, lsy, "pro_seq", "prediction_"+name)
protein_sequences = update_rest_protein(
protein_sequences, "prediction_"+name)
write_fasta(protein_sequences, "prediction_"+name)
print("Protein Sequences Loaded") print("Protein Sequences Loaded")
def main():
arg=sys.argv
dirname=arg[-1]
a_h,d_h=read_prediction_file_folder(dirname)
n=3
a,d,ld,ldg,cmap,cimap=read_genome_maps()#read the genome maps def main():
arg = sys.argv
dirname = arg[-1]
a_h, d_h = read_prediction_file_folder(dirname)
n = 3
a, d, ld, ldg, cmap, cimap = read_genome_maps() # read the genome maps
print("Genome Maps Loaded.") print("Genome Maps Loaded.")
create_synteny_features(a_h,d_h,n,a,d,ld,ldg,cmap,cimap,dirname) create_synteny_features(a_h, d_h, n, a, d, ld, ldg, cmap, cimap, dirname)
if __name__=="__main__":
main()
if __name__ == "__main__":
main()

178
pfam_matrix.py
View file

@ -1,137 +1,151 @@
import pandas as pd import pandas as pd
import numpy as np import numpy as np
import progressbar import progressbar
import os import os
import json import json
import sys import sys
from prepare_synteny_matrix import read_data_homology,load_neighbor_genes from prepare_synteny_matrix import read_data_homology, load_neighbor_genes
from process_data import create_map_list from process_data import create_map_list
from process_negative import read_database_txt from process_negative import read_database_txt
def get_score_overlap(x,y,pfam_db,pfam_map): def get_score_overlap(x, y, pfam_db, pfam_map):
df_1=pfam_db.loc[pfam_map[x]] df_1 = pfam_db.loc[pfam_map[x]]
df_2=pfam_db.loc[pfam_map[y]] df_2 = pfam_db.loc[pfam_map[y]]
l=list(df_2.domain) list_of_domains = list(df_2.domain)
c=0 c = 0
c_1=0 c_1 = 0
for _,row in df_1.iterrows(): for _, row in df_1.iterrows():
if row.domain in l:#check if the domain exists in the list # check if the domain exists in the list
c_1+=1 if row.domain in list_of_domains:
st=int(df_2[df_2["domain"]==row.domain].hmm_from)#get the start c_1 += 1
end=int(df_2[df_2["domain"]==row.domain].hmm_to)#get the end # get the start
if (int(row.hmm_from)>st and int(row.hmm_from)<end) or (int(row.hmm_to)>st and int(row.hmm_from)<end):#check if the domain is a ovelapping domain st = int(df_2[df_2["domain"] == row.domain].hmm_from)
c+=1 # get the end
return c/max(len(df_1),len(df_2)) end = int(df_2[df_2["domain"] == row.domain].hmm_to)
# check if the domain is a ovelapping domain
if (int(
row.hmm_from) > st and int(row.hmm_from) < end) or (int(
row.hmm_to) > st and int(row.hmm_from) < end):
c += 1
return c/max(len(df_1), len(df_2))
def pfam_matrix(g1,g2,n,pfam_db,gmap,pfam_map):
pm=np.zeros((n,n)) def pfam_matrix(g1, g2, n, pfam_db, gmap, pfam_map):
pm = np.zeros((n, n))
for i in range(n): for i in range(n):
if g1[i]=="NULL_GENE": if g1[i] == "NULL_GENE":
continue continue
try: try:
_=gmap[g1[i]] _ = gmap[g1[i]]
except: except Exception:
continue continue
for j in range(n): for j in range(n):
if g2[j]=="NULL_GENE": if g2[j] == "NULL_GENE":
continue continue
try: try:
_=gmap[g2[j]] _ = gmap[g2[j]]
except: except Exception:
continue continue
pm[i][j]=get_score_overlap(g1[i],g2[j],pfam_db,pfam_map) pm[i][j] = get_score_overlap(g1[i], g2[j], pfam_db, pfam_map)
return pm return pm
def create_pfam_matrix(df,lsy,pfam_db,pfam_map):
n=3 def create_pfam_matrix(df, lsy, pfam_db, pfam_map):
glist=list(pfam_db.gene_stable_id) n = 3
gmap=create_map_list(glist) glist = list(pfam_db.gene_stable_id)
pg=[] gmap = create_map_list(glist)
indexes=[] pg = []
for index,row in progressbar.progressbar(df.iterrows()): indexes = []
g1=str(row["gene_stable_id"]) for index, row in progressbar.progressbar(df.iterrows()):
g2=str(row["homology_gene_stable_id"]) g1 = str(row["gene_stable_id"])
x=[] g2 = str(row["homology_gene_stable_id"])
y=[] x = []
y = []
try: try:
_=lsy[g1] _ = lsy[g1]
_=lsy[g2] _ = lsy[g2]
except: except Exception:
continue continue
for i in range(len(lsy[g1]['b'])-1,-1,-1): for i in range(len(lsy[g1]['b'])-1, -1, -1):
x.append(lsy[g1]['b'][i]) x.append(lsy[g1]['b'][i])
x.append(g1) x.append(g1)
for k in lsy[g1]['f']: for k in lsy[g1]['f']:
x.append(k) x.append(k)
for i in range(len(lsy[g2]['b'])-1,-1,-1): for i in range(len(lsy[g2]['b'])-1, -1, -1):
y.append(lsy[g2]['b'][i]) y.append(lsy[g2]['b'][i])
y.append(g2) y.append(g2)
for k in lsy[g2]['f']: for k in lsy[g2]['f']:
y.append(k) y.append(k)
assert(len(x)==len(y)) assert(len(x) == len(y))
assert(len(x)==(2*n+1)) assert(len(x) == (2*n+1))
pmtemp=pfam_matrix(x,y,2*n+1,pfam_db,gmap,pfam_map) pmtemp = pfam_matrix(x, y, 2*n+1, pfam_db, gmap, pfam_map)
pg.append(pmtemp) pg.append(pmtemp)
indexes.append(index) indexes.append(index)
return np.array(pg),np.array(indexes) return np.array(pg), np.array(indexes)
def create_pfam_map(pfam_db): def create_pfam_map(pfam_db):
pfam_map={} pfam_map = {}
for index,row in progressbar.progressbar(pfam_db.iterrows()): for index, row in progressbar.progressbar(pfam_db.iterrows()):
try: try:
_=pfam_map[row.gene_stable_id] _ = pfam_map[row.gene_stable_id]
except: except Exception:
pfam_map[row.gene_stable_id]=[] pfam_map[row.gene_stable_id] = []
pfam_map[row.gene_stable_id].append(index) pfam_map[row.gene_stable_id].append(index)
return pfam_map return pfam_map
def main_positive(): def main_positive():
if not os.path.isdir("processed/pfam_matrices"): if not os.path.isdir("processed/pfam_matrices"):
os.mkdir("processed/pfam_matrices") os.mkdir("processed/pfam_matrices")
a_h,d_h=read_data_homology("data_homology") a_h, d_h = read_data_homology("data_homology")
lsy=load_neighbor_genes() lsy = load_neighbor_genes()
pfam_db=pd.read_hdf("pfam_db_positive.h5") pfam_db = pd.read_hdf("pfam_db_positive.h5")
pfam_map=create_pfam_map(pfam_db) pfam_map = create_pfam_map(pfam_db)
ndir="processed/pfam_matrices/" ndir = "processed/pfam_matrices/"
nf1="pfam_matrices" nf1 = "pfam_matrices"
nf3="pfam_indexes" nf3 = "pfam_indexes"
for i in range(len(a_h)): for i in range(len(a_h)):
df=a_h[i] df = a_h[i]
print(len(df)) print(len(df))
pfam_matrices,indexes=create_pfam_matrix(df,lsy,pfam_db,pfam_map) pfam_matrices, indexes = create_pfam_matrix(df, lsy, pfam_db, pfam_map)
np.save(ndir+str(d_h[i])+"_"+nf1,pfam_matrices) np.save(ndir+str(d_h[i])+"_"+nf1, pfam_matrices)
np.save(ndir+str(d_h[i])+"_"+nf3,indexes) np.save(ndir+str(d_h[i])+"_"+nf3, indexes)
print(len(indexes)) print(len(indexes))
def read_data_negative(arg): def read_data_negative(arg):
df=read_database_txt(arg[-1]) df = read_database_txt(arg[-1])
name=arg[-1].split(".")[0] name = arg[-1].split(".")[0]
ind=np.load("processed/synteny_matrices/"+name+"_indexes.npy") ind = np.load("processed/synteny_matrices/"+name+"_indexes.npy")
df=df.loc[ind] df = df.loc[ind]
pfam_db=pd.read_hdf("pfam_db_negative.h5") pfam_db = pd.read_hdf("pfam_db_negative.h5")
pfam_map=create_pfam_map(pfam_db) pfam_map = create_pfam_map(pfam_db)
with open("processed/neighbor_genes_negative.json","r") as file: with open("processed/neighbor_genes_negative.json", "r") as file:
lsy=dict(json.load(file)) lsy = dict(json.load(file))
return df,pfam_db,pfam_map,lsy,name return df, pfam_db, pfam_map, lsy, name
def main_negative(arg): def main_negative(arg):
df,pfam_db,pfam_map,lsy,name=read_data_negative(arg) df, pfam_db, pfam_map, lsy, name = read_data_negative(arg)
ndir="processed/pfam_matrices/" ndir = "processed/pfam_matrices/"
nf1="pfam_matrices" nf1 = "pfam_matrices"
nf3="pfam_indexes" nf3 = "pfam_indexes"
print(len(df)) print(len(df))
pfam_matrices,indexes=create_pfam_matrix(df,lsy,pfam_db,pfam_map) pfam_matrices, indexes = create_pfam_matrix(df, lsy, pfam_db, pfam_map)
np.save(ndir+name+"_"+nf1,pfam_matrices) np.save(ndir+name+"_"+nf1, pfam_matrices)
np.save(ndir+name+"_"+nf3,indexes) np.save(ndir+name+"_"+nf3, indexes)
print(len(indexes)) print(len(indexes))
def main(): def main():
arg=sys.argv arg = sys.argv
main_positive() main_positive()
main_negative(arg) main_negative(arg)
if __name__=="__main__":
if __name__ == "__main__":
main() main()

55
pfam_parser.py
View file

@ -1,50 +1,53 @@
import pandas as pd import pandas as pd
import progressbar import progressbar
import gc import gc
import sys import sys
def pfam_parse(filename): def pfam_parse(filename):
rlist=[] rlist = []
try: try:
with open(filename) as file: with open(filename) as file:
file.seek(0,0) file.seek(0, 0)
for line in progressbar.progressbar(file): for line in progressbar.progressbar(file):
if line.startswith("#"): if line.startswith("#"):
continue continue
temp_dict={} temp_dict = {}
x=[y for y in line.split(" ") if y !=''] x = [y for y in line.split(" ") if y != '']
if x[9] !="1": if x[9] != "1":
continue continue
temp_dict["gene_stable_id"]=x[3] temp_dict["gene_stable_id"] = x[3]
temp_dict["accession"]=x[1] temp_dict["accession"] = x[1]
temp_dict["tlen"]=x[2] temp_dict["tlen"] = x[2]
temp_dict["qlen"]=x[5] temp_dict["qlen"] = x[5]
temp_dict["domain"]=x[0] temp_dict["domain"] = x[0]
temp_dict["hmm_from"]=x[15] temp_dict["hmm_from"] = x[15]
temp_dict["hmm_to"]=x[16] temp_dict["hmm_to"] = x[16]
temp_dict["ali_from"]=x[17] temp_dict["ali_from"] = x[17]
temp_dict["ali_to"]=x[18] temp_dict["ali_to"] = x[18]
temp_dict["env_from"]=x[19] temp_dict["env_from"] = x[19]
temp_dict["env_to"]=x[20] temp_dict["env_to"] = x[20]
rlist.append(temp_dict) rlist.append(temp_dict)
except Exception as e: except Exception as e:
print(e) print(e)
return pd.DataFrame() return pd.DataFrame()
print(len(rlist)) print(len(rlist))
tdf=pd.DataFrame(rlist) tdf = pd.DataFrame(rlist)
return tdf return tdf
def main(): def main():
arg=sys.argv arg = sys.argv
fname_1=arg[-2] fname_1 = arg[-2]
fname_2=arg[-1] fname_2 = arg[-1]
df=pfam_parse(fname_1) df = pfam_parse(fname_1)
df.to_hdf("pfam_db_positive.h5",key="pfam_db_positive",mode="w") df.to_hdf("pfam_db_positive.h5", key="pfam_db_positive", mode="w")
df="" df = ""
gc.collect() gc.collect()
df=pfam_parse(fname_2) df = pfam_parse(fname_2)
df.to_hdf("pfam_db_negative.h5",key="pfam_db_negative",mode="w") df.to_hdf("pfam_db_negative.h5", key="pfam_db_negative", mode="w")
print("Pfam Databases Written Successfully :)") print("Pfam Databases Written Successfully :)")
if __name__ == "__main__": if __name__ == "__main__":
main() main()

3
prediction.py
View file

@ -130,7 +130,8 @@ def write_preds(fname, model_name, name, preds, index_dict, df):
"_" + "_" +
name + name +
"_multiple.txt") "_multiple.txt")
with open("prediction_" + fname + "_" + model_name + "_" + name + "_multiple.txt", "w") as file: with open("prediction_" + fname + "_" + model_name + "_" + name +
"_multiple.txt", "w") as file:
for index, row in progressbar.progressbar(df.iterrows()): for index, row in progressbar.progressbar(df.iterrows()):
file.write(str(row[0])) file.write(str(row[0]))
file.write("\t") file.write("\t")

245
prediction_pfam.py
View file

@ -1,211 +1,54 @@
import json import pandas as pd
import gc
import pandas as pd
import numpy as np
import pickle
import tensorflow as tf
import sys
import os
import progressbar import progressbar
from neighbor_genes import read_genome_maps import gc
from process_data import create_data_homology_ls import sys
from read_get_gene_seq import read_gene_sequences
from access_data_rest import update_rest,update_rest_protein
from threads import Procerssrunner
from prepare_synteny_matrix import read_data_synteny
from tree_data import create_tree_data
from process_data import create_map_list
from pfam_parser import pfam_parse
from pfam_matrix import create_pfam_map,create_pfam_matrix
def read_database(fname):
df=pd.read_csv(fname,sep="\t",header=None)
label_dict=dict(ortholog_one2one=1,
other_paralog=0,
non_homolog=2,
ortholog_one2many=1,
ortholog_many2many=1,
within_species_paralog=0,
gene_split=4)
label=[]
for _,row in df.iterrows():
label.append(label_dict[row[7]])
df=df.assign(label=label)
df=df.drop(7,axis=1)
df=df.drop(0,axis=1)
df.columns=["gene_stable_id","species","homology_gene_stable_id","homology_species","goc","wga","label"]
return df
def select_data_by_length(df,st,end): def pfam_parse(filename):
rlist = []
try: try:
if end<len(df): with open(filename) as file:
if st<end: file.seek(0, 0)
df=df.loc[df.index.values[st:end]] for line in progressbar.progressbar(file):
else: if line.startswith("#"):
raise ValueError() continue
except: temp_dict = {}
print("Making Predictions for the complete dataframe:)") x = [y for y in line.split(" ") if y != '']
print(len(df)) if x[9] != "1":
return df continue
temp_dict["gene_stable_id"] = x[3]
temp_dict["accession"] = x[1]
temp_dict["tlen"] = x[2]
temp_dict["qlen"] = x[5]
temp_dict["domain"] = x[0]
temp_dict["hmm_from"] = x[15]
temp_dict["hmm_to"] = x[16]
temp_dict["ali_from"] = x[17]
temp_dict["ali_to"] = x[18]
temp_dict["env_from"] = x[19]
temp_dict["env_to"] = x[20]
rlist.append(temp_dict)
except Exception as e:
print(e)
return pd.DataFrame()
print(len(rlist))
tdf = pd.DataFrame(rlist)
return tdf
def create_synteny_features(a_h,d_h,n,a,d,ld,ldg,cmap,cimap,name):
lsy=create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,cmap,cimap,0)
gene_sequences=read_gene_sequences(a_h,lsy,"geneseq","prediction_"+name)
gene_sequences=update_rest(gene_sequences,"prediction_"+name)
print("Gene Sequences Loaded.")
return lsy,gene_sequences
def threadmaker(nop,df,lsy,gene_sequences,n,name):
part=len(df)//nop
pr=Procerssrunner()
pr.start_processes(nop,df,gene_sequences,lsy,part,n,name)
smg,sml,indexes=read_data_synteny(nop,name)
sml=np.array(sml)
smg=np.array(smg)
indexes=np.array(indexes)
return sml,smg,indexes
def get_prediction(smg,sml,pfam_matrices,indexes,bls,blhs,dis,dps,dphs,model_name,no_of_model,w):
preds=np.zeros((len(smg),no_of_model))
pfam_matrices=pfam_matrices.reshape((len(smg),7,7,1))
for i in range(1,no_of_model+1):
try:
model=tf.train.import_meta_graph(model_name+'_v'+str(i)+'/model.ckpt.meta')
except:
print("Something wrong with the model.")
continue
with tf.Session() as sess:
try:
model.restore(sess,model_name+'_v'+str(i)+"/model.ckpt")
graph = tf.get_default_graph()
synmgt,synmlt,pfamt,blst,blhst,dpst,dphst,dist,lrt,yt=graph.get_collection("input_nodes")
predictions=graph.get_tensor_by_name("Predictions/BiasAdd:0")
print("Model Loaded Successfully :)")
except:
print(":(")
sys.exit()
fd={synmgt:smg,
synmlt:sml,
pfamt:pfam_matrices,
blst:bls,
blhst:blhs,
dpst:dps.reshape((len(blhs),1)),
dist:dis.reshape((len(blhs),1)),
dphst:dphs.reshape((len(blhs),1))}
preds_t_1=sess.run([predictions],feed_dict=fd)
preds_t_1=np.array(preds_t_1)[0]
fd={synmgt:smg.transpose((0,2,1,3)),
synmlt:sml.transpose((0,2,1,3)),
pfamt:pfam_matrices.transpose((0,2,1,3)),
blst:blhs,
blhst:bls,
dpst:dphs.reshape((len(blhs),1)),
dist:dis.reshape((len(blhs),1)),
dphst:dps.reshape((len(blhs),1))}
preds_t_2=sess.run([predictions],feed_dict=fd)
preds_t_2=np.array(preds_t_2)[0]
preds=preds+w[i-1]*(preds_t_1+preds_t_2)/2
tf.reset_default_graph()
preds=np.argmax(preds,axis=1)
print(preds.shape)
return preds
def write_preds(fname,model_name,name,preds,index_dict,df):
print("Writing predcitions to:","prediction_"+fname+"_"+model_name+"_"+name+"_multiple_pfam.txt")
with open("prediction_"+fname+"_"+model_name+"_"+name+"_multiple_pfam.txt","w") as file:
for index,row in progressbar.progressbar(df.iterrows()):
file.write(row["gene_stable_id"])
file.write("\t")
file.write(row["homology_gene_stable_id"])
file.write("\t")
file.write(str(row["label"]))
file.write("\t")
if index in index_dict:
file.write(str(preds[index_dict[index]]))
file.write("\t")
if preds[index_dict[index]]==row["label"]:
file.write(str(1))
else:
file.write(str(0))
else:
file.write("Error")
file.write("\t")
file.write("NaN")
file.write("\n")
def main(): def main():
arg=sys.argv arg = sys.argv
arg=arg[1:] fname_1 = arg[-2]
fname=arg[0] fname_2 = arg[-1]
model_name=arg[1] df = pfam_parse(fname_1)
no_of_model=int(arg[2]) df.to_hdf("pfam_db_positive.h5", key="pfam_db_positive", mode="w")
nop=int(arg[3]) df = ""
st=int(arg[4]) gc.collect()
end=int(arg[5]) df = pfam_parse(fname_2)
name=arg[6] df.to_hdf("pfam_db_negative.h5", key="pfam_db_negative", mode="w")
pfam_fname=arg[7] print("Pfam Databases Written Successfully :)")
weight=arg[8]
if weight=="e":
w=[1]*no_of_model
else:
w=[]
for i in range(no_of_model):
w.append(float(arg[9+i]))
df=read_database(fname)
df=select_data_by_length(df,st,end)
n=3
if os.path.exists("prediction_data/data_"+fname):
with open("prediction_data/data_"+fname,"rb") as file: if __name__ == "__main__":
save_dict=pickle.load(file) main()
smg=save_dict["smg"]
sml=save_dict["sml"]
pfam_matrices=save_dict["pfam"]
indexes=save_dict["indexes"]
bls=save_dict["bls"]
blhs=save_dict["blhs"]
dis=save_dict["dis"]
dps=save_dict["dps"]
dphs=save_dict["dphs"]
else:
a,d,ld,ldg,cmap,cimap=read_genome_maps()#read the genome maps
print("Genome Maps Loaded.")
a_h=[df]
d_h=["prediction"]
lsy,gene_sequences=create_synteny_features(a_h,d_h,n,a,d,ld,ldg,cmap,cimap,name)
sml,smg,indexes=threadmaker(nop,df,lsy,gene_sequences,n,name)
a=""
d=""
ld=""
ldg=""
cmap=""
cimap=""
gene_sequences=""
gc.collect()
df_temp=df.loc[indexes]
pfam_db=pd.read_hdf(pfam_fname.split(".")[0]+"_pfam_db.h5")
with open(pfam_fname.split(".")[0]+"_pfam_map","rb") as file:
pfam_map=pickle.load(file)
pfam_matrices,indexes_pfam=create_pfam_matrix(df_temp,lsy,pfam_db,pfam_map)
pfam_db=""
gc.collect()
bls,blhs,dis,dps,dphs=create_tree_data("species_tree.tree",df_temp)
save_dict=dict(smg=smg,sml=sml,pfam=pfam_matrices,indexes=indexes,bls=bls,blhs=blhs,dis=dis,dps=dps,dphs=dphs)
if not os.path.isdir("prediction_data"):
os.mkdir("prediction_data")
with open("prediction_data/data_"+fname,"wb") as file:
pickle.dump(save_dict,file)
index_dict=create_map_list(indexes)
preds=get_prediction(smg,sml,pfam_matrices,indexes,bls,blhs,dis,dps,dphs,model_name,no_of_model,w)
write_preds(fname,model_name,name,preds,index_dict,df)
if __name__=="__main__":
main()

132
prepare_synteny_matrix.py
View file

@ -1,5 +1,4 @@
import numpy as np import numpy as np
import pandas as pd
import json import json
import os import os
import sys import sys
@ -7,95 +6,108 @@ import pickle
from select_data import read_db_homology from select_data import read_db_homology
from threads import Procerssrunner from threads import Procerssrunner
from read_get_gene_seq import read_gene_sequences from read_get_gene_seq import read_gene_sequences
from access_data_rest import update_rest,update_rest_protein from access_data_rest import update_rest, update_rest_protein
from Bio import SeqIO from Bio import SeqIO
from Bio.Seq import Seq from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC from Bio.Alphabet import IUPAC
def write_fasta(sequences,name):
with open(name+".fa","w") as file:
for seq in sequences:
if sequences[seq]=="":
continue
record=SeqRecord(Seq(sequences[seq],IUPAC.protein),id=seq)
SeqIO.write(record,file,"fasta")
def read_data_synteny(nop,name): def write_fasta(sequences, name):
smg=[] with open(name+".fa", "w") as file:
sml=[] for seq in sequences:
indexes=[] if sequences[seq] == "":
continue
record = SeqRecord(Seq(sequences[seq], IUPAC.protein), id=seq)
SeqIO.write(record, file, "fasta")
def read_data_synteny(nop, name):
smg = []
sml = []
indexes = []
for i in range(nop): for i in range(nop):
try: try:
with open("temp_"+name+"/thread_"+str(i+1)+"_smg.temp","rb") as file: with open("temp_"+name+"/thread_"+str(i+1) +
smg=smg+pickle.load(file) "_smg.temp", "rb") as file:
with open("temp_"+name+"/thread_"+str(i+1)+"_sml.temp","rb") as file: smg = smg+pickle.load(file)
sml=sml+pickle.load(file) with open("temp_"+name+"/thread_"+str(i+1) +
with open("temp_"+name+"/thread_"+str(i+1)+"_indexes.temp","rb") as file: "_sml.temp", "rb") as file:
indexes=indexes+pickle.load(file) sml = sml+pickle.load(file)
with open("temp_"+name+"/thread_"+str(i+1) +
"_indexes.temp", "rb") as file:
indexes = indexes+pickle.load(file)
except Exception as e: except Exception as e:
print("Problem with thread",i+1,"detected for",name,e) print("Problem with thread", i+1, "detected for", name, e)
continue continue
print(len(indexes)) print(len(indexes))
return smg,sml,indexes return smg, sml, indexes
def load_neighbor_genes(): def load_neighbor_genes():
with open("processed/neighbor_genes.json","r") as file: with open("processed/neighbor_genes.json", "r") as file:
lsy=dict(json.load(file)) lsy = dict(json.load(file))
print(len(lsy)) print(len(lsy))
print("Neighbor Genes Loaded") print("Neighbor Genes Loaded")
return lsy return lsy
def read_data_homology(dirname): def read_data_homology(dirname):
lf=os.listdir(dirname) lf = os.listdir(dirname)
if len(lf)==0: if len(lf) == 0:
print("No Files in the Directory!!!!!!!") print("No Files in the Directory!!!!!!!")
sys.exit(1) sys.exit(1)
a_h=[] a_h = []
d_h=[] d_h = []
for x in lf: for x in lf:
df,n=read_db_homology(dirname,x) df, n = read_db_homology(dirname, x)
n=n.split()[0] n = n.split()[0]
try: try:
indexes=np.load("processed/"+n+"_selected_indexes.npy") indexes = np.load("processed/"+n+"_selected_indexes.npy")
except: except Exception:
print("Incomplete data for:",n) print("Incomplete data for:", n)
df=df.loc[indexes] df = df.loc[indexes]
a_h.append(df) a_h.append(df)
d_h.append(n) d_h.append(n)
return a_h,d_h return a_h, d_h
def main(): def main():
arg=sys.argv arg = sys.argv
nop=int(arg[-1]) nop = int(arg[-1])
n=3 n = 3
a_h,d_h=read_data_homology("data_homology") a_h, d_h = read_data_homology("data_homology")
print("Data Read") print("Data Read")
lsy=load_neighbor_genes() lsy = load_neighbor_genes()
gene_sequences=read_gene_sequences(a_h,lsy,"geneseq","gene_seq_positive") gene_sequences = read_gene_sequences(
gene_sequences=update_rest(gene_sequences,"gene_seq_positive") a_h, lsy, "geneseq", "gene_seq_positive")
gene_sequences = update_rest(gene_sequences, "gene_seq_positive")
print("Gene Sequences Loaded.") print("Gene Sequences Loaded.")
if not os.path.isdir("processed/synteny_matrices"): if not os.path.isdir("processed/synteny_matrices"):
os.mkdir("processed/synteny_matrices") os.mkdir("processed/synteny_matrices")
ndir="processed/synteny_matrices/" ndir = "processed/synteny_matrices/"
nf1="synteny_matrices_global" nf1 = "synteny_matrices_global"
nf2="synteny_matrices_local" nf2 = "synteny_matrices_local"
nf3="indexes" nf3 = "indexes"
for i in range(len(a_h)): for i in range(len(a_h)):
df=a_h[i] df = a_h[i]
part=len(df)//nop part = len(df)//nop
pr=Procerssrunner() pr = Procerssrunner()
pr.start_processes(nop,df,gene_sequences,lsy,part,n,d_h[i]) pr.start_processes(nop, df, gene_sequences, lsy, part, n, d_h[i])
smg,sml,indexes=read_data_synteny(nop,d_h[i]) smg, sml, indexes = read_data_synteny(nop, d_h[i])
print(len(indexes)) print(len(indexes))
np.save(ndir+str(d_h[i])+"_"+nf1,smg) np.save(ndir+str(d_h[i])+"_"+nf1, smg)
np.save(ndir+str(d_h[i])+"_"+nf2,sml) np.save(ndir+str(d_h[i])+"_"+nf2, sml)
np.save(ndir+str(d_h[i])+"_"+nf3,indexes) np.save(ndir+str(d_h[i])+"_"+nf3, indexes)
a_h[i]=df.loc[indexes] a_h[i] = df.loc[indexes]
print("Synteny Matrices Created Successfully :)") print("Synteny Matrices Created Successfully :)")
protein_sequences=read_gene_sequences(a_h,lsy,"pro_seq","pro_seq_positive") protein_sequences = read_gene_sequences(
protein_sequences=update_rest_protein(protein_sequences,"pro_seq_positive") a_h, lsy, "pro_seq", "pro_seq_positive")
write_fasta(protein_sequences,"protein_seq_positive") protein_sequences = update_rest_protein(
protein_sequences, "pro_seq_positive")
write_fasta(protein_sequences, "protein_seq_positive")
print("Protein Sequences Loaded.") print("Protein Sequences Loaded.")
if __name__=="__main__":
if __name__ == "__main__":
main() main()

24
process_data.py
View file

@ -3,10 +3,10 @@ import progressbar
from save_data import write_dict_json from save_data import write_dict_json
def create_map_list(l): # this function maps the indexes to values def create_map_list(lister): # this function maps the indexes to values
t = {} t = {}
for i in range(len(l)): for i in range(len(lister)):
t[l[i]] = i t[lister[i]] = i
return t return t
@ -58,10 +58,15 @@ def get_nearest_neighbors(g, gs, n, a, d, ld, ldg, cmap, cimap):
try: try:
_ = ld[gi] # see if the corresponding gene map exists _ = ld[gi] # see if the corresponding gene map exists
except BaseException: except BaseException:
# print("Length of Dataframes:{} \t Length of Loaded Genes:{} \t Length of Loaded Genomes Dictionaries:{}".format(len(a),len(ld),len(ldg))) # print("Length of Dataframes:{}
# Length of Loaded Genes:{}
# Length of Loaded Genomes Dictionaries:
# {}".format(len(a),len(ld),len(ldg)))
return ne, nr return ne, nr
sldg = ldg[gi] # select the corresponding map sldg = ldg[gi] # select the corresponding map
if g not in sldg: # if the gene is not present in the dataframe return empty lists # if the gene is not present in the dataframe
# return empty lists
if g not in sldg:
# print(g,"\t",gs) # print(g,"\t",gs)
return ne, nr return ne, nr
i = sldg[g] # find the index of the gene i = sldg[g] # find the index of the gene
@ -80,9 +85,11 @@ def get_nearest_neighbors(g, gs, n, a, d, ld, ldg, cmap, cimap):
end = list(sldf.end) end = list(sldf.end)
end = np.array(end) end = np.array(end)
assert(len(end) == len(sldf)) assert(len(end) == len(sldf))
end = end - start # subtract start from it so as to get relative position # subtract start from it so as to get relative position
end = end - start
end_s = np.argsort(end) # sort them by the order of distance 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 # if all the genes end ahead of the one in consideration
if end[end_s[0]] >= 0:
flag = 1 # increment the pointer flag = 1 # increment the pointer
ne.append("NULL_GENE") # append the NULL_GENE value ne.append("NULL_GENE") # append the NULL_GENE value
continue continue
@ -144,7 +151,8 @@ def create_data_homology_ls(a_h, d_h, n, a, d, ld, ldg, cmap, cimap, to_write):
xl, xr = get_nearest_neighbors( xl, xr = get_nearest_neighbors(
x, xs, n, a, d, ld, ldg, cmap, cimap) x, xs, n, a, d, ld, ldg, cmap, cimap)
if len( if len(
xl) != 0: # check if neighboring genes were successfully found # check if neighboring genes were successfully found
xl) != 0:
lsy[x] = dict(b=xl, f=xr) lsy[x] = dict(b=xl, f=xr)
lsytemp[x] = dict(b=xl, f=xr) lsytemp[x] = dict(b=xl, f=xr)
except BaseException: except BaseException:

82
process_negative.py
View file

@ -1,63 +1,63 @@
import pandas as pd import pandas as pd
import numpy as np import numpy as np
import os
import sys
import progressbar
import json
import sys import sys
from neighbor_genes import read_genome_maps from neighbor_genes import read_genome_maps
from process_data import create_data_homology_ls from process_data import create_data_homology_ls
from threads import Procerssrunner from threads import Procerssrunner
from read_get_gene_seq import read_gene_sequences from read_get_gene_seq import read_gene_sequences
from access_data_rest import update_rest,update_rest_protein from access_data_rest import update_rest, update_rest_protein
from prepare_synteny_matrix import read_data_synteny,write_fasta from prepare_synteny_matrix import read_data_synteny, write_fasta
from save_data import write_dict_json from save_data import write_dict_json
from access_data_rest import update_rest_protein
def read_database_txt(filename): def read_database_txt(filename):
df=pd.read_csv(filename,sep="\t",header=None) df = pd.read_csv(filename, sep="\t", header=None)
df=df.drop(0,axis=1) df = df.drop(0, axis=1)
df.columns=["gene_stable_id","species","homology_gene_stable_id","homology_species","wga","goc","homology_type"] df.columns = ["gene_stable_id", "species", "homology_gene_stable_id",
"homology_species", "wga", "goc", "homology_type"]
return df return df
def main(): def main():
arg=sys.argv arg = sys.argv
a,d,ld,ldg,cmap,cimap=read_genome_maps() a, d, ld, ldg, cmap, cimap = read_genome_maps()
print("Genome Maps Loaded.") print("Genome Maps Loaded.")
df=read_database_txt(arg[-2]) df = read_database_txt(arg[-2])
nop=int(arg[-1]) nop = int(arg[-1])
print("Data Read.") print("Data Read.")
a_h=[] a_h = []
d_h=[] d_h = []
a_h.append(df) a_h.append(df)
d_h.append(arg[-2].split(".")[0]) d_h.append(arg[-2].split(".")[0])
n=3 n = 3
lsy=create_data_homology_ls(a_h,d_h,n,a,d,ld,ldg,cmap,cimap,0) lsy = create_data_homology_ls(a_h, d_h, n, a, d, ld, ldg, cmap, cimap, 0)
write_dict_json("neighbor_genes_negative","processed",lsy) write_dict_json("neighbor_genes_negative", "processed", lsy)
print("Neighbor Genes Found and Saved Successfully:)") print("Neighbor Genes Found and Saved Successfully:)")
gene_sequences=read_gene_sequences(a_h,lsy,"geneseq","gene_seq_negative") gene_sequences = read_gene_sequences(
gene_sequences=update_rest(gene_sequences,"gene_seq_negative") a_h, lsy, "geneseq", "gene_seq_negative")
ndir="processed/synteny_matrices/" gene_sequences = update_rest(gene_sequences, "gene_seq_negative")
nf1="synteny_matrices_global" ndir = "processed/synteny_matrices/"
nf2="synteny_matrices_local" nf1 = "synteny_matrices_global"
nf3="indexes" nf2 = "synteny_matrices_local"
nf3 = "indexes"
for i in range(len(a_h)): for i in range(len(a_h)):
df=a_h[i] df = a_h[i]
part=len(df)//nop part = len(df)//nop
pr=Procerssrunner() pr = Procerssrunner()
pr.start_processes(nop,df,gene_sequences,lsy,part,n,d_h[i]) pr.start_processes(nop, df, gene_sequences, lsy, part, n, d_h[i])
smg,sml,indexes=read_data_synteny(nop,d_h[i]) smg, sml, indexes = read_data_synteny(nop, d_h[i])
print(len(indexes)) print(len(indexes))
np.save(ndir+str(d_h[i])+"_"+nf1,smg) np.save(ndir+str(d_h[i])+"_"+nf1, smg)
np.save(ndir+str(d_h[i])+"_"+nf2,sml) np.save(ndir+str(d_h[i])+"_"+nf2, sml)
np.save(ndir+str(d_h[i])+"_"+nf3,indexes) np.save(ndir+str(d_h[i])+"_"+nf3, indexes)
a_h[i]=df.loc[indexes] a_h[i] = df.loc[indexes]
print("Synteny Matrices Created Successfully :)") print("Synteny Matrices Created Successfully :)")
protein_sequences=read_gene_sequences(a_h,lsy,"pro_seq","pro_seq_negative") protein_sequences = read_gene_sequences(
protein_sequences=update_rest_protein(protein_sequences,"pro_seq_negative") a_h, lsy, "pro_seq", "pro_seq_negative")
write_fasta(protein_sequences,"pro_seq_negative") protein_sequences = update_rest_protein(
protein_sequences, "pro_seq_negative")
write_fasta(protein_sequences, "pro_seq_negative")
if __name__=="__main__":
if __name__ == "__main__":
main() main()

3
read_data.py
View file

@ -63,7 +63,8 @@ def read_data_genome(dir_name, a, dict_ind_genome):
print(e) print(e)
continue continue
# print(data_gene[0:10]) # print(data_gene[0:10])
data_gene = data_gene[(data_gene['gene_biotype'] == 'protein_coding') | ( data_gene = data_gene[(
data_gene['gene_biotype'] == 'protein_coding') | (
data_gene['gene_source'] == 'protein_coding')] data_gene['gene_source'] == 'protein_coding')]
# print(data_gene[data_gene["gene_id"]=="ENSNGAG00000000407"]) # print(data_gene[data_gene["gene_id"]=="ENSNGAG00000000407"])
a.append(data_gene) a.append(data_gene)

20
read_get_gene_seq.py
View file

@ -27,7 +27,8 @@ def create_dict(keys, values, dictionary):
return dictionary return dictionary
# this function maps all the genes to their respective species. # 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) # (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):
@ -64,7 +65,9 @@ def read_gene_seq(dirname, s, genes_by_species):
ftr = [] ftr = []
for f in lof: for f in lof:
if f.split(".")[ 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 # check whether the species is present in the species to read list.
# Will skip those species which are not present in the dataframe
0] in s:
ftr.append(f) ftr.append(f)
data = {} data = {}
for f in progressbar.progressbar(ftr): for f in progressbar.progressbar(ftr):
@ -73,17 +76,18 @@ def read_gene_seq(dirname, s, genes_by_species):
record = SeqIO.parse(file, "fasta") record = SeqIO.parse(file, "fasta")
for r in record: for r in record:
gid, gbt = description_cleaner(r.description) gid, gbt = description_cleaner(r.description)
if str(gid) not in data and str( if str(gid) not in data and str(gid) in \
gid) in genes_by_species[species] and gbt == "protein_coding": genes_by_species[species] and gbt == "protein_coding":
data[gid] = str(r.seq) data[gid] = str(r.seq)
return data return data
def read_gene_sequences(hdf, 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. """The basic idea here is to create a list/dictionary of all the genes by
Once the mapping is done, all the respective fasta sequence files are read by Species their species. Once the mapping is done, all the respective fasta sequence
and the CDNA sequences for each gene in the species record are read and stored. files are read by Species and the CDNA sequences for each gene in the
Thus we don't have to read the same file multiple times.""" species record are read and stored.
Thus we don't have to read the same file multiple times."""
grouped_genes = {} grouped_genes = {}
gene_by_species_dict = {} gene_by_species_dict = {}
for df in progressbar.progressbar(hdf): for df in progressbar.progressbar(hdf):

9
save_data.py
View file

@ -15,9 +15,12 @@ def write_dict_json(name, dir, d):
def write_data_synteny(smg, sml, indexes, i, name): def write_data_synteny(smg, sml, indexes, i, name):
if not os.path.exists("temp_" + name): if not os.path.exists("temp_" + name):
os.mkdir("temp_" + name) os.mkdir("temp_" + name)
with open("temp_" + name + "/thread_" + str(i + 1) + "_smg.temp", "wb") as file: with open("temp_" + name + "/thread_" + str(i + 1) +
"_smg.temp", "wb") as file:
pickle.dump(smg, file) pickle.dump(smg, file)
with open("temp_" + name + "/thread_" + str(i + 1) + "_sml.temp", "wb") as file: with open("temp_" + name + "/thread_" + str(i + 1) +
"_sml.temp", "wb") as file:
pickle.dump(sml, file) pickle.dump(sml, file)
with open("temp_" + name + "/thread_" + str(i + 1) + "_indexes.temp", "wb") as file: with open("temp_" + name + "/thread_" + str(i + 1) +
"_indexes.temp", "wb") as file:
pickle.dump(indexes, file) pickle.dump(indexes, file)

9
threads.py
View file

@ -52,10 +52,12 @@ class Thread_objects():
norm_len = max(len(gene_seq[g1[i]]), len(gene_seq[g2[j]])) norm_len = max(len(gene_seq[g1[i]]), len(gene_seq[g2[j]]))
try: try:
result = ed.align( result = ed.align(
gene_seq[g1[i]], gene_seq[g2[j]], mode="NW", task="distance") gene_seq[g1[i]], gene_seq[g2[j]],
mode="NW", task="distance")
sm[i][j][0] = result["editDistance"] / (norm_len) sm[i][j][0] = result["editDistance"] / (norm_len)
result = ed.align( result = ed.align(
gene_seq[g1[i]], gene_seq[g2[j]][::-1], mode="NW", task="distance") gene_seq[g1[i]], gene_seq[g2[j]][::-1], mode="NW",
task="distance")
sm[i][j][1] = result["editDistance"] / (norm_len) sm[i][j][1] = result["editDistance"] / (norm_len)
_, result, _ = local_pairwise_align_ssw( _, result, _ = local_pairwise_align_ssw(
DNA(gene_seq[g1[i]]), DNA(gene_seq[g2[j]])) DNA(gene_seq[g1[i]]), DNA(gene_seq[g2[j]]))
@ -120,7 +122,8 @@ class Procerssrunner():
def start_thread(self, obj, i, thread_alive, n, name): def start_thread(self, obj, i, thread_alive, n, name):
t = Process(target=obj.synteny_matrix, args=( t = Process(target=obj.synteny_matrix, args=(
obj.gene_sequences, obj.df, obj.lsy, n), name="Thread_" + str(i + 1)) obj.gene_sequences, obj.df, obj.lsy, n),
name="Thread_" + str(i + 1))
print("Thread ", (i + 1), " started for ", name, ".") print("Thread ", (i + 1), " started for ", name, ".")
thread_alive.append(t) thread_alive.append(t)

293
train.py
View file

@ -4,87 +4,119 @@ import sys
import tensorflow as tf import tensorflow as tf
from model import create_model from model import create_model
def create_branch_length_padding(bl): def create_branch_length_padding(bl):
maxlen=0 maxlen = 0
for x in bl: for x in bl:
if len(x)>maxlen: if len(x) > maxlen:
maxlen=len(x) maxlen = len(x)
for x in range(len(bl)): for x in range(len(bl)):
temp=bl[x] temp = bl[x]
for i in range(len(temp),maxlen): for i in range(len(temp), maxlen):
temp=np.append(temp,[0]) temp = np.append(temp, [0])
bl[x]=temp bl[x] = temp
return bl return bl
def train(train_synteny_matrices_global,train_synteny_matrices_local,train_pfam_matrices,train_branch_length_species,train_branch_length_homology_species,train_dist_p_s,train_dist_p_hs,train_distance,train_labels,v,num_epochs,learning_rate,decay,size_train,batch_size,model_name):
print("Going to train model {} for:\n Batch Size:{} \n Learning Rate:{} \n Decay:{}\n On {} Samples".format(v,batch_size,learning_rate,decay,len(train_synteny_matrices_global))) def train(train_synteny_matrices_global, train_synteny_matrices_local,
graph,saver=create_model() train_pfam_matrices, train_branch_length_species,
synmg,synml,pfam,bls,blhs,dps,dphs,dis,lr,y=graph.get_collection("input_nodes") train_branch_length_homology_species, train_dist_p_s,
loss,t_op,accuracy,init,summary=graph.get_collection("output_nodes") train_dist_p_hs, train_distance, train_labels, v,
num_epochs, learning_rate, decay, size_train,
batch_size, model_name):
print("Going to train model {} for:\n Batch Size:{} \
\n Learning Rate:{} \n Decay:{}\n On {} Samples".format(
v, batch_size, learning_rate, decay,
len(train_synteny_matrices_global)
))
graph, saver = create_model()
synmg, synml, pfam, bls, blhs, dps, \
dphs, dis, lr, y = graph.get_collection("input_nodes")
loss, t_op, accuracy, init, summary = graph.get_collection("output_nodes")
with tf.Session(graph=graph) as sess: with tf.Session(graph=graph) as sess:
writer = tf.summary.FileWriter('./'+model_name+'_v'+str(v), sess.graph) writer = tf.summary.FileWriter('./'+model_name+'_v'+str(v), sess.graph)
sess.run(init) sess.run(init)
learn=learning_rate learn = learning_rate
for j in range(num_epochs): for j in range(num_epochs):
for i in range(size_train//batch_size): for i in range(size_train//batch_size):
feed_dict={ feed_dict = {
synmg:train_synteny_matrices_global[i*batch_size:(i+1)*batch_size], synmg: train_synteny_matrices_global
synml:train_synteny_matrices_local[i*batch_size:(i+1)*batch_size], [i*batch_size:(i+1)*batch_size],
pfam:train_pfam_matrices[i*batch_size:(i+1)*batch_size].reshape((batch_size,7,7,1)), synml: train_synteny_matrices_local
bls:train_branch_length_species[i*batch_size:(i+1)*batch_size], [i*batch_size:(i+1)*batch_size],
blhs:train_branch_length_homology_species[i*batch_size:(i+1)*batch_size], pfam: train_pfam_matrices
dps:train_dist_p_s[i*batch_size:(i+1)*batch_size].reshape((batch_size,1)), [i*batch_size:(i+1)*batch_size]
dphs:train_dist_p_hs[i*batch_size:(i+1)*batch_size].reshape((batch_size,1)), .reshape((batch_size, 7, 7, 1)),
dis:train_distance[i*batch_size:(i+1)*batch_size].reshape((batch_size,1)), bls: train_branch_length_species
lr:learn, [i*batch_size:(i+1)*batch_size],
y:train_labels[i*batch_size:(i+1)*batch_size] blhs: train_branch_length_homology_species
[i*batch_size:(i+1)*batch_size],
dps: train_dist_p_s
[i*batch_size:(i+1)*batch_size]
.reshape((batch_size, 1)),
dphs: train_dist_p_hs
[i*batch_size:(i+1)*batch_size]
.reshape((batch_size, 1)),
dis: train_distance
[i*batch_size:(i+1)*batch_size]
.reshape((batch_size, 1)),
lr: learn,
y: train_labels[i*batch_size:(i+1)*batch_size]
} }
sess.run(t_op,feed_dict=feed_dict) sess.run(t_op, feed_dict=feed_dict)
test_dict={ test_dict = {
synmg:train_synteny_matrices_global[size_train:], synmg: train_synteny_matrices_global[size_train:],
synml:train_synteny_matrices_local[size_train:], synml: train_synteny_matrices_local[size_train:],
pfam:train_pfam_matrices[size_train:].reshape((len(train_synteny_matrices_global)-size_train,7,7,1)), pfam: train_pfam_matrices[size_train:]
bls:train_branch_length_species[size_train:], .reshape((
blhs:train_branch_length_homology_species[size_train:], len(train_synteny_matrices_global)-size_train, 7, 7, 1)),
dps:train_dist_p_s[size_train:].reshape((len(train_synteny_matrices_global)-size_train,1)), bls: train_branch_length_species[size_train:],
dphs:train_dist_p_hs[size_train:].reshape((len(train_synteny_matrices_global)-size_train,1)), blhs: train_branch_length_homology_species[size_train:],
dis:train_distance[size_train:].reshape((len(train_synteny_matrices_global)-size_train,1)), dps: train_dist_p_s[size_train:]
y:train_labels[size_train:], .reshape((len(train_synteny_matrices_global)-size_train, 1)),
lr:learn dphs: train_dist_p_hs[size_train:]
} .reshape((len(train_synteny_matrices_global)-size_train, 1)),
accuracy_test,loss_test,summary_write=sess.run([accuracy,loss,summary],feed_dict=test_dict) dis: train_distance[size_train:]
writer.add_summary(summary_write,i+1) .reshape((len(train_synteny_matrices_global)-size_train, 1)),
print("Epoch:{} Test Accuracy:{} Test Loss:{}".format(j+1,accuracy_test*100,loss_test)) y: train_labels[size_train:],
learn*=decay lr: learn
saver.save(sess,model_name+"_v"+str(v)+"/model.ckpt") }
accuracy_test, loss_test, summary_write = sess.run(
[accuracy, loss, summary], feed_dict=test_dict)
writer.add_summary(summary_write, i+1)
print("Epoch:{} Test Accuracy:{} Test Loss:{}".format(
j+1, accuracy_test*100, loss_test))
learn *= decay
saver.save(sess, model_name+"_v"+str(v)+"/model.ckpt")
writer.close() writer.close()
def read_positive(len_p,bls,blhs,dis,dps,dphs,sml,smg,pfam,label):
rowsh=[] def read_positive(len_p, bls, blhs, dis, dps, dphs, sml, smg, pfam, label):
with open("dataset","rb") as file: rowsh = []
rowsh=pickle.load(file) with open("dataset", "rb") as file:
shi=np.random.permutation(len(rowsh)) rowsh = pickle.load(file)
rows_shuffled=[] shi = np.random.permutation(len(rowsh))
rows_shuffled = []
for i in range(len(shi)): for i in range(len(shi)):
rows_shuffled.append(rowsh[shi[i]]) rows_shuffled.append(rowsh[shi[i]])
rowsh=rows_shuffled rowsh = rows_shuffled
spco={} spco = {}
spcp={} spcp = {}
for row in rowsh: for row in rowsh:
if row["species"] not in spco: if row["species"] not in spco:
spco[row["species"]]=0 spco[row["species"]] = 0
spcp[row["species"]]=0 spcp[row["species"]] = 0
maxcount_o=int((len_p)*0.3/14) maxcount_o = int((len_p)*0.3/14)
maxcount_p=int((len_p)*0.7/14) maxcount_p = int((len_p)*0.7/14)
for row in rowsh: for row in rowsh:
if row["label"]==2: if row["label"] == 2:
continue continue
if row["label"]==1 and spco[row["species"]]>maxcount_o: if row["label"] == 1 and spco[row["species"]] > maxcount_o:
continue continue
if row["label"]==0 and spcp[row["species"]]>maxcount_p: if row["label"] == 0 and spcp[row["species"]] > maxcount_p:
continue continue
bls.append(np.array(row["bls"])) bls.append(np.array(row["bls"]))
blhs.append(np.array(row["blhs"])) blhs.append(np.array(row["blhs"]))
@ -95,23 +127,23 @@ def read_positive(len_p,bls,blhs,dis,dps,dphs,sml,smg,pfam,label):
smg.append(row["global_alignment_matrix"]) smg.append(row["global_alignment_matrix"])
pfam.append(row["pfam_matrix"]) pfam.append(row["pfam_matrix"])
label.append(row["label"]) label.append(row["label"])
if row["label"]==1: if row["label"] == 1:
spco[row["species"]]+=1 spco[row["species"]] += 1
if row["label"]==0: if row["label"] == 0:
spcp[row["species"]]+=1 spcp[row["species"]] += 1
def read_negative(len_n,bls,blhs,dis,dps,dphs,sml,smg,pfam,label): def read_negative(len_n, bls, blhs, dis, dps, dphs, sml, smg, pfam, label):
rows=[] rows = []
with open("dataset","rb") as file: with open("dataset", "rb") as file:
rows=pickle.load(file) rows = pickle.load(file)
rows=[row for row in rows if row["label"]==2] rows = [row for row in rows if row["label"] == 2]
shi=np.random.permutation(len(rows)) shi = np.random.permutation(len(rows))
rows_shuffled=[] rows_shuffled = []
for i in range(len(shi)): for i in range(len(shi)):
rows_shuffled.append(rows[shi[i]]) rows_shuffled.append(rows[shi[i]])
rows=rows_shuffled rows = rows_shuffled
rows=rows[:len_n] rows = rows[:len_n]
for row in rows: for row in rows:
bls.append(np.array(row["bls"])) bls.append(np.array(row["bls"]))
blhs.append(np.array(row["blhs"])) blhs.append(np.array(row["blhs"]))
@ -123,70 +155,73 @@ def read_negative(len_n,bls,blhs,dis,dps,dphs,sml,smg,pfam,label):
label.append(row["label"]) label.append(row["label"])
pfam.append(row["pfam_matrix"]) pfam.append(row["pfam_matrix"])
def train_models(model_name,start,end,num_epochs,learn_rate,decay,size_train,batch_size):
k=1 def train_models(model_name, start, end, num_epochs, learn_rate,
for i in range(start//10,end//10+1): decay, size_train, batch_size):
bls=[] k = 1
blhs=[] for i in range(start//10, end//10+1):
dis=[] bls = []
dps=[] blhs = []
dphs=[] dis = []
sml=[] dps = []
smg=[] dphs = []
pfam=[] sml = []
label=[] smg = []
portion=float(i/10) pfam = []
len_n=int(size_train*portion) label = []
len_p=int(size_train*(1-portion)) portion = float(i/10)
read_positive(len_p,bls,blhs,dis,dps,dphs,sml,smg,pfam,label) len_n = int(size_train*portion)
read_negative(len_n,bls,blhs,dis,dps,dphs,sml,smg,pfam,label) len_p = int(size_train*(1-portion))
bls=create_branch_length_padding(bls) read_positive(len_p, bls, blhs, dis, dps, dphs, sml, smg, pfam, label)
blhs=create_branch_length_padding(blhs) read_negative(len_n, bls, blhs, dis, dps, dphs, sml, smg, pfam, label)
bls=np.array(bls) bls = create_branch_length_padding(bls)
blhs = create_branch_length_padding(blhs)
bls = np.array(bls)
print(bls.shape) print(bls.shape)
blhs=np.array(blhs) blhs = np.array(blhs)
print(blhs.shape) print(blhs.shape)
dis=np.array(dis) dis = np.array(dis)
print(dis.shape) print(dis.shape)
dps=np.array(dps) dps = np.array(dps)
print(dps.shape) print(dps.shape)
dphs=np.array(dphs) dphs = np.array(dphs)
print(dphs.shape) print(dphs.shape)
sml=np.array(sml) sml = np.array(sml)
print(sml.shape) print(sml.shape)
smg=np.array(smg) smg = np.array(smg)
print(smg.shape) print(smg.shape)
pfam=np.array(pfam) pfam = np.array(pfam)
print(pfam.shape) print(pfam.shape)
label=np.array(label) label = np.array(label)
print(label.shape) print(label.shape)
shi=np.random.permutation(len(label)) shi = np.random.permutation(len(label))
labels=label[shi] labels = label[shi]
bls=bls[shi] bls = bls[shi]
blhs=blhs[shi] blhs = blhs[shi]
dis=dis[shi] dis = dis[shi]
dps=dps[shi] dps = dps[shi]
dphs=dphs[shi] dphs = dphs[shi]
sml=sml[shi] sml = sml[shi]
smg=smg[shi] smg = smg[shi]
pfam=pfam[shi] pfam = pfam[shi]
train(smg,sml,pfam,bls,blhs,dps,dphs,dis,labels,k,num_epochs,learn_rate,decay,int(0.9*size_train),batch_size,model_name) train(smg, sml, pfam, bls, blhs, dps, dphs, dis, labels, k, num_epochs,
k+=1 learn_rate, decay, int(0.9*size_train), batch_size, model_name)
k += 1
def main(): def main():
arg=sys.argv arg = sys.argv
model_name=arg[-8] model_name = arg[-8]
start_p=int(arg[-7]) start_p = int(arg[-7])
end_p=int(arg[-6]) end_p = int(arg[-6])
num_epochs=int(arg[-5]) num_epochs = int(arg[-5])
learn_rate=float(arg[-4]) learn_rate = float(arg[-4])
decay=float(arg[-3]) decay = float(arg[-3])
size_train=float(arg[-2]) size_train = float(arg[-2])
batch_size=int(arg[-1]) batch_size = int(arg[-1])
train_models(model_name,start_p,end_p,num_epochs,learn_rate,decay,size_train,batch_size) train_models(model_name, start_p, end_p, num_epochs,
learn_rate, decay, size_train, batch_size)
if __name__=="__main__":
if __name__ == "__main__":
main() main()