In [2]:
import matplotlib.pyplot as plt
import basepair
from basepair.plots import regression_eval
from basepair.losses import MultichannelMultinomialNLL
from basepair.exp.dnase.models import KmerBiasModel
from kipoi.readers import HDF5Reader
from basepair.utils import read_pkl
from keras.models import load_model
from basepair.cli.schemas import DataSpec, TaskSpec
from basepair.datasets import StrandedProfile
from basepair.datasets import chip_exo_nexus
from pathlib import Path
import pandas as pd
import numpy as np
from basepair.config import create_tf_session, get_data_dir
ddir = get_data_dir()
SEQ_WIDTH = 1000
In [3]:
ls {ddir}/processed/dnase/exp/models/single-task/
In [4]:
ls {ddir}/processed/dnase/exp/models/single-task-background-corrected/models
In [5]:
ls {ddir}/processed/dnase/exp/models/single-task-background-corrected/
In [6]:
!cat {ddir}/processed/dnase/exp/models/single-task-background-corrected/tryout.model.txt
In [7]:
mdir_root = Path(f"{ddir}//processed/dnase/exp/models/single-task-background-corrected/models")
mdirs = list(mdir_root.glob("*=*"))
mdirs
Out[7]:
In [8]:
# performance comparison
def get_performance(mdir):
try:
dfh = pd.read_csv(f"{mdir}/history.csv")
return dict(dfh.iloc[dfh.val_loss.idxmin()])
except:
return {}
df_exp = pd.DataFrame([{**get_performance(md), "name":md.name}
for md in mdirs])
df_exp.sort_values("val_counts/DNase_loss")
Out[8]:
In [9]:
!cat {str(mdir_root / "tconv_kernel_size=51/hparams.yaml")} | head
In [10]:
mdir_root
Out[10]:
In [11]:
!cat {str(mdir_root / "tconv_kernel_size=51/hparams.yaml")} | head
In [12]:
from keras.models import load_model
In [13]:
create_tf_session(0)
Out[13]:
In [14]:
mdir = mdir_root / "tconv_kernel_size=51"
model = load_model(mdir / "model.h5",
custom_objects={"mc_multinomial_nll": MultichannelMultinomialNLL(1),
"MultichannelMultinomialNLL": MultichannelMultinomialNLL(1)})
In [194]:
mdir = mdir_root / "../../single-task/models/tconv_kernel_size=51"
model = load_model(mdir / "model.h5",
custom_objects={"mc_multinomial_nll": MultichannelMultinomialNLL(1),
"MultichannelMultinomialNLL": MultichannelMultinomialNLL(1)})
In [28]:
from basepair.cli.modisco import compute_importance_scores
In [35]:
compute_importance_scores(mdir, "/tmp/imp.h5")
In [36]:
from kipoi.readers import HDF5Reader
r = HDF5Reader("/tmp/imp.h5")
r.open()
In [47]:
d = HDF5Reader.load("/tmp/imp.h5")
In [37]:
r.ls()
Out[37]:
In [42]:
%time d = r.load_all()
In [46]:
r.close()
In [42]:
from basepair.config import test_chr
from basepair.preproc import AppendCounts
ds = DataSpec.load(mdir / "dataspec.yaml")
ds.task_specs['DNase'].ignore_strand = True # fix
dl_test = StrandedProfile(ds,
incl_chromosomes=test_chr,
peak_width=SEQ_WIDTH,
shuffle=False,
target_transformer=AppendCounts())
In [43]:
test = dl_test.load_all(batch_size=256, num_workers=10)
In [44]:
from basepair.BPNet import BiasModel
In [45]:
ds.task_specs['DNase'].bias_model = "/users/avsec/workspace/basepair/data/processed/dnase/exp/models/background/models/w=1000/model.h5"
In [46]:
x,y_true = BiasModel(ds).predict((test["inputs"], test["targets"]), batch_size=512)
In [47]:
y_true.keys()
Out[47]:
In [48]:
from basepair.plots import regression_eval
from basepair.cli.evaluate import eval_profile
from basepair.math import softmax
In [92]:
y_pred = model.predict(x, verbose=1)
In [77]:
# Not bias-corrected
regression_eval(y_true['counts/DNase'].mean(-1),
y_pred[ds.task2idx("DNase", 'counts')].mean(-1),
alpha=0.1)
In [53]:
# bias-corrected
regression_eval(y_true['counts/DNase'].mean(-1),
y_pred[ds.task2idx("DNase", 'counts')].mean(-1),
alpha=0.1)
In [78]:
task = "DNase"
In [79]:
yp = softmax(y_pred[ds.task2idx(task, "profile")])
yt = y_true["profile/" + task]
In [80]:
df = eval_profile(yt, yp)
In [81]:
# not bias-corrected
df
Out[81]:
In [52]:
# bias-corrected
df
Out[52]:
In [101]:
from basepair.plots import plot_profiles
from basepair import samplers
In [56]:
tasks = list(ds.task_specs)
In [95]:
idx_list = samplers.top_sum_count(y_true['profile/DNase'], 10)
In [93]:
idx_list = samplers.top_max_count(y_true['profile/DNase'], 10)
In [97]:
idx_list = samplers.random(y_true['profile/DNase'], 10)
In [98]:
plot_profiles(y_true, y_pred, tasks, idx_list, preproc=None, figsize=(20, 2), tic_freq=50)
In [102]:
# get high bias sites
In [117]:
high_bias_idx = list(pd.Series(x['bias/counts/DNase'][:,0]).sort_values(ascending=False)[:10].index)
In [118]:
high_bias_idx
Out[118]:
In [ ]:
# TODO - plot the importance scores
In [119]:
plot_profiles(y_true, y_pred, tasks, high_bias_idx, preproc=None, figsize=(20, 2), tic_freq=50)
In [120]:
from basepair.BPNet import BPNetPredictor
In [195]:
bp = BPNetPredictor(model, ds.fasta_file, tasks, bias_model=BiasModel(ds))
In [122]:
from kipoi.data_utils import get_dataset_item
In [149]:
bp.model.inputs[0]
Out[149]:
In [130]:
assert bp.model.inputs[0].name[:3] == "seq"
In [136]:
bp.model.input_names
Out[136]:
In [132]:
bp.model.inputs
Out[132]:
In [150]:
bp.input_grad(get_dataset_item(x, high_bias_idx), strand='pos').shape
Out[150]:
In [1]:
from kipoi.writers import HDF5BatchWriter
In [165]:
from pybedtools import BedTool
In [176]:
bt = list(BedTool.from_dataframe(dl_test.dfm.iloc[high_bias_idx]))
In [159]:
list((1,2,3))
Out[159]:
In [163]:
list(({"seq": 1}, )[1:])
Out[163]:
In [196]:
# non bias-corrected
bp.plot_predict_grad(bt)
In [193]:
# bias-corrected
bp.plot_predict_grad(bt)
In [156]:
(1,) + (1,)[1:]
Out[156]:
In [146]:
bt.head()
In [ ]:
bp.plot_predict_grad()