<jupyter_start><jupyter_text>Hidden Markov Tree<jupyter_code>%load_ext rpy2.ipython
%matplotlib inline
import sys
import h5py
import warnings
import pickle
import numpy as np
import pandas as pd
import stream as st
import matplotlib.pyplot as plt
from anndata import AnnData
from sklearn.metrics import confusion_matrix

import utility_functions as uf
import matrices as mat
import hmt_expectation_maximization as hmtem

warnings.filterwarnings("ignore", category=RuntimeWarning)<jupyter_output><empty_output><jupyter_text>Load dataset<jupyter_code>data_file = 'datafile.hdf5'
dev_abb = ['PMY', 'MY', 'MMY', 'BNE', 'SNE', 'BL', 'BA', 'EO', 'LY', 'PC']
dev_map = {name: i for i, name in enumerate(dev_abb)}

f = h5py.File(data_file, 'r')
X = np.array(f['train']['X'])
Y = np.array(f['train']['y']).flatten()
ID = f.attrs['filenames']
f.close()

classes = np.array([dev_abb[y] for y in Y])
print(X.shape, Y.shape)<jupyter_output>(1821, 50) (1821,)<jupyter_text>The dataset consists of 1821 feature vectors of the ten classes promyelocyte (PMY), myelocyte (MY), metamyelocyte (MMY), band neutrophil (BNE), segmented neutrophil (SNE), blast (BL), basophil (BA), eosinophil (EO), lymphocyte (LY), and plasma cells (PC).  Calculate pseudotime using STREAM https://github.com/pinellolab/STREAM<jupyter_code>X_adata = AnnData(X)
X_adata.uns['label_color'] = {'PMY': '#4d76ff', 'BL': '#f14cc1', 'MY': '#ff7c00','BA': '#a3a3a3', 'MMY': '#1ac938', 
                              'EO': '#ffc400', 'BNE': '#e8000b', 'LY': '#33e0ff', 'SNE': '#a862ea', 'PC': '#FFDEAD'}

col = ['#4d76ff', '#ff7c00', '#1ac938', '#e8000b', '#a862ea', '#f14cc1', 
       '#a3a3a3','#ffc400', '#33e0ff', '#FFDEAD']
colors = np.repeat('#color#', 1821)
for i in range(len(dev_abb)):
    colors[Y == i] = col[i]

X_adata.obs['label_color'] = colors
X_adata.obs['label'] = classes
X_adata.uns['experiment'] = 'rna-seq'
st.log_transform(X_adata)
st.select_variable_genes(X_adata, loess_frac=0.1)
st.dimension_reduction(X_adata, n_components=3)
st.plot_visualization_2D(X_adata, fig_legend_ncol=5)
st.seed_elastic_principal_graph(X_adata, clustering='ap')
st.elastic_principal_graph(X_adata, epg_alpha=0.007)
st.subwaymap_plot(X_adata, percentile_dist=99,root='S7',fig_legend_ncol=5,preference=['S6', 'S4']) 
uf.hmt_boxplot_ptime(X_adata, dev_abb)
# Information about the tree has been stored in the variables X_obs and tree_obs.
# tree_obs is a dictionary containing information about the parents, children, root node, and leaf nodes.
X_obs = pd.read_pickle('X_obs')
pkl_file = open('tree_obs', 'rb')
tree_obs = pickle.load(pkl_file)
pkl_file.close()
X_obs.head()
# Pseudotime difference between a node and the parent of the node
pt_diff = []
for i in range(len(X_obs.index)):
        parent = tree_obs['parents'][str(i)]
        if parent != None:
            pt_diff.append(X_obs.loc[str(i)]['S7_pseudotime'] - X_obs.loc[parent]['S7_pseudotime'])
pt_diff = np.array(pt_diff)<jupyter_output><empty_output><jupyter_text>EM algorithm<jupyter_code># Rescale the pseudotime differences.
pt_diff_max = 1
pt_diff_rescaled = pt_diff_max * (pt_diff - np.max(pt_diff)) / (np.max(pt_diff) - 
                                                                np.min(pt_diff)) + pt_diff_max
# Rescaled lmbda_fix
lmbda_fix_re = []
for abb in dev_abb:
    ptime_class = np.sort(X_obs.loc[X_obs['label'] == abb]['S7_pseudotime'])
    diff = ptime_class[1:] - ptime_class[:-1]
    diff_rescaled = pt_diff_max * (diff - np.max(diff)) / (np.max(diff) - np.min(diff)) + pt_diff_max
    lmbda_fix_re.append(np.log(2)/np.median(diff_rescaled))

ptime_class = np.sort(X_obs.loc[X_obs['label'] == 'LY']['S7_pseudotime'])
diff = ptime_class[1:] - ptime_class[:-1]
diff_rescaled = pt_diff_max * (diff - np.max(diff)) / (np.max(diff) - np.min(diff)) + pt_diff_max
lmbda_fix_re[8] = 1/np.mean(diff_rescaled)
lmbda_fix_re = np.array(lmbda_fix_re)
p, lmbda, p_hist, lmbda_hist, loglik_hist = hmtem.hmt_EM_algorithm(lmbda_fix_re, 100, X_obs, tree_obs, 
                                                                   pt_diff_rescaled, dev_map)

plt.figure()
plt.plot(range(len(loglik_hist)) , loglik_hist)
plt.xlabel('Iteration', fontsize=14)
plt.ylabel('Log-likelihood', fontsize=14)
plt.show()<jupyter_output>Iteration 9<jupyter_text>Results<jupyter_code>tm = mat.hmt_transition_matrix(p, lmbda, lmbda_fix_re, X_obs, tree_obs)
M, X_obs = uf.hmt_viterbi(X_obs, tree_obs, tm, dev_abb, dev_map)
acc =  len(np.where(X_obs['label'] == X_obs['Z'])[0])/1821
print('Accuracy: ', acc)
branches = ['branch S7-S5', 'branch S6-S5', 'branch S5-S3', 'branch S4-S3', 
            'branch S1-S3', 'branch S2-S1', 'branch S0-S1']
transition_idx = dict()
for b in branches:
    branch = list(reversed(tree_obs[b]))
    parent = tree_obs['parents'][branch[0]]
    if parent != None:
        if X_obs.loc[parent].Z != X_obs.loc[branch[0]].Z:
            transition_idx[X_obs.loc[parent].Z + '-' + X_obs.loc[branch[0]].Z] = branch[0]
    Z = X_obs.loc[branch].Z
    for i in range(len(Z)-1):
        if Z[i] != Z[i+1]:
            transition_idx[Z[i] + '-' + Z[i+1]] = branch[i+1]
print('Transition indices: ', transition_idx)
X_obs.head()
transition_times = dict()
for key in transition_idx.keys():
    time_parent = X_obs.loc[tree_obs['parents'][transition_idx[key]]].S7_pseudotime
    time_child = X_obs.loc[transition_idx[key]].S7_pseudotime
    ttime = (time_child + time_parent) / 2
    transition_times[key] = ttime   
print('Transition times: ', transition_times)
labels = np.array(X_obs.label)
for abb in dev_abb:
    labels[np.where(labels == abb)] = dev_map[abb]
labels = labels.astype(int)

Z = np.array(X_obs.Z)
for abb in dev_abb:
    Z[np.where(Z == abb)] = dev_map[abb]
Z = Z.astype(int)

plt.figure(figsize=(8,8))
uf.plot_confusion_matrix(confusion_matrix(Z, labels), classes=dev_abb, normalize=True)
plt.title('Confusion matrix', fontsize=15)
plt.show()
X_adata.obs['label'] = np.array(X_obs.Z)
colors = np.repeat('#color#', 1821)
for i in range(len(dev_abb)):
    colors[Z == i] = col[i]
X_adata.obs['label_color'] = colors
# Classification according to the hidden Markov tree
st.subwaymap_plot(X_adata, percentile_dist=99,root='S7',fig_legend_ncol=5,preference=['S6', 'S4'])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16161616161616163, 'fraction_numerical': 0.025493025493025494, 'mean_word_length': 4.929657794676806, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 19, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '123012', 'license': 'Apache-2.0', 'n_tokens_mistral': 2568, 'n_tokens_neox': 2451, 'n_words': 470, 'path': 'demo/hmt.ipynb', 'repo_name': 'liu-yushan/label-consistency', 'size': '6237'}