Goal¶
- benchmark the k-mer bias model
Tasks¶
- [x] get the training data
- get the default 6-mer model
- [x] make an abstract bias prediction class
- [x] evaluate the default model
- [ ] write a comprehensive evaluation function
- [ ] train your own bias-prediction 6-mer model
- [ ] evaluate
- [ ] do they match? If now. Why?
- [ ] Evaluate your own model
- [ ] Imrove your own model
get the training data¶
In [1]:
import matplotlib.pyplot as plt
import basepair
from basepair.cli.schemas import DataSpec, TaskSpec
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()
In [2]:
mkdir -p {ddir}/raw/dnase/6mer
In [460]:
!wget https://raw.githubusercontent.com/jpiper/pyDNase/master/pyDNase/data/IMR90_6mer.pickle -O {ddir}/raw/dnase/6mer/IMR90_6mer.pickle
In [3]:
from basepair.utils import read_pkl
In [4]:
# Load the kmer file
kmer_model = read_pkl(f"{ddir}/raw/dnase/6mer/IMR90_6mer.pickle")
# dict -> pd.DataFrame
df = pd.DataFrame([{"motif": k, **v} for k,v in kmer_model.items()])
df['forward'] = df['forward'].astype(float)
df['reverse'] = df['reverse'].astype(float)
In [5]:
kmer_model['AAAAAA']
Out[5]:
In [6]:
len(df) # same as 4**6
Out[6]:
In [7]:
df.head()
Out[7]:
In [4]:
from basepair.utils import reverse_complement
In [9]:
df['rc_motif'] = df.motif.map(reverse_complement)
In [10]:
df.head()
Out[10]:
In [12]:
df.forward.plot.hist(bins=50);
plt.xlabel("Forward");
In [13]:
df.plot.scatter('forward', 'reverse', alpha=0.3);
Why are there forward and reverse strands?
In [14]:
dfm = pd.merge(df[['motif', 'forward']], df[['rc_motif', 'reverse']], left_on='motif', right_on='rc_motif')
In [15]:
dfm.plot.scatter('forward', 'reverse', alpha=0.3);
In [16]:
np.allclose(dfm.forward, dfm.reverse)
Out[16]:
Both forward and reverse scores are symmetrically stored
make an abstract bias prediction class¶
In [143]:
from basepair.exp.dnase.models import KmerBiasModel
In [87]:
kmb.counts.get("asd", 0)
Out[87]:
In [88]:
kmb = KmerBiasModel()
In [89]:
len("ACGTAGTCGAG")
Out[89]:
In [90]:
kmb.fit_on_batch(encodeDNA(["ACGTAGTCGAG"]*10), np.ones((10,11)))
In [93]:
kmb.normalize(False)
In [94]:
kmb.d
Out[94]:
In [91]:
kmb.normalize()
In [92]:
kmb.d
Out[92]:
In [57]:
{k:v for k,v in kmb.counts.items()}
Out[57]:
In [55]:
sum(kmb.counts.values())
Out[55]:
In [51]:
d = {k: float(v['forward']) for k,v in kmer_model.items()}
In [21]:
km = KmerBiasModel(d)
In [22]:
km.predict_seq("ACGTAGTCGAG", pad=True)
Out[22]:
In [23]:
km.predict_seq("ACGTAGTCGAG", pad=False)
Out[23]:
In [24]:
km.predict(encodeDNA(["ACGTAGTCGAG"])[0], pad=False)
Out[24]:
In [25]:
km.predict_on_batch(encodeDNA(["ACGTAGTCGAG"]), pad=False)
Out[25]:
evaluate the default model¶
- total count prediction in a window
- single base prediction
In [98]:
dpath = Path("/mnt/lab_data/users/avsec/dnase_bias")
SEQ_WIDTH = 1000
ds = DataSpec(
task_specs={"DNase": TaskSpec(task="DNase",
pos_counts=dpath / "IMR90_Naked_DNase.pos.bw",
neg_counts=dpath / "IMR90_Naked_DNase.neg.bw",
peaks=dpath / "genome.stride200.w200.nonblacklist.bed")},
fasta_file=dpath / "GRCh38.genome.fa",
)
In [99]:
valid_chr=['chr2', 'chr3', 'chr4']
test_chr=['chr1', 'chr8', 'chr9']
In [101]:
from basepair.datasets import StrandedProfile
class AppendCounts:
def transform(self, x):
counts = {k.replace("profile/", "counts/"): np.log(1 + x[k].sum(0))
for k in x}
return {**x, **counts}
In [102]:
SEQ_WIDTH=1000
In [113]:
dl = StrandedProfile(ds,
excl_chromosomes=valid_chr+test_chr,
peak_width=SEQ_WIDTH, shuffle=True,
target_transformer=AppendCounts())
In [32]:
dl_valid = StrandedProfile(ds,
incl_chromosomes=valid_chr,
peak_width=SEQ_WIDTH,
shuffle=True,
target_transformer=AppendCounts())
In [103]:
dl_test = StrandedProfile(ds,
incl_chromosomes=test_chr,
peak_width=SEQ_WIDTH,
shuffle=True,
target_transformer=AppendCounts())
In [34]:
from kipoi.data_utils import numpy_collate_concat
In [36]:
def touch_file(file, verbose=True):
if verbose:
add = "v"
else:
add = ""
!vmtouch -{add}tf {file}
def touch_all_files(ds, verbose=True):
touch_file(ds.fasta_file, verbose)
for ts in ds.task_specs.values():
touch_file(ts.pos_counts, verbose)
touch_file(ts.neg_counts, verbose)
#touch_file(ts.peaks)
In [37]:
touch_all_files(ds)
In [3]:
from tqdm import tqdm
In [110]:
from kipoi.writers import HDF5BatchWriter
In [111]:
!mkdir -p {ddir}/processed/dnase/6mer
In [42]:
#!rm {ddir}/processed/dnase/6mer/IMR90_6mer.preds.h5
In [107]:
dl_test[0]['targets']['profile/DNase'].shape
Out[107]:
In [124]:
%load_ext line_profiler
In [126]:
from numba import jit
In [133]:
new_kmer_model = KmerBiasModel(k=6)
In [136]:
%lprun -f new_kmer_model.fit_on_batch new_kmer_model.fit_on_batch(batch['inputs'], batch['targets']['profile/DNase'].sum(axis=-1))
Write a comprehensive evaluation function¶
- make a correlation plot of base-pair level counts
- make a correlation plot of 1kb level counts
- show some example profiles
In [33]:
from kipoi.readers import HDF5Reader
In [31]:
ls {ddir}/processed/dnase/6mer/
In [52]:
r = HDF5Reader(f"{ddir}/processed/dnase/6mer/IMR90_6mer.preds.h5")
In [53]:
r = HDF5Reader(f"{ddir}/processed/dnase/6mer/retrained.IMR90_6mer.w200.normalized.preds.h5")
r.open()
In [54]:
r.ls()
Out[54]:
In [70]:
preds = r.f['/preds'][:]
targets = r.f['/targets/profile/DNase'][:].sum(axis=-1)
y_true = targets[:,3:(1000-2)]
In [57]:
#y_true = targets[:,3:(200-1)]
In [72]:
assert preds.shape == y_true.shape
In [73]:
y_true.shape
Out[73]:
In [68]:
preds.shape
Out[68]:
In [74]:
y_true_l = np.ravel(y_true)
In [75]:
y_pred_l = np.ravel(preds)
In [76]:
y_true_l[:5]
Out[76]:
In [61]:
y_true_l.max()
Out[61]:
In [62]:
import random
In [63]:
idx = np.unique(random.sample(range(0, len(y_pred_l)), 1000000))
In [47]:
plt.scatter(y_true_l[idx], y_pred_l[idx], alpha=0.1)
plt.xlabel("true values")
plt.ylabel("predictions")
Out[47]:
Aggregated plot¶
In [77]:
from basepair.plots import regression_eval
In [78]:
yt_counts = np.log(1+y_true.sum(axis=-1))
yp_counts = np.log(1+preds.sum(axis=-1))
In [79]:
regression_eval(yt_counts, yp_counts, alpha=0.1)
In [50]:
regression_eval(yt_counts, yp_counts, alpha=0.1)
In [48]:
regression_eval(yt_counts, yp_counts, alpha=0.1)
In [ ]:
plt.scatter(yt_counts, yp_counts, alpha=0.1)
plt.xlabel("true values")
plt.ylabel("predictions")
In [34]:
from scipy.stats import pearsonr, spearmanr
In [35]:
yt = y_true_l
yp =y_pred_l
In [ ]:
rp = pearsonr(yt, yp)[0]
rs = spearmanr(yt, yp)[0]
In [ ]:
print(f"pearsonr: {rp}")
print(f"spearmanr: {rs}")
train your own bias-prediction 6-mer model¶
- done in the script
evaluate¶
do they match? If now. Why?¶
They don't match for some reason. Check again how they compute the vector.
In [2]:
from basepair.utils import read_pkl
In [3]:
km_new = read_pkl(f"{ddir}/processed/dnase/6mer/new_model.train.pkl")
In [4]:
ls {ddir}/processed/dnase/6mer/
In [20]:
km_new_pos3 = read_pkl(f"{ddir}/processed/dnase/6mer/new_model.train-pos3.w200.normalized.pkl")
In [21]:
from basepair.exp.dnase.models import KmerBiasModel
In [22]:
kmer_model = read_pkl(f"{ddir}/raw/dnase/6mer/IMR90_6mer.pickle")
km_old = KmerBiasModel({k: float(v['forward']) for k,v in kmer_model.items()})
In [23]:
df_new = pd.DataFrame([{"kmer": k, "value": v} for k,v in km_new.d.items()])
In [24]:
df_new_pos3 = pd.DataFrame([{"kmer": k, "value": v} for k,v in km_new_pos3.d.items()])
In [25]:
df_old = pd.DataFrame([{"kmer": k, "value": v} for k,v in km_old.d.items()])
In [26]:
dfeval = pd.merge(df_new, df_old, on='kmer', suffixes=('_new', '_old'))
In [27]:
dfeval2 = pd.merge(df_new_pos3, df_old, on='kmer', suffixes=('_new', '_old'))
In [28]:
plt.scatter(dfeval.value_new*100, dfeval.value_old, alpha=0.1)
#plt.xlim([0, 0.3])
#plt.ylim([0, 0.3])
plt.xlabel("value_new")
plt.ylabel("value_old")
Out[28]:
In [29]:
plt.scatter(np.log(1+dfeval2.value_new*100), np.log(1+dfeval.value_old), alpha=0.1)
plt.xlabel("value_new")
plt.ylabel("value_old")
Out[29]:
In [30]:
plt.scatter(dfeval2.value_new*100, dfeval.value_old, alpha=0.1)
#plt.xlim([0, 0.3])
#plt.ylim([0, 0.3])
plt.xlabel("value_new")
plt.ylabel("value_old")
Out[30]:
In [76]:
plt.scatter(dfeval.value_new*100, dfeval.kmer.map({k: float(v['forward']) for k,v in kmer_model.items()}), alpha=0.1)
#plt.xlim([0, 0.3])
#plt.ylim([0, 0.3])
plt.xlabel("value_new")
plt.ylabel("value_old")
Out[76]:
In [75]:
plt.scatter(dfeval.value_new*100, dfeval.kmer.map({k: float(v['reverse']) for k,v in kmer_model.items()}), alpha=0.1)
#plt.xlim([0, 0.3])
#plt.ylim([0, 0.3])
plt.xlabel("value_new")
plt.ylabel("value_old")
Out[75]:
In [69]:
dfeval.value_new.min()
Out[69]: