BACKGROUND
the brain consists of an enormous number of cells organized into structures  <cit> .
the structured organization of cells is the key to the functional efficiency of the brain  <cit> .
hence, a natural first step toward understanding the brain function would be to address basic research questions at the structure level.
how cells are organized into structures  <cit> ?
what are the functions of structures  <cit> ?
how the structures are connected to each other  <cit> ?
however, a fundamental difficulty of understanding brain functions at the structure level lies in that there is no universally agreed division of cells into structures  <cit> .
from a developmental perspective, the delicate organization of brain into structures is the consequence of stringent spatiotemporal patterning controlled by the molecular signals during development.
in this process, cells at different spatial locations read different morphogenetic positional signals produced by the graded distribution of signaling molecules.
these signals control the expression of a relatively small set of transcription factors, which in turn regulate the expression of a larger number of genes.
this sequential cascade of expression control ultimately leads to cell differentiation and the emergence of connections and functional properties  <cit> .
the discovery that certain marker genes are expressed in regionally restricted patterns in the developing brain has either led to the introduction of new structural boundaries or made it possible to re-define existing boundaries at a higher resolution  <cit> .
currently, studies on the molecular-structural associations are beginning to reveal how the spatiotemporal gene expression patterns are related to cellular differentiation and structural development  <cit> .
in this article, we study the relationship between brain anatomy and spatiotemporal gene expression patterns in the developing mouse brain.
this global study of developing neuroanatomy is made possible by the high-resolution, three-dimensional  gene expression patterns provided by the allen brain atlas   <cit> .
as part of the aba, the allen developing mouse brain atlas provides spatiotemporal in situ hybridization  gene expression pattern images across four embryonic and three postnatal developmental ages  <cit> , yielding effectively a four-dimensional brain atlas.
to establish a common coordinate framework for analyzing the ish data, the ish image series are aligned to the allen developing mouse brain reference atlas.
this enables the global, computational study of the spatiotemporal gene expression patterns of many genes and comparison of the results with developmental anatomy.
to enable visual explorations of the gene expression patterns and correlate the results with classically defined neuroanatomy, we first map the high-dimensional, voxel-level gene expression data to low-dimensional space in which data visualization can be readily achieved.
numerous multivariate analysis methods can be used for this purpose.
however, traditional methods either retain the global structures or the local structures in computing the mapping, producing results that are not satisfactory.
to preserve both the local and the global structures in the spatial gene expression space, we employ a recent method known as the t-distributed stochastic neighbor embedding   <cit>  for mapping the high-dimensional data.
this method is able to capture the local similarities in the high-dimensional space, while retaining the global structures as much as possible.
we map the high-dimensional gene expression data to 2-d space using t-sne and visualize the reduced data at multiple levels of the allen developing mouse brain reference atlas ontology, which was created based on the “prosomeric model”  <cit> .
this models proposes that the neural tube is divided into grid-like pattern of longitudinal and transverse regions.
our results show that the brain anatomy is largely preserved in the low-dimensional gene expression space at multiple levels.
to provide a quantitative comparison of the relationship between gene expression patterns and neuroanatomy, we cluster the brain voxels into groups based on gene expression data in the original high-dimensional space and in the dimensionality-reduced space.
our results show that the clustering results in the low-dimensional space are more consistent with developmental anatomy than those in the original high-dimensional space.
methods
allen developing mouse brain atlas
the allen developing mouse brain atlas  contains spatiotemporal in situ hybridization  gene expression data across multiple stages of mouse brain development  <cit> .
the primary data consist of 3-d, cellular resolution ish expression patterns of approximately  <dig> genes in sagittal plane across four embryonic  and three early postnatal ages .
the ish image series are processed by an informatics pipeline at the allen institute for brain science  <cit> .
to establish a common coordinate framework for analyzing the ish data, the ish image series are aligned to the reference atlas in 3-d space.
after the ish image series are mapped to the reference space, a gridding module is applied to divide the 3-d reference space into regular grid.
the resolution of the data grids varies with age and are shown in table  <dig> along with the sizes of each dimension.
for each grid voxel, an expression energy value is extracted.
all downstream analysis functions provided by the allen brain atlas, such as the anatomic search, gene search, and temporal search, are based on the expression energy.
our analysis in this work is also based on the grid-level expression energy.
the 3-d reference space is in pir orientation where x axis corresponds to anterior-to-posterior, y axis corresponds to superior-to-inferior, and z axis corresponds to left-to-right.
the reference atlas ontology was created based on the prosomeric model, which proposes that the developing brain is divided along the transversal and longitudinal boundaries, giving rise to a grid-like pattern .
the ontology was designed to capture the progressive development and regionalization of the nervous system.
an ontological term at each level has multiple child terms at the next level, reflecting the subdivision of the corresponding structure into multiple substructures.
by this construction, the reference atlas ontology forms a 13-level hierarchy in which the root corresponds to the undivided neural plate.
the ontology was colorized so that anatomically and developmentally related structures are coded with similar colors.
the ontology from level  <dig> to level  <dig> is shown in figures  <dig> and  <dig>
dimensionality reduction and visualization
dimensionality reduction is the procedure of mapping high-dimensional data points to low-dimensional space by optimizing certain criterion.
such techniques facilitate visual exploration of the high-dimensional data when they are mapped to 2-d or 3-d space.
traditional techniques for dimensionality reduction include linear method such as principal component analysis , multidimensional scaling , and nonlinear approaches such as local linear embedding   <cit> .
these techniques either capture the global structure of the original data or try to retain the local structure within the neighborhood of each data point.
in order to capture both the local structure and the global structure such as the presence of clusters, a class of methods, known as the stochastic neighbor embedding , have been developed  <cit> .
to simplify the optimization and overcome the so-called “crowding problem”, sne is extended to t-distributed sne  in  <cit> .
given n high-dimensional data points {xi}i=1n where xi∈ℝd, t-sne computes n low-dimensional data points {yi}i=1n, known as map points, by trying to preserve the pairwise similarities in the high-dimensional space.
to this end, t-sne computes an n × n similarity matrix in both the original data space and in the low-dimensional space.
the similarity matrix in the high-dimensional space is obtained based on symmetrized gaussian conditional distributions, while that in the low-dimensional space is computed from student t-distributions.
the map points are learned by minimizing the kullback-leibler  divergence between the probability distributions in the original data space and the embedding space.
to map our ish gene expression data, xi represents the high-dimensional gene expression vector of the ith voxel, and yi represents its representation in the low-dimensional space.
in the high-dimensional space, we define the similarity of data point xj to data point xi as the conditional probability pj|i, which captures the probability that xi chooses to have xj as its neighbor when neighbors are selected in proportion to their probability density under a gaussian distribution centered at xi.
formally, the conditional probabilities pj|i and pi|j are defined as
pj|i=exp−xi−xj2/2σi2∑k≠iexp−xk−xi2/2σi <dig>
pi|j=exp−xj−xi2/2σj2∑k≠jexp−xk−xj2/2σj <dig>
where σi <dig> denotes the variance of the gaussian distribution centered at xi.
the variance for each data point is tuned separately based on a fixed perplexity specified by the user.
the perplexity can be interpreted as a smooth measure of the effective number of neighbors, and it has been shown that the performance of t-sne is robust to changes in the perplexity  <cit> .
the recommended range for perplexity value is  <dig> to  <dig>  and the default value of  <dig> is used in the experiments.
note that pj|i≠pi|j, and pi|i =  <dig> for all i, since only pairwise similarities are of interests.
then the pairwise similarities in the high-dimensional space are symmetrized as
pij=pj|i+pi|j <dig>
the original sne method employs gaussian distributions to derive the pairwise similarities in the low-dimensional space.
this, however, leads to the crowding problem  <cit> .
to overcome this limitation, the distances in the low-dimensional space are converted into probabilities using a heavy-tailed student t-distribution in t-sne as
qij=1+yi−yj2−1∑k≠ℓ1+yk−yℓ2− <dig>
to learn the map points {yi}i=1n, t-sne minimizes the kl divergence between the probability distribution p and q in the high-dimensional and low-dimensional spaces as
minqkl=∑i∑jpijlogpijqij.
because kl divergence is not symmetric, different types of mismatches contribute differently to the overall cost.
specifically, a large cost will be induced if distant map points are used to represent nearby original data points, while a small cost is incurred if distant original data points are mapped to nearby map points.
this indicates that t-sne is able to preserve the local structure of the high-dimensional data points.
it has been shown that the objective function of t-sne is particularly straightforward to optimize in comparison to the original sne objective.
the original algorithm in  <cit>  for computing the low-dimensional map points has a time and space complexity of o, where n is the number of data points.
in  <cit> , a more efficient algorithm, known as the barnes-hut-sne, is developed, and it has o time and o space complexity.
this enables the application of t-sne to the large-scale allen developing mouse brain atlas data.
the implementations of t-sne can be found at http://homepage.tudelft.nl/19j49/t-sne.html.
clustering
to study the relationship between spatial gene expression patterns and classical neuroanatomy in the adult mouse brain, bohland et al.
<cit>  use the allen mouse brain atlas data  <cit>  and apply principal component analysis  to reduce the data dimensionality before the k-means algorithm is used to cluster the brain voxels into groups.
to visualize the spatial gene expression patterns, they also map the high-dimensional gene expression data to 3-d space using pca and visualize the data using scatter plots.
following  <cit> , we apply the k-means clustering algorithms to group brain voxels into clusters based on the gene expression data in both the original high-dimensional space and the dimensionality-reduced space.
since the results of the k-means algorithm depend on the initial cluster centers that are randomly selected, we repeat this algorithm  <dig> times and use the results with the smallest within-cluster sum of squares error.
the number of clusters in k-means is set to be equal to the number of brain structures at each particular ontology level.
we reduce the high-dimensional gene expression data to 2-d and 10-d spaces using t-sne and pca and then apply the k-means algorithm to cluster the voxels based on these low-dimensional representations.
we then quantitatively compare the consistency between voxel clusters and the neuroanatomy at multiple levels in the reference atlas developmental ontology.
we employ four performance measures, including the normalized mutual information , adjusted rand index , purity, and s-index, to evaluate the consistency between clustering results and developmental neuroanatomy.
the first three measures have been commonly used in the clustering community as external criteria for evaluating clustering results  <cit> , and the ari and s-index have been used for comparing different brain parcellation schemes  <cit> .
we treat the voxel annotations as their class labels and compare them with the clustering results.
in computing purity, each cluster is assigned to the most frequent class in the cluster, and then the final measure is the proportion of correctly assigned samples.
one disadvantage of purity is that it cannot trade off the quality of the clustering against the number of clusters  <cit> .
this limitation can be overcome by the nmi, which measures the amount of  information by which our knowledge about the classes increases when we are given the clustering results.
the ari computes the normalized fraction of all possible pairs of voxels that  have the same class label and are assigned to the same cluster or  have different class labels and are assigned to different clusters.
the s-index was specifically designed to compare different brain parcellations, and it “penalizes” class-to-cluster relationships that are overlapping, but that are not pure subset relationships  <cit> .
different measures capture different aspects of class-to-cluster consistency, and thus the trend of performance by different measures might not always be the same.
RESULTS
we retrieve the ish expression energy grid data, the reference atlas ontology and annotation data for all seven developmental ages from the allen brain atlas api  <cit> .
we remove voxels in the spinal cord for all developmental ages, as our primary goal is to study the brain gene expression and anatomy.
we also remove voxels that are annotated to a level less than level  <dig>  since we are interested in studying the relationship between the spatial gene expression patterns and the transversal and longitudinal grid-like domains that correspond to level  <dig> annotations in the reference atlas ontology.
in the current release of the data , most of the annotations have been worked down to levels between  <dig> and 8; hence only a small number of voxels were removed in this step.
after this step, all voxels are annotated with structures at levels between  <dig> and  <dig>  to study the developing mouse brain anatomy at multiple levels of granularity, we propagate the annotation of each voxel up to level  <dig>  level  <dig>  and level  <dig>  resulting in three annotated structures for each voxel that correspond to ancestor-child relations in the reference atlas ontology.
the statistics of the data sets that are used in this work are shown in table  <dig>  the input data to the t-sne method for each developmental stage is a data matrix of size n×d, where n is the number of voxels, and d is the number of genes.
the data sets span seven developmental ages, and the number of genes and the number of voxels are reported.
data visualization at multiple ontology levels
to visually explore the relationship between spatial gene expression patterns and brain neuroanatomy, we project the high-dimensional, voxel-level gene expression vectors onto 2-d space using t-sne and pca.
in pca, the data matrices are centered by subtracting the mean.
to investigate this relationship at multiple levels of the ontology, we display each projected data point using its level  <dig>  level  <dig>  and level  <dig> annotations, where the structure abbreviation is used as the marker that is color-coded according to its reference atlas ontology color.
the full names of structures can be found in figure  <dig>  we show the results generated by t-sne and pca using level  <dig>  level  <dig>  and level  <dig> annotations in figures  <dig> and  <dig> for ages e <dig>  and p <dig>  respectively.
the complete set of visualization results for all other ages are included in the additional file  <dig>
we observe that t-sne is better at visualizing the high-dimensional gene expression data than pca.
specifically, we can observe that, at all developmental ages, the three major brain structures at level  <dig>  are very well separated.
the results by t-sne preserve the brain anatomy more faithfully than those by pca at this level.
the second rows of figures  <dig> and  <dig> show the results by t-sne and pca displayed using the level  <dig> annotations, which identify the major transversal segments.
we can observe that both the global and local brain structures at this level are largely preserved in the dimensionality-reduced gene expression data space.
the third rows of figures  <dig> and  <dig> show the scatter plots of reduced data displayed using the level  <dig> annotations, which identify the four longitudinal zones in addition to the transversal segments.
we can observe that within each of the transversal segments, voxels belong to the same longitudinal zones are usually placed close to each other.
however, voxels in the same longitudinal zone but belong to different transversal segments are not necessarily placed at nearby locations.
we can observe from figures  <dig> and  <dig> that t-sne is able to map high-dimensional data to 2-d space in which the neuroanatomy can be largely recovered.
for example, in figures  <dig> and  <dig> the overall organization of the three brain structures at level  <dig> are largely preserved, where the midbrain voxels are placed between the forebrain and hindbrain voxels.
these results indicate that t-sne is able to preserve both the local and the global structures of the data simultaneously.
in addition, the shapes of the structures are also preserved to some extent.
for example, it is known that the midbrain is a wedge-shaped structure due to the sharp flexion of the neuraxis in this region  <cit> .
we can see from figures  <dig> and  <dig> that this is largely preserved in most plots.
this is especially clear from plot for the developmental age e <dig> .
this is presumably due to the much larger number of voxels in late ages , which prevent some global structures from being fully incorporated.
at level  <dig> shown in figures  <dig> and  <dig>  the transversal segment structures are also largely preserved.
in particular, p <dig> voxels are almost always close to the midbrain voxels, while p <dig> voxels are usually on the secondary prosencephalon side.
m <dig> voxels are mostly placed closely with p <dig> voxels, while m <dig> voxels are nearby with hindbrain voxels.
in the hindbrain, prepontine hindbrain voxels  are mostly close to midbrain voxels; medullary hindbrain voxels  are placed on the far side; pontine hindbrain  and pontomedullary hindbrain  voxels are somewhere in between.
we also observe that the global brain structures are less well preserved at late developmental ages.
this might be due to the increasingly larger number of brain voxels at late ages, which makes it increasingly difficult to preserve both the global and the local structures.
in this case, t-sne tends to focus more on retaining the local structure due to the asymmetric nature of the kl divergence.
clustering and comparison with neuroanatomy
since our visual explorations have shown that the brain anatomy is largely preserved in the dimensionality-reduced space, we expect that grouping of voxels into clusters based on the low-dimensional representations might lead to voxel clusters that are more consistent with neuroanatomy than those obtained from the original high-dimensional representations.
we use t-sne and pca to reduce the data to 2-d and 10-d spaces and then apply the k-means clustering algorithm to group the low-dimensional representations.
we employ four performance measures to evaluate the consistency between clustering results and neuroanatomy.
the results at developmental ontology levels  <dig>   <dig>  and  <dig>  respectively, are reported in tables  <dig>   <dig>  and  <dig>
“k-means” denotes applications of the k-means to the original high-dimensional gene expression data; “pca2” and “pca10” denote applications of k-means to pca reduced data, where the first  <dig> and  <dig> principal components are used, respectively; “t-sne2” and “t-sne10” denote applications of k-means to data reduced to  <dig> and  <dig> dimensions, respectively, by t-sne.
four different measures, including normalized mutual information , s-index, adjusted rand index , and purity, are used to compare the clustering results to the reference atlas annotations.
see the footnote of table  <dig> for detailed explanations.
see the footnote of table  <dig> for detailed explanations.
we can observe from table  <dig> that the results from low-dimensional representations computed by t-sne are much more consistent with neuroanatomy than those from the original representations at level  <dig>  on average, the performance measured by nmi and s-index has been more than doubled, and that by adjusted rand index has been increased from  <dig>  to  <dig> .
on the other hand, the results from pca-reduced data are similar to those by the original data.
this is consistent with the visualization results that pca-reduced data fail to separate voxels from different brain structures clearly at this level.
we also observe that the results of pca are similar to those by the original data for measures nmi, ari, and purity.
for s-index, these two sets of results are not similar.
this might indicate that s-index measures class-to-cluster consistency in a different way than other measures.
as has been mentioned in section “clustering”, s-index penalizes class-to-cluster relationships that are overlapping, but that are not pure subset relationships  <cit> .
the other three measures are not specifically designed to capture such relationship.
at levels  <dig> and  <dig>  we can observe from tables  <dig> and  <dig> that, on average, the clustering results based on the t-sne reduced data are more consistent with the neuroanatomy than those by the original data.
in addition, the t-sne results are more consistent with the neuroanatomy than those by pca for measures nmi, ari, and purity.
the pca-reduced data give better performance than the original and the t-sne reduced data for measure s-index.
this again indicates that s-index measures consistency in a different way compared with the other three measures.
we can conclude from the above results that, although t-sne gives better visualization results than pca at all levels, the clustering results based on pca-reduced data could yield higher consistency with the neuroanatomy than those based on t-sne for certain measure.
these results are consistent with the results reported in  <cit> .
dimensionality reduction by t-sne and pca
we observe that t-sne gives the best results in terms of preserving both the local and the global structures in the high-dimensional gene expression space in comparison with pca.
we also observe that when the data sets are very large, such as those in late developmental ages of the allen developing mouse brain atlas, preserving both the local and the global structures might be very hard or even impossible.
in these cases, t-sne tries to preserve local structures at the price of losing some global structures.
this tradeoff is achieved by giving different costs to different types of errors in computing the mapping.
in particular, because kl divergence is not symmetric, different types of mismatches contribute differently to the overall cost.
a large cost will be induced if distant map points are used to represent nearby original data points.
this large cost will ensure that the local structures are faithfully preserved.
in contrast, a relatively small cost is incurred if distant original data points are mapped to nearby map points.
hence, a small cost will be incurred if the global structures are not preserved accurately.
this asymmetric property makes t-sne especially useful in reducing and visualizing large-scale brain data sets in comparison to other traditional techniques, which preserve either the global or the local structures.
longitudinal zones versus transversal segments
in developmental neuroanatomy, two primary models have been proposed to explain the neural plate and tube regionalization based on gene expression and morphological information  <cit> .
these are the topographic “columnar” model  <cit> , and the topological “segmental” model known as the “prosomeric model”  <cit> .
recent experimental data have shown that the prosomeric model is more consistent with morphological and molecular evidences.
this leads to the adoption of this model in the allen developing mouse brain reference atlas.
the columnar model mainly focuses on dividing the neural plate and tube along the longitudinal dimension, while the segmental model favors division into transversal domains.
in the prosomeric model , the developing nervous system is divided into a grid-like pattern of longitudinal and transversal histogenetic domains.
along the longitudinal axis, four zones, known as the floor plate, basal plate, alar plate, and roof plate, are specified by dv patterning mechanisms.
along the transversal axis, the ap patterning signals subdivide the brain wall into a constant set of segments known as neuromeres.
to provide in-depth visual exploration of the genetic neuroanatomy along the longitudinal and transversal dimensions, we display in figure  <dig> the e <dig>  and p <dig> data sets according to the longitudinal zone that each voxel belongs to.
these results can be compared with the level  <dig> visualizations in figures  <dig> and  <dig>  which displays the reduced data according to the transversal segment that each voxel belongs to.
we can observe from the t-sne results that voxels from the same longitudinal zones do not form clear clusters in comparison to the clustering patterns along the transversal dimension.
in general, voxels belongs to the alar plate and basal plate form clear clusters, while those in the roof plate and floor plate tend to be widely distributed.
however, we can observe that voxels in the roof and alar plates are usually close to each other, and those in the basal and floor plates tend to form clusters.
this shows that our computational results are more consistent with the segmental model, which is also supported by recent experimental evidences.
manifold structures in developmental gene expression
we have observed that clustering of the low-dimension representations generated by t-sne leads to more consistent results with neuroanatomy than those by the original and the pca-reduced representations.
this might indicate that the original gene expression data lie on a low-dimensional manifold in the high-dimensional space.
in addition, a general trend that we have observed in comparing the clustering results with neuroanatomy is that clustering using the low-dimensional representations gives very significant performance improvement at level  <dig> in comparison to those by the original and the pca-reduced representations.
this improvement decreases as we move to level  <dig> and level  <dig>  such trend is consistent with our hypothesis that the original gene expression data lie on a manifold in the high-dimensional space, because the level  <dig> structures are simpler and thus are easier to capture by low-dimensional representations than those at level  <dig> and level  <dig>  hence, embedding of the simple manifold into low-dimensional space facilitates the faithful characterization of the underlying structures.
on the other hand, reducing relatively complex manifold structures to low-dimensional space might not lead to better representations.
CONCLUSIONS
we employ global computational analysis to study the relationship between gene expression patterns and neuroanatomy in the developing mouse brain.
to enable visual explorations, we map the high-dimensional ish gene expression data to low-dimensional space by preserving both the local and the global structures.
this unsupervised, data-driven mapping of spatial gene expression data leads to low-dimensional representations that can be easily visualized.
our results show that the developmental neuroanatomy is largely preserved in the low-dimensional gene expression data space.
to provide quantitative results, we cluster both the original high-dimensional data and the low-dimensional mapped data and compare the results with the developmental neuroanatomy.
our results show that the clusters in the low-dimensional space are more consistent with developmental neuroanatomy than those in the high-dimensional space.
in this work, the data set at each developmental age is analyzed separately.
since development is a continuous process, it would be interesting to map and cluster the data by incorporating temporal smoothness constraints  <cit> .
we will explore time-varying dimensionality reduction and clustering algorithms in the future.
our results have shown that, although majority of the voxels are mapped to locations that are consistent with their anatomical annotations, there do exist some exceptions.
we will investigate these cases in the future.
competing interests
the author declares no competing interests.
supplementary material
additional file 1
the additional file contains the complete set of visualization results for developmental ages not shown in the main texts.
click here for file
acknowledgements
we thank the allen institute for brain science for making the allen developing mouse brain atlas data available.
we thank chinh dang, david feng, terri gilbert, michael hawrylycz, luis puelles, and carol thompson for assistance in interpreting the data and results.
this work was supported by research grants from the national science foundation  and old dominion university office of research.

 the brain consists of an enormous number of cells organized into structures. the structured organization of cells is the key to the functional efficiency of the brain. hence, a natural first step toward understanding the brain function would be to address basic research questions at the structure level. how cells are organized into structures? what are the functions of structures? how the structures are connected to each other? however, a fundamental difficulty of understanding brain functions at the structure level lies in that there is no universally agreed division of cells into structures. from a developmental perspective, the delicate organization of brain into structures is the consequence of stringent spatiotemporal patterning controlled by the molecular signals during development. in this process, cells at different spatial locations read different morphogenetic positional signals produced by the graded distribution of signaling molecules. these signals control the expression of a relatively small set of transcription factors, which in turn regulate the expression of a larger number of genes. this sequential cascade of expression control ultimately leads to cell differentiation and the emergence of connections and functional properties. the discovery that certain marker genes are expressed in regionally restricted patterns in the developing brain has either led to the introduction of new structural boundaries or made it possible to re-define existing boundaries at a higher resolution. currently, studies on the molecular-structural associations are beginning to reveal how the spatiotemporal gene expression patterns are related to cellular differentiation and structural development. in this article, we study the relationship between brain anatomy and spatiotemporal gene expression patterns in the developing mouse brain. this global study of developing neuroanatomy is made possible by the high-resolution, three-dimensional gene expression patterns provided by the allen brain atlas. as part of the aba, the allen developing mouse brain atlas provides spatiotemporal in situ hybridization gene expression pattern images across four embryonic and three postnatal developmental ages, yielding effectively a four-dimensional brain atlas. to establish a common coordinate framework for analyzing the ish data, the ish image series are aligned to the allen developing mouse brain reference atlas. this enables the global, computational study of the spatiotemporal gene expression patterns of many genes and comparison of the results with developmental anatomy. to enable visual explorations of the gene expression patterns and correlate the results with classically defined neuroanatomy, we first map the high-dimensional, voxel-level gene expression data to low-dimensional space in which data visualization can be readily achieved. numerous multivariate analysis methods can be used for this purpose. however, traditional methods either retain the global structures or the local structures in computing the mapping, producing results that are not satisfactory. to preserve both the local and the global structures in the spatial gene expression space, we employ a recent method known as the t-distributed stochastic neighbor embedding for mapping the high-dimensional data. this method is able to capture the local similarities in the high-dimensional space, while retaining the global structures as much as possible. we map the high-dimensional gene expression data to 2-d space using t-sne and visualize the reduced data at multiple levels of the allen developing mouse brain reference atlas ontology, which was created based on the “prosomeric model”. this models proposes that the neural tube is divided into grid-like pattern of longitudinal and transverse regions. our results show that the brain anatomy is largely preserved in the low-dimensional gene expression space at multiple levels. to provide a quantitative comparison of the relationship between gene expression patterns and neuroanatomy, we cluster the brain voxels into groups based on gene expression data in the original high-dimensional space and in the dimensionality-reduced space. our results show that the clustering results in the low-dimensional space are more consistent with developmental anatomy than those in the original high-dimensional space. methods allen developing mouse brain atlas the allen developing mouse brain atlas contains spatiotemporal in situ hybridization gene expression data across multiple stages of mouse brain development. the primary data consist of 3-d, cellular resolution ish expression patterns of approximately genes in sagittal plane across four embryonic and three early postnatal ages. the ish image series are processed by an informatics pipeline at the allen institute for brain science. to establish a common coordinate framework for analyzing the ish data, the ish image series are aligned to the reference atlas in 3-d space. after the ish image series are mapped to the reference space, a gridding module is applied to divide the 3-d reference space into regular grid. the resolution of the data grids varies with age and are shown in table along with the sizes of each dimension. for each grid voxel, an expression energy value is extracted. all downstream analysis functions provided by the allen brain atlas, such as the anatomic search, gene search, and temporal search, are based on the expression energy. our analysis in this work is also based on the grid-level expression energy. the 3-d reference space is in pir orientation where x axis corresponds to anterior-to-posterior, y axis corresponds to superior-to-inferior, and z axis corresponds to left-to-right. the reference atlas ontology was created based on the prosomeric model, which proposes that the developing brain is divided along the transversal and longitudinal boundaries, giving rise to a grid-like pattern. the ontology was designed to capture the progressive development and regionalization of the nervous system. an ontological term at each level has multiple child terms at the next level, reflecting the subdivision of the corresponding structure into multiple substructures. by this construction, the reference atlas ontology forms a 13-level hierarchy in which the root corresponds to the undivided neural plate. the ontology was colorized so that anatomically and developmentally related structures are coded with similar colors. the ontology from level to level is shown in figures and dimensionality reduction and visualization dimensionality reduction is the procedure of mapping high-dimensional data points to low-dimensional space by optimizing certain criterion. such techniques facilitate visual exploration of the high-dimensional data when they are mapped to 2-d or 3-d space. traditional techniques for dimensionality reduction include linear method such as principal component analysis, multidimensional scaling, and nonlinear approaches such as local linear embedding. these techniques either capture the global structure of the original data or try to retain the local structure within the neighborhood of each data point. in order to capture both the local structure and the global structure such as the presence of clusters, a class of methods, known as the stochastic neighbor embedding, have been developed. to simplify the optimization and overcome the so-called “crowding problem”, sne is extended to t-distributed sne in. given n high-dimensional data points {xi}i=1n where xi∈ℝd, t-sne computes n low-dimensional data points {yi}i=1n, known as map points, by trying to preserve the pairwise similarities in the high-dimensional space. to this end, t-sne computes an n × n similarity matrix in both the original data space and in the low-dimensional space. the similarity matrix in the high-dimensional space is obtained based on symmetrized gaussian conditional distributions, while that in the low-dimensional space is computed from student t-distributions. the map points are learned by minimizing the kullback-leibler divergence between the probability distributions in the original data space and the embedding space. to map our ish gene expression data, xi represents the high-dimensional gene expression vector of the ith voxel, and yi represents its representation in the low-dimensional space. in the high-dimensional space, we define the similarity of data point xj to data point xi as the conditional probability pj|i, which captures the probability that xi chooses to have xj as its neighbor when neighbors are selected in proportion to their probability density under a gaussian distribution centered at xi. formally, the conditional probabilities pj|i and pi|j are defined as pj|i=exp−xi−xj2/2σi2∑k≠iexp−xk−xi2/2σi pi|j=exp−xj−xi2/2σj2∑k≠jexp−xk−xj2/2σj where σi denotes the variance of the gaussian distribution centered at xi. the variance for each data point is tuned separately based on a fixed perplexity specified by the user. the perplexity can be interpreted as a smooth measure of the effective number of neighbors, and it has been shown that the performance of t-sne is robust to changes in the perplexity. the recommended range for perplexity value is to and the default value of is used in the experiments. note that pj|i≠pi|j, and pi|i =  for all i, since only pairwise similarities are of interests. then the pairwise similarities in the high-dimensional space are symmetrized as pij=pj|i+pi|j the original sne method employs gaussian distributions to derive the pairwise similarities in the low-dimensional space. this, however, leads to the crowding problem. to overcome this limitation, the distances in the low-dimensional space are converted into probabilities using a heavy-tailed student t-distribution in t-sne as qij=1+yi−yj2−1∑k≠ℓ1+yk−yℓ2− to learn the map points {yi}i=1n, t-sne minimizes the kl divergence between the probability distribution p and q in the high-dimensional and low-dimensional spaces as minqkl=∑i∑jpijlogpijqij. because kl divergence is not symmetric, different types of mismatches contribute differently to the overall cost. specifically, a large cost will be induced if distant map points are used to represent nearby original data points, while a small cost is incurred if distant original data points are mapped to nearby map points. this indicates that t-sne is able to preserve the local structure of the high-dimensional data points. it has been shown that the objective function of t-sne is particularly straightforward to optimize in comparison to the original sne objective. the original algorithm in for computing the low-dimensional map points has a time and space complexity of o, where n is the number of data points. in, a more efficient algorithm, known as the barnes-hut-sne, is developed, and it has o time and o space complexity. this enables the application of t-sne to the large-scale allen developing mouse brain atlas data. the implementations of t-sne can be found at http://homepage.tudelft.nl/19j49/t-sne.html. clustering to study the relationship between spatial gene expression patterns and classical neuroanatomy in the adult mouse brain, bohland et al. use the allen mouse brain atlas data and apply principal component analysis to reduce the data dimensionality before the k-means algorithm is used to cluster the brain voxels into groups. to visualize the spatial gene expression patterns, they also map the high-dimensional gene expression data to 3-d space using pca and visualize the data using scatter plots. following, we apply the k-means clustering algorithms to group brain voxels into clusters based on the gene expression data in both the original high-dimensional space and the dimensionality-reduced space. since the results of the k-means algorithm depend on the initial cluster centers that are randomly selected, we repeat this algorithm times and use the results with the smallest within-cluster sum of squares error. the number of clusters in k-means is set to be equal to the number of brain structures at each particular ontology level. we reduce the high-dimensional gene expression data to 2-d and 10-d spaces using t-sne and pca and then apply the k-means algorithm to cluster the voxels based on these low-dimensional representations. we then quantitatively compare the consistency between voxel clusters and the neuroanatomy at multiple levels in the reference atlas developmental ontology. we employ four performance measures, including the normalized mutual information, adjusted rand index, purity, and s-index, to evaluate the consistency between clustering results and developmental neuroanatomy. the first three measures have been commonly used in the clustering community as external criteria for evaluating clustering results, and the ari and s-index have been used for comparing different brain parcellation schemes. we treat the voxel annotations as their class labels and compare them with the clustering results. in computing purity, each cluster is assigned to the most frequent class in the cluster, and then the final measure is the proportion of correctly assigned samples. one disadvantage of purity is that it cannot trade off the quality of the clustering against the number of clusters. this limitation can be overcome by the nmi, which measures the amount of information by which our knowledge about the classes increases when we are given the clustering results. the ari computes the normalized fraction of all possible pairs of voxels that have the same class label and are assigned to the same cluster or have different class labels and are assigned to different clusters. the s-index was specifically designed to compare different brain parcellations, and it “penalizes” class-to-cluster relationships that are overlapping, but that are not pure subset relationships. different measures capture different aspects of class-to-cluster consistency, and thus the trend of performance by different measures might not always be the same.  we retrieve the ish expression energy grid data, the reference atlas ontology and annotation data for all seven developmental ages from the allen brain atlas api. we remove voxels in the spinal cord for all developmental ages, as our primary goal is to study the brain gene expression and anatomy. we also remove voxels that are annotated to a level less than level since we are interested in studying the relationship between the spatial gene expression patterns and the transversal and longitudinal grid-like domains that correspond to level annotations in the reference atlas ontology. in the current release of the data, most of the annotations have been worked down to levels between and 8; hence only a small number of voxels were removed in this step. after this step, all voxels are annotated with structures at levels between and to study the developing mouse brain anatomy at multiple levels of granularity, we propagate the annotation of each voxel up to level level and level resulting in three annotated structures for each voxel that correspond to ancestor-child relations in the reference atlas ontology. the statistics of the data sets that are used in this work are shown in table the input data to the t-sne method for each developmental stage is a data matrix of size n×d, where n is the number of voxels, and d is the number of genes. the data sets span seven developmental ages, and the number of genes and the number of voxels are reported. data visualization at multiple ontology levels to visually explore the relationship between spatial gene expression patterns and brain neuroanatomy, we project the high-dimensional, voxel-level gene expression vectors onto 2-d space using t-sne and pca. in pca, the data matrices are centered by subtracting the mean. to investigate this relationship at multiple levels of the ontology, we display each projected data point using its level level and level annotations, where the structure abbreviation is used as the marker that is color-coded according to its reference atlas ontology color. the full names of structures can be found in figure we show the results generated by t-sne and pca using level level and level annotations in figures and for ages e and p respectively. the complete set of visualization results for all other ages are included in the additional file we observe that t-sne is better at visualizing the high-dimensional gene expression data than pca. specifically, we can observe that, at all developmental ages, the three major brain structures at level are very well separated. the results by t-sne preserve the brain anatomy more faithfully than those by pca at this level. the second rows of figures and show the results by t-sne and pca displayed using the level annotations, which identify the major transversal segments. we can observe that both the global and local brain structures at this level are largely preserved in the dimensionality-reduced gene expression data space. the third rows of figures and show the scatter plots of reduced data displayed using the level annotations, which identify the four longitudinal zones in addition to the transversal segments. we can observe that within each of the transversal segments, voxels belong to the same longitudinal zones are usually placed close to each other. however, voxels in the same longitudinal zone but belong to different transversal segments are not necessarily placed at nearby locations. we can observe from figures and that t-sne is able to map high-dimensional data to 2-d space in which the neuroanatomy can be largely recovered. for example, in figures and the overall organization of the three brain structures at level are largely preserved, where the midbrain voxels are placed between the forebrain and hindbrain voxels. these results indicate that t-sne is able to preserve both the local and the global structures of the data simultaneously. in addition, the shapes of the structures are also preserved to some extent. for example, it is known that the midbrain is a wedge-shaped structure due to the sharp flexion of the neuraxis in this region. we can see from figures and that this is largely preserved in most plots. this is especially clear from plot for the developmental age e. this is presumably due to the much larger number of voxels in late ages, which prevent some global structures from being fully incorporated. at level shown in figures and the transversal segment structures are also largely preserved. in particular, p voxels are almost always close to the midbrain voxels, while p voxels are usually on the secondary prosencephalon side. m voxels are mostly placed closely with p voxels, while m voxels are nearby with hindbrain voxels. in the hindbrain, prepontine hindbrain voxels are mostly close to midbrain voxels; medullary hindbrain voxels are placed on the far side; pontine hindbrain and pontomedullary hindbrain voxels are somewhere in between. we also observe that the global brain structures are less well preserved at late developmental ages. this might be due to the increasingly larger number of brain voxels at late ages, which makes it increasingly difficult to preserve both the global and the local structures. in this case, t-sne tends to focus more on retaining the local structure due to the asymmetric nature of the kl divergence. clustering and comparison with neuroanatomy since our visual explorations have shown that the brain anatomy is largely preserved in the dimensionality-reduced space, we expect that grouping of voxels into clusters based on the low-dimensional representations might lead to voxel clusters that are more consistent with neuroanatomy than those obtained from the original high-dimensional representations. we use t-sne and pca to reduce the data to 2-d and 10-d spaces and then apply the k-means clustering algorithm to group the low-dimensional representations. we employ four performance measures to evaluate the consistency between clustering results and neuroanatomy. the results at developmental ontology levels and respectively, are reported in tables and “k-means” denotes applications of the k-means to the original high-dimensional gene expression data; “pca2” and “pca10” denote applications of k-means to pca reduced data, where the first and principal components are used, respectively; “t-sne2” and “t-sne10” denote applications of k-means to data reduced to and dimensions, respectively, by t-sne. four different measures, including normalized mutual information, s-index, adjusted rand index, and purity, are used to compare the clustering results to the reference atlas annotations. see the footnote of table for detailed explanations. see the footnote of table for detailed explanations. we can observe from table that the results from low-dimensional representations computed by t-sne are much more consistent with neuroanatomy than those from the original representations at level on average, the performance measured by nmi and s-index has been more than doubled, and that by adjusted rand index has been increased from to. on the other hand, the results from pca-reduced data are similar to those by the original data. this is consistent with the visualization results that pca-reduced data fail to separate voxels from different brain structures clearly at this level. we also observe that the results of pca are similar to those by the original data for measures nmi, ari, and purity. for s-index, these two sets of results are not similar. this might indicate that s-index measures class-to-cluster consistency in a different way than other measures. as has been mentioned in section “clustering”, s-index penalizes class-to-cluster relationships that are overlapping, but that are not pure subset relationships. the other three measures are not specifically designed to capture such relationship. at levels and we can observe from tables and that, on average, the clustering results based on the t-sne reduced data are more consistent with the neuroanatomy than those by the original data. in addition, the t-sne results are more consistent with the neuroanatomy than those by pca for measures nmi, ari, and purity. the pca-reduced data give better performance than the original and the t-sne reduced data for measure s-index. this again indicates that s-index measures consistency in a different way compared with the other three measures. we can conclude from the above results that, although t-sne gives better visualization results than pca at all levels, the clustering results based on pca-reduced data could yield higher consistency with the neuroanatomy than those based on t-sne for certain measure. these results are consistent with the results reported in. dimensionality reduction by t-sne and pca we observe that t-sne gives the best results in terms of preserving both the local and the global structures in the high-dimensional gene expression space in comparison with pca. we also observe that when the data sets are very large, such as those in late developmental ages of the allen developing mouse brain atlas, preserving both the local and the global structures might be very hard or even impossible. in these cases, t-sne tries to preserve local structures at the price of losing some global structures. this tradeoff is achieved by giving different costs to different types of errors in computing the mapping. in particular, because kl divergence is not symmetric, different types of mismatches contribute differently to the overall cost. a large cost will be induced if distant map points are used to represent nearby original data points. this large cost will ensure that the local structures are faithfully preserved. in contrast, a relatively small cost is incurred if distant original data points are mapped to nearby map points. hence, a small cost will be incurred if the global structures are not preserved accurately. this asymmetric property makes t-sne especially useful in reducing and visualizing large-scale brain data sets in comparison to other traditional techniques, which preserve either the global or the local structures. longitudinal zones versus transversal segments in developmental neuroanatomy, two primary models have been proposed to explain the neural plate and tube regionalization based on gene expression and morphological information. these are the topographic “columnar” model, and the topological “segmental” model known as the “prosomeric model”. recent experimental data have shown that the prosomeric model is more consistent with morphological and molecular evidences. this leads to the adoption of this model in the allen developing mouse brain reference atlas. the columnar model mainly focuses on dividing the neural plate and tube along the longitudinal dimension, while the segmental model favors division into transversal domains. in the prosomeric model, the developing nervous system is divided into a grid-like pattern of longitudinal and transversal histogenetic domains. along the longitudinal axis, four zones, known as the floor plate, basal plate, alar plate, and roof plate, are specified by dv patterning mechanisms. along the transversal axis, the ap patterning signals subdivide the brain wall into a constant set of segments known as neuromeres. to provide in-depth visual exploration of the genetic neuroanatomy along the longitudinal and transversal dimensions, we display in figure the e and p data sets according to the longitudinal zone that each voxel belongs to. these results can be compared with the level visualizations in figures and which displays the reduced data according to the transversal segment that each voxel belongs to. we can observe from the t-sne results that voxels from the same longitudinal zones do not form clear clusters in comparison to the clustering patterns along the transversal dimension. in general, voxels belongs to the alar plate and basal plate form clear clusters, while those in the roof plate and floor plate tend to be widely distributed. however, we can observe that voxels in the roof and alar plates are usually close to each other, and those in the basal and floor plates tend to form clusters. this shows that our computational results are more consistent with the segmental model, which is also supported by recent experimental evidences. manifold structures in developmental gene expression we have observed that clustering of the low-dimension representations generated by t-sne leads to more consistent results with neuroanatomy than those by the original and the pca-reduced representations. this might indicate that the original gene expression data lie on a low-dimensional manifold in the high-dimensional space. in addition, a general trend that we have observed in comparing the clustering results with neuroanatomy is that clustering using the low-dimensional representations gives very significant performance improvement at level in comparison to those by the original and the pca-reduced representations. this improvement decreases as we move to level and level such trend is consistent with our hypothesis that the original gene expression data lie on a manifold in the high-dimensional space, because the level structures are simpler and thus are easier to capture by low-dimensional representations than those at level and level hence, embedding of the simple manifold into low-dimensional space facilitates the faithful characterization of the underlying structures. on the other hand, reducing relatively complex manifold structures to low-dimensional space might not lead to better representations.  we employ global computational analysis to study the relationship between gene expression patterns and neuroanatomy in the developing mouse brain. to enable visual explorations, we map the high-dimensional ish gene expression data to low-dimensional space by preserving both the local and the global structures. this unsupervised, data-driven mapping of spatial gene expression data leads to low-dimensional representations that can be easily visualized. our results show that the developmental neuroanatomy is largely preserved in the low-dimensional gene expression data space. to provide quantitative results, we cluster both the original high-dimensional data and the low-dimensional mapped data and compare the results with the developmental neuroanatomy. our results show that the clusters in the low-dimensional space are more consistent with developmental neuroanatomy than those in the high-dimensional space. in this work, the data set at each developmental age is analyzed separately. since development is a continuous process, it would be interesting to map and cluster the data by incorporating temporal smoothness constraints. we will explore time-varying dimensionality reduction and clustering algorithms in the future. our results have shown that, although majority of the voxels are mapped to locations that are consistent with their anatomical annotations, there do exist some exceptions. we will investigate these cases in the future. competing interests the author declares no competing interests. supplementary material additional file 1 the additional file contains the complete set of visualization results for developmental ages not shown in the main texts. click here for file acknowledgements we thank the allen institute for brain science for making the allen developing mouse brain atlas data available. we thank chinh dang, david feng, terri gilbert, michael hawrylycz, luis puelles, and carol thompson for assistance in interpreting the data and results. this work was supported by research grants from the national science foundation and old dominion university office of research.
to enable visual explorations of the high-dimensional data, we map the in situ hybridization gene expression data to a two-dimensional space by preserving both the global and the local structures. the structured organization of cells in the brain plays a key role in its functional efficiency. our results show that the developing brain anatomy is largely preserved in the reduced gene expression space. our results show that the clusters in the low-dimensional space are more consistent with neuroanatomy than those in the original space. to provide a quantitative analysis, we cluster the reduced data into groups and measure the consistency with neuroanatomy at multiple levels. currently, studies on the molecular-structural association are beginning to reveal how the spatiotemporal gene expression patterns are related to cellular differentiation and structural development. gene expression patterns and developing brain anatomy are closely related. in this article, we aim at a global, data-driven study of the relationship between gene expressions and neuroanatomy in the developing mouse brain. this delicate organization is the consequence of unique molecular identity of each cell gradually established by precise spatiotemporal gene expression control during development. dimensionality reduction and visual exploration facilitate the study of this relationship.
 the structured organization of cells in the brain plays a key role in its functional efficiency. this delicate organization is the consequence of unique molecular identity of each cell gradually established by precise spatiotemporal gene expression control during development. currently, studies on the molecular-structural association are beginning to reveal how the spatiotemporal gene expression patterns are related to cellular differentiation and structural development.  in this article, we aim at a global, data-driven study of the relationship between gene expressions and neuroanatomy in the developing mouse brain. to enable visual explorations of the high-dimensional data, we map the in situ hybridization gene expression data to a two-dimensional space by preserving both the global and the local structures. our results show that the developing brain anatomy is largely preserved in the reduced gene expression space. to provide a quantitative analysis, we cluster the reduced data into groups and measure the consistency with neuroanatomy at multiple levels. our results show that the clusters in the low-dimensional space are more consistent with neuroanatomy than those in the original space.  gene expression patterns and developing brain anatomy are closely related. dimensionality reduction and visual exploration facilitate the study of this relationship. 