In [1]:
import os
import sys
import pickle
import numpy as np
from math import exp
from scipy import stats
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt
from vizsequence.viz_sequence import plot_weights_given_ax
from scipy.special import softmax
import keras
import keras.losses
from keras.models import Model, Sequential, load_model
from keras import backend as K
import numpy.random as rng
import seaborn as sns
from collections import OrderedDict
from basepair.losses import twochannel_multinomial_nll
import modisco
import modisco.tfmodisco_workflow.workflow
import h5py
import modisco.util
from collections import Counter
from modisco.visualization import viz_sequence
import modisco.affinitymat.core
import modisco.cluster.phenograph.core
import modisco.cluster.phenograph.cluster
import modisco.cluster.core
import modisco.aggregator
%matplotlib inline
In [2]:
folder = "pho4_pbexo"
task_names = ["pho4"]
for task_id, task in enumerate(task_names):
tfname = task
print("\n\n"+tfname+" COUNTS\n\n")
counts_hdf5_results = h5py.File("data/modisco/"+folder+"/"+tfname+"_counts/results.hdf5","r")
activity_patterns = np.array(counts_hdf5_results['metaclustering_results']['attribute_vectors'])[
np.array(
[x[0] for x in sorted(
enumerate(counts_hdf5_results['metaclustering_results']['metacluster_indices']),
key=lambda x: x[1])])]
plt.figure()
sns.heatmap(activity_patterns, center=0)
plt.show()
# plt.savefig("data/modisco/"+folder+'/counts/activity_patterns.png')
metacluster_names = [
x.decode("utf-8") for x in
list(counts_hdf5_results["metaclustering_results"]
["all_metacluster_names"][:])]
all_patterns = []
for metacluster_name in metacluster_names:
print("\n"+metacluster_name+"\n")
metacluster_grp = (counts_hdf5_results["metacluster_idx_to_submetacluster_results"]
[metacluster_name])
all_pattern_names = [x.decode("utf-8") for x in
list(metacluster_grp["seqlets_to_patterns_result"]
["patterns"]["all_pattern_names"][:])]
for pattern_name in all_pattern_names:
all_patterns.append((metacluster_name, pattern_name))
pattern = metacluster_grp["seqlets_to_patterns_result"]["patterns"][pattern_name]
print("\n"+pattern_name+": "+str(len(pattern["seqlets_and_alnmts"]["seqlets"]))+" seqlets\n")
fig, ax = plt.subplots(4, 1, figsize=(30,20))
background = np.array([0.27, 0.23, 0.23, 0.27])
plot_weights_given_ax(ax[0], pattern[tfname+"_hypothetical_contribs"]["fwd"])
ax[0].set_title("hypothetical scores")
plot_weights_given_ax(ax[1], pattern[tfname+"_contrib_scores"]["fwd"])
ax[1].set_title("actual importance scores")
plot_weights_given_ax(ax[2], viz_sequence.ic_scale(np.array(pattern["sequence"]["fwd"]),
background=background))
ax[2].set_title("onehot fwd")
plot_weights_given_ax(ax[3], viz_sequence.ic_scale(np.array(pattern["sequence"]["rev"]),
background=background))
ax[3].set_title("onehot rev")
plt.show()
# fig.savefig("data/modisco/"+folder+"/counts_"+\
# str(metacluster_name)+'_'+str(pattern_name)+'_activity_'+\
# str(metacluster_grp["activity_pattern"][:])+'_seqlets_'+\
# str(len(pattern["seqlets_and_alnmts"]["seqlets"]))+'.png')
counts_hdf5_results.close()
print("\n\n"+tfname+" PROFILE\n\n")
profile_hdf5_results = h5py.File("data/modisco/"+folder+"/"+tfname+"_profile/results.hdf5","r")
activity_patterns = np.array(profile_hdf5_results['metaclustering_results']['attribute_vectors'])[
np.array(
[x[0] for x in sorted(
enumerate(profile_hdf5_results['metaclustering_results']['metacluster_indices']),
key=lambda x: x[1])])]
plt.figure()
sns.heatmap(activity_patterns, center=0)
plt.show()
# plt.savefig("data/modisco/"+folder+'/profile/activity_patterns.png')
metacluster_names = [
x.decode("utf-8") for x in
list(profile_hdf5_results["metaclustering_results"]
["all_metacluster_names"][:])]
all_patterns = []
for metacluster_name in metacluster_names:
print("\n"+metacluster_name+"\n")
metacluster_grp = (profile_hdf5_results["metacluster_idx_to_submetacluster_results"]
[metacluster_name])
all_pattern_names = [x.decode("utf-8") for x in
list(metacluster_grp["seqlets_to_patterns_result"]
["patterns"]["all_pattern_names"][:])]
for pattern_name in all_pattern_names:
all_patterns.append((metacluster_name, pattern_name))
pattern = metacluster_grp["seqlets_to_patterns_result"]["patterns"][pattern_name]
print("\n"+pattern_name+": "+str(len(pattern["seqlets_and_alnmts"]["seqlets"]))+" seqlets\n")
fig, ax = plt.subplots(4, 1, figsize=(30,20))
background = np.array([0.27, 0.23, 0.23, 0.27])
plot_weights_given_ax(ax[0], pattern[tfname+"_hypothetical_contribs"]["fwd"])
ax[0].set_title("hypothetical scores")
plot_weights_given_ax(ax[1], pattern[tfname+"_contrib_scores"]["fwd"])
ax[1].set_title("actual importance scores")
plot_weights_given_ax(ax[2], viz_sequence.ic_scale(np.array(pattern["sequence"]["fwd"]),
background=background))
ax[2].set_title("onehot fwd")
plot_weights_given_ax(ax[3], viz_sequence.ic_scale(np.array(pattern["sequence"]["rev"]),
background=background))
ax[3].set_title("onehot rev")
plt.show()
# fig.savefig("data/modisco/"+folder+"/profile_"+\
# str(metacluster_name)+'_'+str(pattern_name)+'_activity_'+\
# str(metacluster_grp["activity_pattern"][:])+'_seqlets_'+\
# str(len(pattern["seqlets_and_alnmts"]["seqlets"]))+'.png')
profile_hdf5_results.close()
In [3]:
print("DONE")