<jupyter_start><jupyter_text>Example WorkflowIn this notebook we will look at an example workflow of running topslam on a dataset Y. In order to show the expected behaviours and outcomes, we will use a simulated differentiation pattern. The simulation is included in the package and can be accessed as follows:<jupyter_code>from topslam.simulation import qpcr_simulation
seed_differentiation = 5001
seed_gene_expression = 0
Xsim, simulate_new, t, c, labels, seed = qpcr_simulation(seed=seed_differentiation)
np.random.seed(seed_gene_expression)
Y = simulate_new()<jupyter_output><empty_output><jupyter_text>Dimensionality reduction initializationIn order to fit the Baysian GPLVM, we will use five other dimensionality reduction techniques as initialization. This means, we present the BaysianGPLVM with the solutions generated from these techniques as a **starting point**. BayesianGPLVM will then proceed from there and find the most likeliy fit and reduced space dimensionality. (See ARD parameters later)<jupyter_code>from sklearn.manifold import TSNE, LocallyLinearEmbedding, SpectralEmbedding, Isomap
from sklearn.decomposition import FastICA, PCA

methods = {'t-SNE':TSNE(n_components=2),
           'PCA':PCA(n_components=2), 
           'Spectral': SpectralEmbedding(n_components=2, n_neighbors=10), 
           'Isomap': Isomap(n_components=2, n_neighbors=10), 
           'ICA': FastICA(n_components=2)
           }
from topslam.optimization import run_methods

X_init, dims = run_methods(Y, methods)<jupyter_output><empty_output><jupyter_text>BayesianGPLVMWe run BayesianGPLVM using the optimization routines given in the optimization package of topslam. You can run the BayesianGPLVM yourself for more control (see for example https://github.com/SheffieldML/notebook/blob/master/compbio/SingleCellDataWithGPy.ipynb). We do not expect any linear relation in the data and will therefore not model a linear dimenion:<jupyter_code>from topslam.optimization import create_model, optimize_model
m = create_model(Y, X_init, linear_dims=0)
m.optimize(messages=1, max_iters=5e3)<jupyter_output><empty_output><jupyter_text>ARD: Discovery of significant dimensionsWe can assign the number of dimensions to be used in the reduced space by the ARD parameters of the kernel.<jupyter_code>m.kern.plot_ARD()
print("Most significant input dimensions: {}".format(m.kern.get_most_significant_input_dimensions()[:2]))<jupyter_output>Most significant input dimensions: (8, 6)<jupyter_text>Lower Dimensional RepresentationSo now we know which dimensions are important and we can look at those, by plotting them using the labels generated by the simulation:<jupyter_code>m.plot_scatter(labels=labels)<jupyter_output><empty_output><jupyter_text>Landscape extractionEverything before this stage was initialization and running of the model. We can now use the learnt model to extract the probabilistic landscape from this model, and plot the extracted timeline. We will also compare it to the other dimensionality reduction techniques and their extracted timelines (for ICA and t-SNE).<jupyter_code>from topslam import ManifoldCorrectionTree
m_topslam = ManifoldCorrectionTree(m)<jupyter_output><empty_output><jupyter_text>We can plot the extracted landscape and a graph representing the minimal spanning tree along which we extract the time. We simulated the time for the samples in order, so the start for the extracted time is 0. The startin cell plays a crucial role in the extraction of the timeline and should always be supplied.<jupyter_code>ax = m_topslam.plot_waddington_landscape()
m_topslam.plot_time_graph(labels=labels, start=0, ax=ax)
from topslam.plotting import plot_comparison
with sns.axes_style('white'):
    ax = plot_comparison(m_topslam, X_init, dims, labels, np.unique(labels), 0)<jupyter_output><empty_output><jupyter_text>PseudotimeWe can now extract the pseudotime and compare it to simulated times. We will compare topslam to ICA and t-SNE extracted landscapes using pearson correlation between extracted and simulated times:<jupyter_code>from scipy.stats import linregress
from scipy.spatial.distance import pdist, squareform
from topslam.simulation.graph_extraction import extract_manifold_distances_mst
pt_topslam = m_topslam.get_pseudo_time(start=0)
D_ica, mst_ica = extract_manifold_distances_mst(squareform(pdist(X_init[:, dims['ICA']])))
pt_ica = D_ica[0] # only take the distances from the start (which is 0 in this case)
D_tsne, mst_tsne = extract_manifold_distances_mst(squareform(pdist(X_init[:, dims['t-SNE']])))
pt_tsne = D_tsne[0] # only take the distances from the start (which is 0 in this case)
result_df = pd.DataFrame(columns=['slope', 'intercept', 'r_value', 'p_value', 'std_err'])
result_df.loc['topslam'] = linregress(pt_topslam, t.flat)[:]
result_df.loc['ICA'] = linregress(pt_ica, t.flat)[:]
result_df.loc['t-SNE'] = linregress(pt_tsne, t.flat)
result_df
time_df = pd.DataFrame(t, columns=['simulated'])
time_df['topslam'] = pt_topslam
time_df['ICA'] = pt_ica
time_df['t-SNE'] = pt_tsne
sns.jointplot('simulated', 'topslam', data=time_df, kind='reg')
sns.jointplot('simulated', 'ICA', data=time_df, kind='reg')
sns.jointplot('simulated', 't-SNE', data=time_df, kind='reg')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08492366412213741, 'fraction_numerical': 0.005534351145038168, 'mean_word_length': 6.881203007518797, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 32, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '646416', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1682, 'n_tokens_neox': 1595, 'n_words': 507, 'path': 'notebooks/ExampleWorkflow.ipynb', 'repo_name': 'mzwiessele/cellSLAM', 'size': '5240'}