<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Blastocyst Development in Mice
Step1:  Next we load in the data. We've provided a convenience function for loading in the data with GPy. It is loaded in as a pandas DataFrame. This allows us to summarize it with the describe attribute.
Step2:  Principal Component Analysis
Now we follow Guo et al in performing PCA on the data. Rather than calling a 'PCA routine', here we break the algorithm down into its steps
Step3:  Here we will define a little routine, which puts the legend on the right side of the plot, so that the legend does not overwright the data plot
Step4:  PCA Result
Now, using the helper function we can plot the results with appropriate labels.
Step5:  GP-LVM on the Data
<a href="http
Step6:  Bayesian GP-LVM
Finally we show the new code that uses the Bayesian GP-LVM to fit the data. This means we can automatically determine the dimensionality of the model whilst fitting a non-linear dimensionality reduction. The approximations we use also mean that it is faster than the original GP-LVM.
Step7:  This gives a really nice result. Broadly speaking two latent dimensions dominate the representation. When we visualize using these two dimensions we can see the entire cell phylogeny laid out nicely in the two dimensions. Additionally we can see the missclassification of the some cells, using the 'standard' approach of repeated k-means clustering and PCA on sub clustered (This was used to get the sample colors of the 64 cellstage).
Other Nonlinear Approaches
How do other non-linear approaches compare with the GP-LVM? We can use scikit-learn to try a few. First we consider Isomap, which often performs really well on data.
Isomap
Isomap first builds a neighbourhood graph, and then uses distances along this graph to approximate the geodesic distance between points. These distances are then visualized by performing classical multidimensional scaling (which involves computing the eigendecomposition of the centred distance matrix). As the neighborhood size is increased to match the data, principal component analysis is recovered (or strictly speaking, principal coordinate analysis). The fewer the neighbors, the more 'non-linear' the isomap embeddings are.
Step8:  Locally Linear Embedding
Next we try locally linear embedding. In locally linear embedding a neighborhood is also computed. Each point is then reconstructed by it's neighbors using a linear weighting. This implies a locally linear patch is being fitted to the data in that region. These patches are assimilated into a large $n\times n$ matrix and a lower dimensional data set which reflects the same relationships is then sought. Quite a large number of neighbours needs to be selected for the data to not collapse in on itself. When a large number of neighbours is selected the embedding is more linear and begins to look like PCA. However, the algorithm does not converge to PCA in the limit as the number of neighbors approaches $n$.
Step9:  t-SNE (t-distributed Stochastic Neighbour Embedding)
t-SNE does a great job in unravelling the structure of the data itself. It also shows the misslabelling of the later data quite clearly. It is not able to put the right order in the TE side, though.
<ASSISTANT_TASK:>
Python Code:
import pods, GPy, itertools
%matplotlib inline
from matplotlib import pyplot as plt
Explanation: Blastocyst Development in Mice: Single Cell TaqMan Arrays
presented at the EBI BioPreDyn Course 'The Systems Biology Modelling Cycle'
Max Zwiessele, Oliver Stegle, Neil Lawrence 12th May 2014 University of Sheffield and EBI.
In this notebook we follow  <a href="http://bioinformatics.oxfordjournals.org/content/28/18/i626.full.pdf">Buettner and Theis (2012)</a> and use the GP-LVM to analyze some single cell data from <a href='http://dx.doi.org/10.1016/j.devcel.2010.02.012'>Guo et al (2010)</a>. They performed qPCR TaqMan array on single cells from the developing blastocyst in mouse. The data is taken from the early stages of development when the Blastocyst is forming. At the 32 cell stage the data is already separated into the trophectoderm (TE) which goes onto form the placenta and the inner cellular mass (ICM). The ICM further differentiates into the epiblast (EPI)---which gives rise to the endoderm, mesoderm and ectoderm---and the primitive endoderm (PE) which develops into the amniotic sack. Guo et al selected 48 genes for expression measurement. They labelled the resulting cells and their labels are included as an aide to visualization.
In the paper, they first visualized their data using principal component analysis. In the first two principal components this fails to separate the domains. This is perhaps because the principal components are dominated by the variation in the 64 cell systems. This in turn may be because there are more cells from the data set in that regime, or additionally it could be that the natural variation across the 64 cell systems is greater. Both are probable causes of the dominance of these cells in the first two principal components. The first thing we do is plot the principal coordinates of the cells on the first two dominant directions. 
This Notebook
In this notebook we will perform PCA on the original data, showing that the different regimes do not separate. Then, we follow <a href="http://bioinformatics.oxfordjournals.org/content/28/18/i626.full.pdf">Buettner and Theis (2012)</a> in applying the GP-LVM to the data. There is a slight pathology in the result, one which they fixed by using priors that were dependent on the developmental stage. We then show how the Bayesian GP-LVM doesn't exhibit those pathologies and gives nice results that seems to show the lineage of the cells. 
First we perform some set up.
End of explanation
s = pods.datasets.singlecell()
Ydf = s['Y']
Y = Ydf.values
labels = s['labels']
marker = '<>^vsd'
Ydf.describe()
Explanation: Next we load in the data. We've provided a convenience function for loading in the data with GPy. It is loaded in as a pandas DataFrame. This allows us to summarize it with the describe attribute.
End of explanation
import numpy as np
# obtain a centred version of data.
centredY = Y - Y.mean()
# compute inner product matrix
C = np.dot(centredY,centredY.T)
# perform eigendecomposition
V, U = np.linalg.eig(C)
# sort eigenvalues and vectors according to size
ind = V.argsort()  
ev = V[ind[::-1]]
U = U[:, ind[::-1]]
Explanation: Principal Component Analysis
Now we follow Guo et al in performing PCA on the data. Rather than calling a 'PCA routine', here we break the algorithm down into its steps: compute the data covariance, compute the eigenvalues and eigenvectors and sort according to magnitude of eigenvalue. Because we want to visualize the data, we've chose to compute the eigenvectors of the inner product matrix rather than the covariance matrix. This allows us to plot the eigenvalues directly. However, this is less efficient (in this case because the number of genes is smaller than the number of data) than computing the eigendecomposition of the covariance matrix.
End of explanation
def legend_to_the_right(ax):
    box = ax.get_position()
    ax.set_position([box.x0, box.y0, box.width, box.height])
    _ = ax.legend(loc='center left', bbox_to_anchor=(1, 0.5), numpoints=1)
    
def plot_latent(ax, x, y, marker, labels):
    GPy.plotting.Tango.reset()
    # make sure labels are in order of input:
    ulabels = []
    for lab in labels:
        if not lab in ulabels:
            ulabels.append(lab)
    for i, [label, m] in enumerate(zip(ulabels, itertools.cycle(marker))):
        symbol = marker[i % len(marker)]
        ind = labels == label
        ax.plot(x[ind], y[ind], marker=symbol, 
                   markerfacecolor=GPy.plotting.Tango.nextMedium(), 
                linestyle='',   
                label=label, mew=.2, alpha=1)
Explanation: Here we will define a little routine, which puts the legend on the right side of the plot, so that the legend does not overwright the data plot
End of explanation
fig, [ax1, ax2] = plt.subplots(1,2,figsize=(13,4))
fracs = (ev.real/ev.real.sum())[:30] # Take only the first 30 eigenvalues
ax1.bar(np.arange(len(fracs)), fracs)
plot_latent(ax2, U[:,0], U[:,1], marker, labels)
legend_to_the_right(ax2)
Explanation: PCA Result
Now, using the helper function we can plot the results with appropriate labels.
End of explanation
m = GPy.models.GPLVM(Y, 2, kernel=GPy.kern.RBF(2, ARD=1)+GPy.kern.Bias(2))
m.optimize(messages=1)
fig, [ax1, ax2] = plt.subplots(1,2,figsize=(9,4))
m.kern.plot_ARD(ax=ax1)
m.plot_latent(labels=labels, marker=marker, legend=False, ax=ax2)
_ = legend_to_the_right(ax2)
fig, [ax1, ax2] = plt.subplots(1,2,figsize=(9,4))
m.plot_latent(labels=labels, marker=marker, legend=False, ax=ax1)
m.plot_magnification(labels=labels, marker=marker, legend=False, ax=ax2)
_ = legend_to_the_right(ax2)
Explanation: GP-LVM on the Data
<a href="http://bioinformatics.oxfordjournals.org/content/28/18/i626.long">Buettner and Theis</a> applied a Gaussian process latent variable model to the data, but with a modified prior to ensure that small differences between cells at the same differential stage were preserved. Here we apply a standard GP-LVM (no modified prior) to the data.
End of explanation
m = GPy.models.BayesianGPLVM(Y, 5, num_inducing=15, kernel=GPy.kern.RBF(5,ARD=1)+GPy.kern.Bias(5))
m.optimize('bfgs', max_iters=5e3, messages=1)
fig, [ax1, ax2] = plt.subplots(1,2,figsize=(9,4))
m.kern.plot_ARD(ax=ax1)
m.plot_latent(labels=s['labels'], marker='<>^vsd', legend=False, ax=ax2)
legend_to_the_right(ax2)
Explanation: Bayesian GP-LVM
Finally we show the new code that uses the Bayesian GP-LVM to fit the data. This means we can automatically determine the dimensionality of the model whilst fitting a non-linear dimensionality reduction. The approximations we use also mean that it is faster than the original GP-LVM.
End of explanation
n_neighbors = 40
import sklearn.manifold
m = sklearn.manifold.Isomap(n_neighbors=n_neighbors, n_components=2)
X = m.fit_transform(Ydf)
fig, ax = plt.subplots(1)
plot_latent(ax, X[:, 0], X[:, 1], marker, labels)
_ = legend_to_the_right(ax)
Explanation: This gives a really nice result. Broadly speaking two latent dimensions dominate the representation. When we visualize using these two dimensions we can see the entire cell phylogeny laid out nicely in the two dimensions. Additionally we can see the missclassification of the some cells, using the 'standard' approach of repeated k-means clustering and PCA on sub clustered (This was used to get the sample colors of the 64 cellstage).
Other Nonlinear Approaches
How do other non-linear approaches compare with the GP-LVM? We can use scikit-learn to try a few. First we consider Isomap, which often performs really well on data.
Isomap
Isomap first builds a neighbourhood graph, and then uses distances along this graph to approximate the geodesic distance between points. These distances are then visualized by performing classical multidimensional scaling (which involves computing the eigendecomposition of the centred distance matrix). As the neighborhood size is increased to match the data, principal component analysis is recovered (or strictly speaking, principal coordinate analysis). The fewer the neighbors, the more 'non-linear' the isomap embeddings are.
End of explanation
n_neighbors = 30
m = sklearn.manifold.LocallyLinearEmbedding(n_neighbors=n_neighbors, n_components=2)
X = m.fit_transform(Y)
fig, ax = plt.subplots(1)
plot_latent(ax, X[:, 0], X[:, 1], marker, labels)
_ = legend_to_the_right(ax)
Explanation: Locally Linear Embedding
Next we try locally linear embedding. In locally linear embedding a neighborhood is also computed. Each point is then reconstructed by it's neighbors using a linear weighting. This implies a locally linear patch is being fitted to the data in that region. These patches are assimilated into a large $n\times n$ matrix and a lower dimensional data set which reflects the same relationships is then sought. Quite a large number of neighbours needs to be selected for the data to not collapse in on itself. When a large number of neighbours is selected the embedding is more linear and begins to look like PCA. However, the algorithm does not converge to PCA in the limit as the number of neighbors approaches $n$.
End of explanation
m = sklearn.manifold.TSNE(n_components=2)
X = m.fit_transform(Y)
fig, ax = plt.subplots(1)
plot_latent(ax, X[:, 0], X[:, 1], marker, labels)
_ = legend_to_the_right(ax)
Explanation: t-SNE (t-distributed Stochastic Neighbour Embedding)
t-SNE does a great job in unravelling the structure of the data itself. It also shows the misslabelling of the later data quite clearly. It is not able to put the right order in the TE side, though.
End of explanation
<END_TASK>