In [1]:
%load_ext autoreload
%autoreload 2
import pickle
from collections import OrderedDict
import numpy as np
import pandas as pd
import os
import subprocess
import math
import copy
import pickle
import h5py
import glob
from matlas.genome_data import *
from matlas.pwms import shorten_tfname, get_motifDB_id_maps
from matlas.matches import fig2vdom, vdom_pssm
#from basepair.modisco.results import ModiscoResult
from matlas.modisco_results import ModiscoResult
project_root = "/mnt/lab_data/kundaje/msharmin/mouse_hem/with_tfd"
oak_root = "/oak/stanford/groups/akundaje/msharmin/mouse_hem/from_labcluster/with_tfd"
srv_root = "/srv/scratch/msharmin/mouse_hem/with_tfd"
data_name = "full_mouse50"
modisco_root = "{0}/{1}/Naive_modisco2019_msharmin".format(oak_root, data_name)
modisco_root = "{0}/{1}/Naive_modisco2019".format(srv_root, data_name)
homer_root = "{0}/{1}/Naive_scans".format(srv_root, data_name)
model_root = "{0}/{1}/fineFactorized".format(srv_root, data_name)
basset_root = "{0}/{1}/Naive_basset_motifs".format(srv_root, data_name)
reportfile= "/mnt/lab_data/kundaje/msharmin/mouse_hem/with_tfd/full_mouse50/filtering samples_MS2.xlsx"
sheetname = "filter23"
In [2]:
from matlas.dlutils import load_deeplift_data
from matlas.fimo_hits import (dichotomize_chipseq_for_gold_standard,
assign_significance_on_fimo_hits,
retain_fimo_hits_within_nb,
optimize_performance_info,
check_performance_competition,
collect_motif_df,
plot_various_instance_scores,
decide_optimum_threshold,
collect_distance_info,
compute_distances_from_peaks,
plot_positional_information
)
def prepare_verification_data_per_cell(
task_idx=273,
tf_biosample = "MEL cell line",
the_k=0.05627951420283896,
atac_root="/mnt/lab_data2/msharmin/oc-atlas/mtf/atac_specific_hits",
instance_root="/mnt/lab_data2/msharmin/oc-atlas/mtf/atac_various_hits",
deeplift_root="/mnt/lab_data2/msharmin/oc-atlas/deeplift2019_h5",
background_files_root="/mnt/lab_data2/msharmin/oc-atlas/mtf/t_stats"
):
tf_targets = sorted(set(CHIP_PEAK_BEDS[tf_biosample].keys()).intersection(set(CHIP_TO_MOTIF.keys())))
motifs = [CHIP_TO_MOTIF[tf] for tf in tf_targets]
# divide chipseq peaks into gold_standard and nongold for each tf under a cell type
fimodir = "{0}/task_{1}".format(atac_root, task_idx)
atacpeaks = "{}/peaks.bed".format(fimodir)
#print(fimodir)
#print(atacpeaks)
for tf, motif in zip(tf_targets, motifs):
print(tf_biosample, tf, motif)
chippeaks = "{0}/{1}".format(CHIP_PEAK_BED_ROOT, CHIP_PEAK_BEDS[tf_biosample][tf])
instancedir="{}/task_{}/{}".format(instance_root, task_idx, motif)
#print(chippeaks)
#print(instancedir)
#print()
best_key = decide_optimum_threshold(instancedir)
distances = collect_distance_info(best_key=best_key, instancedir=instancedir,
null_key='deeplift_cdf', key_to_split='deeplift_cdf',
fimobed="{}/{}.bed.gz".format(fimodir, motif),
verbose=False
)
plot_positional_information(distances, motif)
# compute_joint_chip_matrix and plot_joint_chip_matrix
return None
from ast import literal_eval
def prepare_verification_data():
total_zeros = 0
for key, tf_biosamples in OC_TO_TF_SAMPLES.items():
task_idx, sample_id, oc_sample = literal_eval(key)
print(task_idx, oc_sample, tf_biosamples, len(tf_biosamples))
total_zeros = prepare_verification_data_per_cell(task_idx=task_idx,
tf_biosample=tf_biosamples[0],
the_k=PER_BASE_BACKGROUND_IMPORTANCES[task_idx]
)
#break
return None
prepare_verification_data()
In [ ]:
# def prepare_verification_data():
# total_zeros = 0
# for key, tf_biosamples in OC_TO_TF_SAMPLES.items():
# task_idx, sample_id, oc_sample = literal_eval(key)
# print(task_idx, oc_sample, tf_biosamples, len(tf_biosamples))
# total_zeros = prepare_verification_data_per_cell(task_idx=task_idx,
# tf_biosample=tf_biosamples[0],
# the_k=PER_BASE_BACKGROUND_IMPORTANCES[task_idx]
# )
# #break
# return None
# prepare_verification_data()
In [12]:
from matlas.fimo_hits import compute_distances_from_peaks
#mdf = compute_distances_from_peaks()
mdf
Out[12]:
In [5]:
# determine the threshold and send it to both display as params
instancedir = "/mnt/lab_data2/msharmin/oc-atlas/mtf/atac_various_hits/task_140/Ctcf"
tmp = np.load("{}/cdf_cdf.npz".format(instancedir), allow_pickle=True)
dl_performdict, nondl_performdict, performdict = tmp.f.a.tolist(), tmp.f.b.tolist(), tmp.f.c.tolist()
performdict
Out[5]:
In [13]:
precisions = [(v)['Precision'].values[0] for v in dl_performdict.values()]
df = pd.DataFrame({'keys': list(dl_performdict.keys()),
'Precision': [(v)['Precision'].values[0] for v in dl_performdict.values()],
'Recall': [(v)['Recall'].values[0] for v in dl_performdict.values()],
'auPRC_scored': [(v)['auPRC_scored'].values[0] for v in dl_performdict.values()],
'F1-score': [(v)['F1-score'].values[0] for v in dl_performdict.values()],
'Balanced Accuracy': [(v)['Balanced Accuracy'].values[0] for v in dl_performdict.values()]
})
nondf = pd.DataFrame({'keys': list(nondl_performdict.keys()),
'Precision': [(v)['Precision'].values[0] for v in nondl_performdict.values()],
'Recall': [(v)['Recall'].values[0] for v in nondl_performdict.values()],
'auPRC_scored': [(v)['auPRC_scored'].values[0] for v in nondl_performdict.values()],
'F1-score': [(v)['F1-score'].values[0] for v in nondl_performdict.values()],
'Balanced Accuracy': [(v)['Balanced Accuracy'].values[0] for v in nondl_performdict.values()]
})
precision_cond = df['Precision']>=performdict['Fimo50']['Precision'].values[0]
df = df[precision_cond]
nondf = nondf[precision_cond]
recall_cond = df['Recall']>=nondf['Recall']
df = df[recall_cond]
nondf = nondf[recall_cond]
pd.concat([df, nondf], axis=1)
recall_cond
Out[13]:
In [36]:
from scipy.stats import *
from scipy import stats
from matlas.fimo_hits import collect_motif_df, PEAK_TYPES
from sklearn.metrics import precision_recall_curve
from sklearn.metrics import f1_score
from sklearn.metrics import auc
from matplotlib import pyplot as plt
from matlas.fimo_hits import compute_instance_calling_performance
perform_df, truth_df = compute_instance_calling_performance(
bedname="Dl_k50.0.bed",
key_to_optimize="deeplift_cdf",
instancedir="/mnt/lab_data2/msharmin/oc-atlas/mtf/atac_various_hits/task_224/Gata4",
atac_bed="/mnt/lab_data2/msharmin/oc-atlas/mtf/atac_specific_hits/task_224/peaks.bed",
verbose=True
)
# fimo_df, truth_df = compute_peak_calling_performance(key_to_optimize='mean_deeplift_0')
# fimo_df.T
In [19]:
dlk = fimo_df[fimo_df['deeplift_cdf']>0.9]
nondlk = fimo_df[fimo_df['deeplift_cdf']<=0.9]
from matplotlib import pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(16,6))
ax_pdf = axes[0]
ax_cdf = axes[1]
key_to_show = 'sum_deeplift_0'
ax_pdf.hist(dlk[key_to_show].values, density=True, histtype='stepfilled',
alpha=0.6, bins=60, label='high')
# plot the scores for an example sequence of a as well
ax_pdf.hist(nondlk[key_to_show].values, density=True, histtype='stepfilled',
alpha=0.6, bins=60, label='low')
ax_pdf.legend(loc='best', frameon=False)
ax_pdf.legend(loc='best', frameon=False)
ax_cdf.hist(dlk[key_to_show].values, cumulative=True, density=True, histtype='bar',
alpha=0.6, bins=60, label='high')
# plot the scores for an example sequence of a as well
ax_cdf.hist(nondlk[key_to_show].values, cumulative=True, density=True, histtype='bar',
alpha=0.6, bins=60, label='low')
ax_cdf.legend(loc='best', frameon=False)
plt.show()
In [ ]:
In [ ]:
In [ ]:
# revise performance measurements with changing the probability metric
# revise fimo50
# revise positional distribution and cdfs for a bed
# zcat chippeaks.bed.gz | sort -k1,1 -k2,2n | bedtools intersect -a ../../../atac_specific_hits/task_140/peaks.bed -b stdin -wb -wa -v > nongold3.bed
# zcat chippeaks.bed.gz | sort -k1,1 -k2,2n | bedtools intersect -a stdin -b ../../../atac_specific_hits/task_140/peaks.bed -wb -wa > gold.bed
# zcat /mnt/lab_data2/msharmin/tf-mouse-atlas/chip_idrs/ENCFF806PDR.bed.gz | sort -k1,1 -k2,2n | bedtools intersect -a stdin -b atac_specific_hits/task_140/peaks.bed -wb > atac_various_hits/task_140/Ctcf/gold.bed
In [117]:
b = np.sum(null_scores, axis=2)
a = np.sum(deeplift_scores, axis=2)
a.shape, b.shape
Out[117]:
In [118]:
np.argmin(cdf_scores[:,5]), np.min(cdf_scores[:,5])
Out[118]:
In [150]:
i = [0, 1]
c = b[0].flatten()
print(c.shape)
shape, mean, stdev = t.fit(c)
kstat, kval = stats.kstest(c, 't', args=[shape, mean, stdev])
# mean, stdev = norm.fit(c)
# kstat, kval = stats.kstest(c, 'norm', args=[ mean, stdev])
# shape, mean, stdev = gamma.fit(c)
# kstat, kval = stats.kstest(c, 'gamma', args=[shape, mean, stdev])
print(kstat, kval)
In [191]:
from matlas.fimo_hits import plot_fitted_dist_for_importance
plot_fitted_dist_for_importance(c, a[0], shape, mean, stdev)
In [11]:
# scores_ = [np.sum(j) for j in scores[0][300:700]]
# shape, mean, stdev = t.fit(scores_)
# kstat, kval = stats.kstest(scores_, 't', args=[shape, mean, stdev])
# print('kstat:', kstat, ' kval:', kval)
from scipy import stats
from scipy.stats import *
shape, mean, stdev = t.fit(b)
print(shape, mean, stdev)
kstat, kval = stats.kstest(b, 't', args=[shape, mean, stdev])
print('t sitribution', kstat, kval)
kstat, kval = stats.kstest(b, 'norm')
print('normal sitribution', kstat, kval)
In [ ]:
In [28]:
#get the scores for all the windows, take per_base_imp distribution
# nonagnostic_out = "/mnt/lab_data/kundaje/msharmin/mouse_hem/mtf/specific_hits/task_224/Ctcf.bed.gz"
# df = pd.read_csv(nonagnostic_out, sep="\t")
# selected_peaks = np.sort(np.array(list(set(df['seqno'].values))))
# selected_null_scores = hyp_scores[selected_peaks,:,:]
# selected_null_scores.shape
# atacdf = pd.read_csv("/mnt/lab_data/kundaje/msharmin/mouse_hem/mtf/specific_hits/task_224/peaks.bed", header=None, sep="\t")
# #selected_atacs = atacdf.loc[selected_peaks,:]
# #selected_atacs.to_csv("/mnt/lab_data/kundaje/msharmin/mouse_hem/mtf/specific_hits/task_224/peaks_Ctcf.bed",
# #header=None, index=None, sep="\t")
In [ ]:
In [42]:
# from matlas.render_report import render_cell_tf_reports_in_parallel
# render_cell_tf_reports_in_parallel()
In [2]:
# call dichotomize_chipseq_for_gold_standard for each tf under MEL
# call get_deeplift_passed_fimo_hits for each tf under MEL
# call get_peaks_within_nb for each tf under MEL
# for a tf: collect_distance_info, plot_distances_from_peaks
# for a tf: collect_performance_info, plot_performance_info
# for a tf: compute_joint_chip_matrix, plot_joint_chip_heatmap
# for a tf: load the png file
In [2]:
from matlas.fimo_hits import combine_motif_dfs
logdir = "/mnt/lab_data/kundaje/msharmin/mouse_hem/mtf"
outdir = "/mnt/lab_data/kundaje/msharmin/mouse_hem/mtf/specific_hits/task_{}".format(273)
nonagnostic_out = "{0}/{1}.bed.gz".format(outdir, "Adnp2")
df_dict = combine_motif_dfs(logdir, outdir,
combine_func="collect_motif_df",
)
df = pd.concat([df_dict[key] for key in df_dict.keys()])
#nonz_dfs = combine_motif_dfs(logdir, outdir, zscored=False)
In [93]:
from matlas.fimo_hits import plot_motif_scores
motifs = ['Spi1', 'Ctcf', 'Gata1']
turn = 0
fig, axes = None, None
for key in motifs:
fig, axes = plot_motif_scores(df_dict[key], title=key, turn=turn, axes=axes, fig=fig, kind='scatter')
turn = (turn+1)%4
In [87]:
import seaborn as sns
combined_df = pd.concat([df_dict[key] for key in df_dict.keys()])
sns.jointplot(x="avg_deeplift", y="fimo", data=combined_df, s=1, marginal_kws=dict(bins=100));
In [19]:
df_dict['Spi1'][:5]
Out[19]:
In [3]:
from matlas.render_report import render_test_plots
render_test_plots(nbname="FimoVsDeepliftReport")
In [71]:
gencode_file = '/mnt/data/annotations/by_release/mm10/GENCODE_ann/gencode.vM6.annotation.gtf'
gencode = pd.read_table(gencode_file, sep='\t', header=None, comment='#')
genes = gencode[(gencode[0]!='chrM') & (gencode[2]=='gene')&(np.array(['protein_coding' in i for i in gencode[8].values]))]
gene_df = genes.loc[:, [0,3,4,8,5,6]]
gene_df = gene_df.sort_values([0,3])
gene_df.to_csv("{0}/genecode.bed.gz".format(modisco_dir), header=False, index=False, compression='gzip', sep="\t")
genes[:5]
Out[71]:
In [ ]: