1.1 --- a/grant.txt Sat Apr 11 21:43:12 2009 -0700 1.2 +++ b/grant.txt Sun Apr 12 02:49:55 2009 -0700 1.3 @@ -1,7 +1,5 @@ 1.4 == Specific aims == 1.5 1.6 -todo test4 1.7 - 1.8 Massive new datasets obtained with techniques such as in situ hybridization (ISH) and BAC-transgenics allow the expression levels of many genes at many locations to be compared. Our goal is to develop automated methods to relate spatial variation in gene expression to anatomy. We want to find marker genes for specific anatomical regions, and also to draw new anatomical maps based on gene expression patterns. We have three specific aims: 1.9 1.10 (1) develop an algorithm to screen spatial gene expression data for combinations of marker genes which selectively target anatomical regions 1.11 @@ -17,12 +15,12 @@ 1.12 == Background and significance == 1.13 1.14 === Aim 1 === 1.15 -==== Machine learning terminology ==== 1.16 +==== Machine learning terminology: supervised learning ==== 1.17 The task of looking for marker genes for anatomical subregions means that one is looking for a set of genes such that, if the expression level of those genes is known, then the locations of the subregions can be inferred. 1.18 1.19 If we define the subregions so that they cover the entire anatomical structure to be divided, then instead of saying that we are using gene expression to find the locations of the subregions, we may say that we are using gene expression to determine to which subregion each voxel within the structure belongs. We call this a __classification task__, because each voxel is being assigned to a class (namely, its subregion). 1.20 1.21 -Therefore, an understanding of the relationship between the combination of their expression levels and the locations of the subregions may be expressed as a function. The input to this function is a voxel, along with the gene expression levels within that voxel; the output is the subregional identity of the target voxel, that is, the subregion to which the target voxel belongs. We call this function a __classifier__. In general, the input to a classifier is called an __instance__, and the output is called a __label__. 1.22 +Therefore, an understanding of the relationship between the combination of their expression levels and the locations of the subregions may be expressed as a function. The input to this function is a voxel, along with the gene expression levels within that voxel; the output is the subregional identity of the target voxel, that is, the subregion to which the target voxel belongs. We call this function a __classifier__. In general, the input to a classifier is called an __instance__, and the output is called a __label__ (or a __class label__). 1.23 1.24 The object of aim 1 is not to produce a single classifier, but rather to develop an automated method for determining a classifier for any known anatomical structure. Therefore, we seek a procedure by which a gene expression dataset may be analyzed in concert with an anatomical atlas in order to produce a classifier. Such a procedure is a type of a machine learning procedure. The construction of the classifier is called __training__ (also __learning__), and the initial gene expression dataset used in the construction of the classifier is called __training data__. 1.25 1.26 @@ -57,7 +55,37 @@ 1.27 1.28 1.29 === Aim 2 === 1.30 -todo 1.31 +==== Machine learning terminology: clustering ==== 1.32 +If one is given a dataset consisting merely of instances, with no class labels, then analysis of the dataset is referred to as __unsupervised learning__ in the jargon of machine learning. One thing that you can do with such a dataset is to group instances together. A set of similar instances is called a __cluster__, and the activity of finding grouping the data into clusters is called clustering or cluster analysis. 1.33 + 1.34 +The task of deciding how to carve up a structure into anatomical subregions can be put into these terms. The instances are once again voxels (or pixels) along with their associated gene expression profiles. We make the assumption that voxels from the same subregion have similar gene expression profiles, at least compared to the other subregions. This means that clustering voxels is the same as finding potential subregions; we seek a partitioning of the voxels into subregions, that is, into clusters of voxels with similar gene expression. 1.35 + 1.36 +It is desirable to determine not just one set of subregions, but also how these subregions relate to each other, if at all; perhaps some of the subregions are more similar to each other than to the rest, suggesting that, although at a fine spatial scale they could be considered separate, on a coarser spatial scale they could be grouped together into one large subregion. This suggests the outcome of clustering may be a hierarchial tree of clusters, rather than a single set of clusters which partition the voxels. This is called hierarchial clustering. 1.37 + 1.38 +==== Similarity scores ==== 1.39 + 1.40 +todo 1.41 + 1.42 +==== Spatially contiguous clusters; image segmentation ==== 1.43 + 1.44 +We have shown that aim 2 is a type of clustering task. In fact, it is a special type of clustering task because we have an additional constraint on clusters; voxels grouped together into a cluster must be spatially contiguous. In Preliminary Results, we show that one can get reasonable results without enforcing this constraint, however, we plan to compare these results against other methods which guarantee contiguous clusters. 1.45 + 1.46 +Perhaps the biggest source of continguous clustering algorithms is the field of computer vision, which has produced a variety of image segmentation algorithms. Image segmentation is the task of partitioning the pixels in a digital image into clusters, usually contiguous clusters. Aim 2 is similar to an image segmentation task. There are two main differences; in our task, there are thousands of color channels (one for each gene), rather than just three. There are imaging tasks which use more than three colors, however, for example multispectral imaging and hyperspectral imaging, which are often used to process satellite imagery. A more crucial difference is that there are various cues which are appropriate for detecting sharp object boundaries in a visual scene but which are not appropriate for segmenting abstract spatial data such as gene expression. Although many image segmentation algorithms can be expected to work well for segmenting other sorts of spatially arranged data, some of these algorithms are specialized for visual images. 1.47 + 1.48 +==== Dimensionality reduction ==== 1.49 + 1.50 +Unlike aim 1, there is no externally-imposed need to select only a handful of informative genes for inclusion in the instances. However, some clustering algorithms perform better on small numbers of features. There are techniques which "summarize" a larger number of features using a smaller number of features; these techniques go by the name of feature extraction or dimensionality reduction. The small set of features that such a technique yields is called the __reduced feature set__. After the reduced feature set is created, the instances may be replaced by __reduced instances__, which have as their features the reduced feature set rather than the original feature set of all gene expression levels. Note that the features in the reduced feature set do not necessarily correspond to genes; each feature in the reduced set may be any function of the set of gene expression levels. 1.51 + 1.52 +Another use for dimensionality reduction is to visualize the relationships between subregions. For example, one might want to make a 2-D plot upon which each subregion is represented by a single point, and with the property that subregions with similar gene expression profiles should be nearby on the plot (that is, the property that distance between pairs of points in the plot should be proportional to some measure of dissimilarity in gene expression). It is likely that no arrangement of the points on a 2-D plan will exactly satisfy this property -- however, dimensionality reduction techniques allow one to find arrangements of points that approximately satisfy that property. Note that in this application, dimensionality reduction is being applied after clustering; whereas in the previous paragraph, we were talking about using dimensionality reduction before clustering. 1.53 + 1.54 +==== Clustering genes rather than voxels ==== 1.55 + 1.56 +Although the ultimate goal is to cluster the instances (voxels or pixels), one strategy to achieve this goal is to first cluster the features (genes). There are two ways that clusters of genes could be used. 1.57 + 1.58 +Gene clusters could be used as part of dimensionality reduction: rather than have one feature for each gene, we could have one reduced feature for each gene cluster. 1.59 + 1.60 +Gene clusters could also be used to directly yield a clustering on instances. This is because many genes have an expression pattern which seems to pick out a single, spatially continguous subregion. Therefore, it seems likely that an anatomically interesting subregion will have multiple genes which each individually pick it out\footnote{This would seem to contradict our finding in aim 1 that some cortical areas are combinatorially coded by multiple genes. However, it is possible that the currently accepted cortical maps divide the cortex into subregions which are unnatural from the point of view of gene expression; perhaps there is some other way to map the cortex for which each subregion can be identified by single genes.}. This suggests the following procedure: cluster together genes which pick out similar subregions, and then to use the more popular common subregions as the final clusters. In the Preliminary Data we show that a number of anatomically recognized cortical regions, as well as some "superregions" formed by lumping together a few regions, are associated with gene clusters in this fashion. 1.61 + 1.62 1.63 1.64 1.65 @@ -85,11 +113,19 @@ 1.66 === Related work === 1.67 todo 1.68 1.69 +vs. AGEA -- i wrote something on this but i'm going to rewrite it 1.70 + 1.71 1.72 == Preliminary work == 1.73 1.74 -=== Justification of principles 1 thur 3 === 1.75 -==== Principle 1: Combinatorial gene expression ==== 1.76 +=== Format conversion between SEV, MATLAB, NIFTI === 1.77 +todo 1.78 + 1.79 +=== Flatmap of cortex === 1.80 +todo 1.81 + 1.82 + 1.83 +==== Using combinations of multiple genes is necessary and sufficient to delineate some cortical areas ==== 1.84 Here we give an example of a cortical area which is not marked by any single gene, but which can be identified combinatorially. according to logistic regression, gene wwc1\footnote{"WW, C2 and coiled-coil domain containing 1"; EntrezGene ID 211652} is the best fit single gene for predicting whether or not a pixel on the cortical surface belongs to the motor area (area MO). The upper-left picture in Figure \ref{MOcombo} shows wwc1's spatial expression pattern over the cortex. The lower-right boundary of MO is represented reasonably well by this gene, however the gene overshoots the upper-left boundary. This flattened 2-D representation does not show it, but the area corresponding to the overshoot is the medial surface of the cortex. MO is only found on the lateral surface (todo). 1.85 1.86 Gnee mtif2\footnote{"mitochondrial translational initiation factor 2"; EntrezGene ID 76784} is shown in figure the upper-right of Fig. \ref{MOcombo}. Mtif2 captures MO's upper-left boundary, but not its lower-right boundary. Mtif2 does not express very much on the medial surface. By adding together the values at each pixel in these two figures, we get the lower-left of Figure \ref{MOcombo}. This combination captures area MO much better than any single gene. 1.87 @@ -102,13 +138,8 @@ 1.88 \caption{Upper left: $wwc1$. Upper right: $mtif2$. Lower left: wwc1 + mtif2 (each pixel's value on the lower left is the sum of the corresponding pixels in the upper row). Within each picture, the vertical axis roughly corresponds to anterior at the top and posterior at the bottom, and the horizontal axis roughly corresponds to medial at the left and lateral at the right. The red outline is the boundary of region MO. Pixels are colored approximately according to the density of expressing cells underneath each pixel, with red meaning a lot of expression and blue meaning little.} 1.89 \end{figure} 1.90 1.91 -==== Principle 2: Only look at combinations of small numbers of genes ==== 1.92 -In order to see how well one can do when looking at all genes at once, we ran a support vector machine to classify cortical surface pixels based on their gene expression profiles. We achieved classification accuracy of about 81%\footnote{Using the Shogun SVM package (todo:cite), with parameters type=GMNPSVM (multiclass b-SVM), kernal = gaussian with sigma = 0.1, c = 10, epsilon = 1e-1 -- these are the first parameters we tried, so presumably performance would improve with different choices of parameters. 5-fold cross-validation.}. As noted above, however, a classifier that looks at all the genes at once isn't practically useful. 1.93 - 1.94 -The requirement to find combinations of only a small number of genes limits us from straightforwardly applying many of the most simple techniques from the field of supervised machine learning. In the parlance of machine learning, our task combines feature selection with supervised learning. 1.95 - 1.96 - 1.97 -==== Principle 3: Use geometry ==== 1.98 + 1.99 +==== Geometric and pointwise scoring methods provide complementary information ==== 1.100 1.101 1.102 To show that local geometry can provide useful information that cannot be detected via pointwise analyses, consider Fig. \ref{AUDgeometry}. The top row of Fig. \ref{AUDgeometry} displays the 3 genes which most match area AUD, according to a pointwise method\footnote{For each gene, a logistic regression in which the response variable was whether or not a surface pixel was within area AUD, and the predictor variable was the value of the expression of the gene underneath that pixel. The resulting scores were used to rank the genes in terms of how well they predict area AUD.}. The bottom row displays the 3 genes which most match AUD according to a method which considers local geometry\footnote{For each gene the gradient similarity (see section \ref{gradientSim}) between (a) a map of the expression of each gene on the cortical surface and (b) the shape of area AUD, was calculated, and this was used to rank the genes.} The pointwise method in the top row identifies genes which express more strongly in AUD than outside of it; its weakness is that this includes many areas which don't have a salient border matching the areal border. The geometric method identifies genes whose salient expression border seems to partially line up with the border of AUD; its weakness is that this includes genes which don't express over the entire area. Genes which have high rankings using both pointwise and border criteria, such as $Aph1a$ in the example, may be particularly good markers. None of these genes are, individually, a perfect marker for AUD; we deliberately chose a "difficult" area in order to better contrast pointwise with geometric methods. 1.103 @@ -125,130 +156,51 @@ 1.104 \caption{The top row shows the three genes which (individually) best predict area AUD, according to logistic regression. The bottom row shows the three genes which (individually) best match area AUD, according to gradient similarity. From left to right and top to bottom, the genes are $Ssr1$, $Efcbp1$, $Aph1a$, $Ptk7$, $Aph1a$ again, and $Lepr$} 1.105 \end{figure} 1.106 1.107 - 1.108 - 1.109 - 1.110 - 1.111 -==== Principle 4: Work in 2-D whenever possible ==== 1.112 - 1.113 -In anatomy, the manifold of interest is usually either defined by a combination of two relevant anatomical axes (todo), or by the surface of the structure (as is the case with the cortex). In the former case, the manifold of interest is a plane, but in the latter case it is curved. If the manifold is curved, there are various methods for mapping the manifold into a plane. 1.114 - 1.115 -The method that we will develop will begin by mapping the data into a 2-D plane. Although the manifold that characterized cortical areas is known to be the cortical surface, it remains to be seen which method of mapping the manifold into a plane is optimal for this application. We will compare mappings which attempt to preserve size (such as the one used by Caret\ref{van_essen_integrated_2001}) with mappings which preserve angle (conformal maps). 1.116 - 1.117 -Although there is much 2-D organization in anatomy, there are also structures whose shape is fundamentally 3-dimensional. If possible, we would like the method we develop to include a statistical test that warns the user if the assumption of 2-D structure seems to be wrong. 1.118 - 1.119 - 1.120 - 1.121 - 1.122 - 1.123 - 1.124 - 1.125 ------- 1.126 - 1.127 - 1.128 - 1.129 -Massive new datasets obtained with techniques such as in situ hybridization (ISH) and BAC-transgenics allow the expression levels of many genes at many locations to be compared. This can be used to find marker genes for specific anatomical structures, as well as to draw new anatomical maps. Our goal is to develop automated methods to relate spatial variation in gene expression to anatomy. We have five specific aims: 1.130 - 1.131 -(1) develop an algorithm to screen spatial gene expression data for combinations of marker genes which selectively target individual anatomical structures 1.132 -(2) develop an algorithm to screen spatial gene expression data for combinations of marker genes which can be used to delineate most of the boundaries between a number of anatomical structures at once 1.133 -(3) develop an algorithm to suggest new ways of dividing a structure up into anatomical subregions, based on spatial patterns in gene expression 1.134 -(4) create a flat (2-D) map of the mouse cerebral cortex that contains a flattened version of the Allen Mouse Brain Atlas ISH dataset, as well as the boundaries of anatomical areas within the cortex. For each cortical layer, a layer-specific flat dataset will be created. A single combined flat dataset will be created which averages information from all of the layers. These datasets will be made available in both MATLAB and Caret formats. 1.135 -(5) validate the methods developed in (1), (2) and (3) by applying them to the cerebral cortex datasets created in (4) 1.136 - 1.137 -All algorithms that we develop will be implemented in an open-source software toolkit. The toolkit, as well as the machine-readable datasets developed in aim (4) and any other intermediate dataset we produce, will be published and freely available for others to use. 1.138 - 1.139 -In addition to developing generally useful methods, the application of these methods to cerebral cortex will produce immediate benefits that are only one step removed from clinical application, while also supporting the development of new neuroanatomical techniques. The method developed in aim (1) will be applied to each cortical area to find a set of marker genes. Currently, despite the distinct roles of different cortical areas in both normal functioning and disease processes, there are no known marker genes for many cortical areas. Finding marker genes will be immediately useful for drug discovery as well as for experimentation because once marker genes for an area are known, interventions can be designed which selectively target that area. 1.140 - 1.141 - 1.142 - 1.143 - 1.144 - 1.145 - 1.146 - 1.147 -The method developed in aim (2) will be used to find a small panel of genes that can find most of the boundaries between areas in the cortex. Today, finding cortical areal boundaries in a tissue sample is a manual process that requires a skilled human to combine multiple visual cues over a large area of the cortical surface. A panel of marker genes will allow the development of an ISH protocol that will allow experimenters to more easily identify which anatomical areas are present in small samples of cortex. 1.148 - 1.149 - 1.150 - 1.151 - 1.152 - 1.153 - 1.154 - 1.155 - 1.156 - 1.157 - 1.158 - 1.159 - 1.160 -For each cortical layer, a layer-specific flat dataset will be created. A single combined flat dataset will be created which averages information from all of the layers. These datasets will be made available in both MATLAB and Caret formats. 1.161 - 1.162 - 1.163 - 1.164 - 1.165 - 1.166 - 1.167 - 1.168 - 1.169 ----- 1.170 - 1.171 - 1.172 - 1.173 -New techniques allow the expression levels of many genes at many locations to be compared. It is thought that even neighboring anatomical structures have different gene expression profiles. We propose to develop automated methods to relate the spatial variation in gene expression to anatomy. We will develop two kinds of techniques: 1.174 - 1.175 -(a) techniques to screen for combinations of marker genes which selectively target anatomical structures 1.176 -(b) techniques to suggest new ways of dividing a structure up into anatomical subregions, based on the shapes of contours in the gene expression 1.177 - 1.178 -The first kind of technique will be helpful for finding marker genes associated with known anatomical features. The second kind of technique will be helpful in creating new anatomical maps, maps which reflect differences in gene expression the same way that existing maps reflect differences in histology. 1.179 - 1.180 -We intend to develop our techniques using the adult mouse cerebral cortex as a testbed. The Allen Brain Atlas has collected a dataset containing the expression level of about 4000 genes* over a set of over 150000 voxels, with a spatial resolution of approximately 200 microns\cite{lein_genome-wide_2007}. 1.181 - 1.182 -We expect to discover sets of marker genes that pick out specific cortical areas. This will allow the development of drugs and other interventions that selectively target individual cortical areas. Therefore our research will lead to application in drug discovery, in the development of other targeted clinical interventions, and in the development of new experimental techniques. 1.183 - 1.184 -The best way to divide up rodent cortex into areas has not been completely determined, as can be seen by the differences in the recent maps given by Swanson on the one hand, and Paxinos and Franklin on the other. It is likely that our study, by showing which areal divisions naturally follow from gene expression data, as opposed to traditional histological data, will contribute to the creation of a better map. While we do not here propose to analyze human gene expression data, it is conceivable that the methods we propose to develop could be used to suggest modifications to the human cortical map as well. 1.185 - 1.186 - 1.187 -In the following, we will only be talking about coronal data. 1.188 - 1.189 -The Allen Brain Atlas provides "Smoothed Energy Volumes", which are 1.190 - 1.191 - 1.192 -One type of artifact in the Allen Brain Atlas data is what we call a "slice artifact". We have noticed two types of slice artifacts in the dataset. The first type, a "missing slice artifact", occurs when the ISH procedure on a slice did not come out well. In this case, the Allen Brain investigators excluded the slice at issue from the dataset. This means that no gene expression information is available for that gene for the region of space covered by that slice. This results in an expression level of zero being assigned to voxels covered by the slice. This is partially but not completely ameliorated by the smoothing that is applied to create the Smoothed Energy Volumes. The usual end result is that a region of space which is shaped and oriented like a coronal slice is marked as having less gene expression than surrounding regions. 1.193 - 1.194 -The second type of slice artifact is caused by the fact that all of the slices have a consistent orientation. Since there may be artifacts (such as how well the ISH worked) which are constant within each slice but which vary between different slices, the result is that ceteris paribus, when one compares the genetic data of a voxel to another voxel within the same coronal plane, one would expect to find more similarity than if one compared a voxel to another voxel displaced along the rostrocaudal axis. 1.195 - 1.196 - 1.197 - 1.198 - 1.199 -We are enthusiastic about the sharing of methods, data, and results, and at the conclusion of the project, we will make all of our data and computer source code publically available. Our goal is that replicating our results, or applying the methods we develop to other targets, will be quick and easy for other investigators. In order to aid in understanding and replicating our results, we intend to include a software program which, when run, will take as input the Allen Brain Atlas raw data, and produce as output all numbers and charts found in publications resulting from the project. 1.200 - 1.201 - 1.202 -To aid in the replication of our results, we will include a script which takes as input the dataset in aim (3) and provides as output all of the tables in figures in our publications . 1.203 - 1.204 - 1.205 - 1.206 - 1.207 -We also expect to weigh in on the debate about how to best partition rodent cortex 1.208 - 1.209 - 1.210 - 1.211 -be useful for drug discovery as well 1.212 - 1.213 - 1.214 - 1.215 -* Another 16000 genes are available, but they do not cover the entire cerebral cortex with high spatial resolution. 1.216 - 1.217 - 1.218 -User-definable ROIs 1.219 -Combinatorial gene expression 1.220 -Negative as well as positive signal 1.221 -Use geometry 1.222 -Search for local boundaries if necessary 1.223 -Flatmapped 1.224 - 1.225 - 1.226 - 1.227 - 1.228 - 1.229 - 1.230 -== Specific aims == 1.231 +==== Areas which can be identified by single genes ==== 1.232 +todo 1.233 + 1.234 + 1.235 +=== Aim 1 (and Aim 3) === 1.236 + 1.237 + 1.238 +==== SVM on all genes at once ==== 1.239 +In order to see how well one can do when looking at all genes at once, we ran a support vector machine to classify cortical surface pixels based on their gene expression profiles. We achieved classification accuracy of about 81%\footnote{Using the Shogun SVM package (todo:cite), with parameters type=GMNPSVM (multiclass b-SVM), kernal = gaussian with sigma = 0.1, c = 10, epsilon = 1e-1 -- these are the first parameters we tried, so presumably performance would improve with different choices of parameters. 5-fold cross-validation.}. As noted above, however, a classifier that looks at all the genes at once isn't practically useful. 1.240 + 1.241 +The requirement to find combinations of only a small number of genes limits us from straightforwardly applying many of the most simple techniques from the field of supervised machine learning. In the parlance of machine learning, our task combines feature selection with supervised learning. 1.242 + 1.243 + 1.244 +==== Decision trees ==== 1.245 +todo 1.246 + 1.247 + 1.248 +=== Aim 2 (and Aim 3) === 1.249 + 1.250 +=== Raw dimensionality reduction results === 1.251 + 1.252 + 1.253 +=== Dimensionality reduction plus K-means or spectral clustering === 1.254 + 1.255 + 1.256 +==== Many areas are captured by clusters of genes ==== 1.257 +todo 1.258 + 1.259 + 1.260 + 1.261 + 1.262 + 1.263 + 1.264 + 1.265 + 1.266 + 1.267 + 1.268 + 1.269 +todo 1.270 + 1.271 +== Research plan == 1.272 + 1.273 +todo 1.274 + 1.275 +amongst other thigns: 1.276 1.277 ==== Develop algorithms that find genetic markers for anatomical regions ==== 1.278 # Develop scoring measures for evaluating how good individual genes are at marking areas: we will compare pointwise, geometric, and information-theoretic measures. 1.279 @@ -275,316 +227,15 @@ 1.280 1.281 1.282 1.283 - 1.284 - 1.285 - 1.286 - 1.287 -gradient similarity is calculated as: 1.288 -\sum_pixels cos(abs(\angle \nabla_1 - \angle \nabla_2)) \cdot \frac{\vert \nabla_1 \vert + \vert \nabla_2 \vert}{2} \cdot \frac{pixel\_value_1 + pixel\_value_2}{2} 1.289 - 1.290 - 1.291 - 1.292 - 1.293 - 1.294 - 1.295 - 1.296 -(todo) Technically, we say that an anatomical structure has a fundamentally 2-D organization when there exists a commonly used, generic, anatomical structure-preserving map from 3-D space to a 2-D manifold. 1.297 - 1.298 - 1.299 -Related work: 1.300 - 1.301 - 1.302 -The Allen Brain Institute has developed an interactive web interface called AGEA which allows an investigator to (1) calculate lists of genes which are selectively overexpressed in certain anatomical regions (ABA calls this the "Gene Finder" function) (2) to visualize the correlation between the genetic profiles of voxels in the dataset, and (3) to visualize a hierarchial clustering of voxels in the dataset \cite{ng_anatomic_2009}. AGEA is an impressive and useful tool, however, it does not solve the same problems that we propose to solve with this project. 1.303 - 1.304 -First we describe AGEA's "Gene Finder", and then compare it to our proposed method for finding marker genes. AGEA's Gene Finder first asks the investigator to select a single "seed voxel" of interest. It then uses a clustering method, combined with built-in knowledge of major anatomical structures, to select two sets of voxels; an "ROI" and a "comparator region"*. The seed voxel is always contained within the ROI, and the ROI is always contained within the comparator region. The comparator region is similar but not identical to the set of voxels making up the major anatomical region containing the ROI. Gene Finder then looks for genes which can distinguish the ROI from the comparator region. Specifically, it finds genes for which the ratio (expression energy in the ROI) / (expression energy in the comparator region) is high. 1.305 - 1.306 -Informally, the Gene Finder first infers an ROI based on clustering the seed voxel with other voxels. Then, the Gene Finder finds genes which overexpress in the ROI as compared to other voxels in the major anatomical region. 1.307 - 1.308 -There are three major differences between our approach and Gene Finder. 1.309 - 1.310 -First, Gene Finder focuses on individual genes and individual ROIs in isolation. This is great for regions which can be picked out from all other regions by a single gene, but not all of them can (todo). There are at least two ways this can miss out on useful genes. First, a gene might express in part of a region, but not throughout the whole region, but there may be another gene which expresses in the rest of the region*. Second, a gene might express in a region, but not in any of its neighbors, but it might express also in other non-neighboring regions. To take advantage of these types of genes, we propose to find combinations of genes which, together, can identify the boundaries of all subregions within the containing region. 1.311 - 1.312 -Second, Gene Finder uses a pointwise metric, namely expression energy ratio, to decide whether a gene is good for picking out a region. We have found better results by using metrics which take into account not just single voxels, but also the local geometry of neighboring voxels, such as the local gradient (todo). In addition, we have found that often the absence of gene expression can be used as a marker, which will not be caught by Gene Finder's expression energy ratio (todo). 1.313 - 1.314 -Third, Gene Finder chooses the ROI based only on the seed voxel. This often does not permit the user to query the ROI that they are interested in. For example, in all of our tests of Gene Finder in cortex, the ROIs chosen tend to be cortical layers, rather than cortical areas. 1.315 - 1.316 -In summary, when Gene Finder picks the ROI that you want, and when this ROI can be easily picked out from neighboring regions by single genes which selectively overexpress in the ROI compared to the entire major anatomical region, Gene Finder will work. However, Gene Finder will not pick cortical areas as ROIs, and even if it could, many cortical areas cannot be uniquely picked out by the overexpression of any single gene. By contrast, we will target cortical areas, we will explore a variety of metrics which can complement the shortcomings of expression energy ratio, and we will use the combinatorial expression of genes to pick out cortical areas even when no individual gene will do. 1.317 - 1.318 - 1.319 -* The terms "ROI" and "comparator region" are our own; the ABI calls them the "local region" and the "larger anatomical context". The ABI uses the term "specificity comparator" to mean the major anatomic region containing the ROI, which is not exactly identical to the comparator region. 1.320 - 1.321 -** In this case, the union of the area of expression of the two genes would suffice; one could also imagine that there could be situations in which the intersection of multiple genes would be needed, or a combination of unions and intersections. 1.322 - 1.323 - 1.324 -Now we describe AGEA's hierarchial clustering, and compare it to our proposal. The goal of AGEA's hierarchial clustering is to generate a binary tree of clusters, where a cluster is a collection of voxels. AGEA begins by computing the Pearson correlation between each pair of voxels. They then employ a recursive divisive (top-down) hierarchial clustering procedure on the voxels, which means that they start with all of the voxels, and then they divide them into clusters, and then within each cluster, they divide that cluster into smaller clusters, etc***. At each step, the collection of voxels is partitioned into two smaller clusters in a way that maximizes the following quantity: average correlation between all possible pairs of voxels containing one voxel from each cluster. 1.325 - 1.326 -There are three major differences between our approach and AGEA's hierarchial clustering. First, AGEA's clustering method separates cortical layers before it separates cortical areas. 1.327 - 1.328 - 1.329 - 1.330 - 1.331 - 1.332 -following procedure is used for the purpose of dividing a collection of voxels into smaller clusters: partition the voxels into two sets, such that the following quantity is maximized: 1.333 - 1.334 -*** depending on which level of the tree is being created, the voxels are subsampled in order to save time 1.335 - 1.336 - 1.337 - 1.338 - 1.339 - 1.340 -does not allow the user to input anything other than a seed voxel; this means that for each seed voxel, there is only one 1.341 - 1.342 - 1.343 - 1.344 -The role of the "local region" is to serve as a region of interest for which marker genes are desired; the role of the "larger anatomical context" is to be the structure 1.345 - 1.346 - 1.347 - 1.348 -There are two kinds of differences between AGEA and our project; differences that relate to the treatment of the cortex, and differences in the type of generalizable methods being developed. As relates 1.349 - 1.350 - 1.351 -indicate an ROI 1.352 - 1.353 -explore simple correlation-based relationships between voxels, genes, and clusters of voxels. 1.354 - 1.355 - 1.356 -There have not yet been any studies which describe the results of applying AGEA to the cerebral cortex; however, we suspect that the AGEA metrics are not optimal for the task of relating genes to cortical areas. A voxel's gene expression profile depends upon both its cortical area and its cortical layer, however, AGEA has no mechanism to distinguish these two. As a result, voxels in the same layer but different areas are often clustered together by AGEA. As part of the project, we will compare the performance of our techniques against AGEA's. 1.357 - 1.358 ---- 1.359 - 1.360 -The Allen Brain Institute has developed interactive tools called AGEA which allow an investigator to explore simple correlation-based relationships between voxels, genes, and clusters of voxels. There have not yet been any studies which describe the results of applying AGEA to the cerebral cortex; however, we suspect that the AGEA metrics are not optimal for the task of relating genes to cortical areas. A voxel's gene expression profile depends upon both its cortical area and its cortical layer, however, AGEA has no mechanism to distinguish these two. As a result, voxels in the same layer but different areas are often clustered together by AGEA. As part of the project, we will compare the performance of our techniques against AGEA's. 1.361 - 1.362 -Another difference between our techniques and AGEA's is that AGEA allows the user to enter only a voxel location, and then to either explore the rest of the brain's relationship to that particular voxel, or explore a partitioning of the brain based on pairwise voxel correlation. If the user is interested not in a single voxel, but rather an entire anatomical structure, AGEA will only succeed to the extent that the selected voxel is a typical representative of the structure. As discussed in the previous paragraph, this poses problems for structures like cortical areas, which (because of their division into cortical layers) do not have a single "typical representative". 1.363 - 1.364 -By contrast, in our system, the user will start by selecting, not a single voxel, but rather, an anatomical superstructure to be divided into pieces (for example, the cerebral cortex). We expect that our methods will take into account not just pairwise statistics between voxels, but also large-scale geometric features (for example, the rapidity of change in gene expression as regional boundaries are crossed) which optimize the discriminability of regions within the selected superstructure. 1.365 - 1.366 - 1.367 ------ 1.368 - 1.369 -screen for combinations of marker genes which selectively target anatomical structures 1.370 -pick delineate the boundaries between neighboring anatomical structures. 1.371 -(b) techniques to screen for marker genes which pick out anatomical structures of interest 1.372 - 1.373 -, techniques which: (a) screen for marker genes , and (b) suggest new anatomical maps based on 1.374 - 1.375 - 1.376 -whose expression partitions the region of interest into its anatomical substructures, and (b) use the natural contours of gene expression to suggest new ways of dividing an organ into 1.377 - 1.378 - 1.379 -The Allen Brain Atlas 1.380 - 1.381 - 1.382 - 1.383 - 1.384 --- 1.385 - 1.386 -to: brooksl@mail.nih.gov 1.387 - 1.388 -Hi, I'm writing to confirm the applicability of a potential research 1.389 -project to the challenge grant topic "New computational and 1.390 -statistical methods for the analysis of large 1.391 -data sets from next-generation sequencing technologies". 1.392 - 1.393 -We want to develop methods for the analysis of gene expression 1.394 -datasets that can be used to uncover the relationships between gene 1.395 -expression and anatomical regions. Specifically, we want to develop 1.396 -techniques to (a) given a set of known anatomical areas, identify 1.397 -genetic markers for each of these areas, and (b) given an anatomical structure 1.398 -whose substructure is unknown, suggest a map, that is, a division of 1.399 -the space into anatomical sub-structures, that represents the 1.400 -boundaries inherent in the gene expression data. 1.401 - 1.402 -We propose to develop our techniques on the Allen Brain 1.403 -Atlas mouse brain gene expression dataset by finding genetic markers 1.404 -for anatomical areas within the cerebral cortex. The Allen Brain Atlas 1.405 -contains a registered 3-D map of gene expression data with 200-micron 1.406 -voxel resolution which was created from in situ hybridization 1.407 -data. The dataset contains about 4000 genes which are available at 1.408 -this resolution across the entire cerebral cortex. 1.409 - 1.410 -Despite the distinct roles of different cortical 1.411 -areas in both normal functioning and disease processes, there are no 1.412 -known marker genes for many cortical areas. This project will be 1.413 -immediately useful for both drug discovery and clinical research 1.414 -because once the markers are known, interventions can be designed 1.415 -which selectively target specific cortical areas. 1.416 - 1.417 -This techniques we develop will be useful because they will be 1.418 -applicable to the analysis of other anatomical areas, both in 1.419 -terms of finding marker genes for known areas, and in terms of 1.420 -suggesting new anatomical subdivisions that are based upon the gene 1.421 -expression data. 1.422 - 1.423 - 1.424 - 1.425 ---- 1.426 1.427 - 1.428 - 1.429 - 1.430 - 1.431 - 1.432 -It is likely that our study, by showing which areal divisions naturally follow from gene expression data, as opposed to traditional histological data, will contribute to the creation of 1.433 - 1.434 -there are clear genetic or chemical markers known for only a few cortical areas. This makes it difficult to target drugs to specific 1.435 - 1.436 -As part of aims (1) and (5), we will discover sets of marker genes that pick out specific cortical areas. This will allow the development of drugs and other interventions that selectively target individual cortical areas. As part of aims (2) and (5), we will also discover small panels of marker genes that can be used to delineate most of the cortical areal map. 1.437 - 1.438 - 1.439 - 1.440 -With aims (2) and (4), we 1.441 - 1.442 -There are five principals 1.443 - 1.444 - 1.445 - 1.446 -In addition to validating the usefulness of the algorithms, the application of these methods to cerebral cortex will produce immediate benefits that are only one step removed from clinical application. 1.447 - 1.448 - 1.449 -todo: remember to check gensat, etc for validation (mention bias/variance) 1.450 - 1.451 - 1.452 - 1.453 -=== Why it is useful to apply these methods to cortex === 1.454 - 1.455 - 1.456 -There is still room for debate as to exactly how the cortex should be parcellated into areas. 1.457 - 1.458 - 1.459 -The best way to divide up rodent cortex into areas has not been completely determined, 1.460 - 1.461 - 1.462 -not yet been accounted for in 1.463 - 1.464 - that the expression of some genes will contain novel spatial patterns which are not account 1.465 - 1.466 - that a genoarchitectonic map 1.467 - 1.468 - 1.469 -This principle is only applicable to aim 1 (marker genes). For aim 2 (partition a structure in into anatomical subregions), we plan to work with many genes at once. 1.470 - 1.471 - 1.472 -tood: aim 2 b+s? 1.473 - 1.474 - 1.475 - 1.476 - 1.477 -==== Principle 5: Interoperate with existing tools ==== 1.478 - 1.479 -In order for our software to be as useful as possible for our users, it will be able to import and export data to standard formats so that users can use our software in tandem with other software tools created by other teams. We will support the following formats: NIFTI (Neuroimaging Informatics Technology Initiative), SEV (Allen Brain Institute Smoothed Energy Volume), and MATLAB. This ensures that our users will not have to exclusively rely on our tools when analyzing data. For example, users will be able to use the data visualization and analysis capabilities of MATLAB and Caret alongside our software. 1.480 - 1.481 -To our knowledge, there is no currently available software to convert between these formats, so we will also provide a format conversion tool. This may be useful even for groups that don't use any of our other software. 1.482 - 1.483 - 1.484 - 1.485 -todo: is "marker gene" even a phrase that we should use at all? 1.486 - 1.487 - 1.488 - 1.489 -note for aim 1 apps: combo of genes is for voxel, not within any single cell 1.490 - 1.491 - 1.492 - 1.493 - 1.494 -, as when genetic markers allow the development of selective interventions; the reason that one can be confident that the intervention is selective is that it is only turned on when a certain combination of genes is turned on and off. The result procedure is what assures us that when that combination is present, the local tissue is probably part of a certain subregion. 1.495 - 1.496 - 1.497 - 1.498 -The basic idea is that we want to find a procedure by 1.499 - 1.500 -The task of finding genes that mark anatomical areas can be phrased in terms of what the field of machine learning calls a "supervised learning" task. The goal of this task is to learn a function (the "classifier") which 1.501 - 1.502 -If a person knows a combination of genes that mark an area, that implies that the person can be told how strong those genes express in any voxel, and the person can use this information to determine how 1.503 - 1.504 -finding how to infer the areal identity of a voxel if given the gene expression profile of that voxel. 1.505 - 1.506 - 1.507 -For each voxel in the cortex, we want to start with data about the gene expression 1.508 - 1.509 - 1.510 - 1.511 -There are various ways to look for marker genes. We will define some terms, and along the way we will describe a few design choices encountered in the process of creating a marker gene finding method, and then we will present four principles that describe which options we have chosen. 1.512 - 1.513 -In developing a procedure for finding marker genes, we are developing a procedure that takes a dataset of experimental observations and produces a result. One can think of the result as merely a list of genes, but really the result is an understanding of a predictive relationship between, on the one hand, the expression levels of genes, and, on the other hand, anatomical subregions. 1.514 - 1.515 -One way to more formally define this understanding is to look at it as a procedure. In this view, the result of the learning procedure is itself a procedure. The result procedure provides a way to use the gene expression profiles of voxels in a tissue sample in order to determine where the subregions are. 1.516 - 1.517 -This result procedure can be used directly, as when an experimenter has a tissue sample and needs to know what subregions are present in it, and, if multiple subregions are present, where they each are. Or it can be used indirectly; imagine that the result procedure tells us that whenever a certain combination of genes are expressed, the local tissue is probably part of a certain subregion. This means that we can then confidentally develop an intervention which is triggered only when that combination of genes are expressed; and to the extent that the result procedure is reliable, we know that the intervention will only be triggered in the target subregion. 1.518 - 1.519 -We said that the result procedure provides "a way to use the gene expression profiles of voxels in a tissue sample" in order to "determine where the subregions are". 1.520 - 1.521 - 1.522 -Does the result procedure get as input all of the gene expression profiles of each voxel in the entire tissue sample, and produce as output all of the subregional boundaries all at once? 1.523 - 1.524 - 1.525 - 1.526 - 1.527 - 1.528 - 1.529 - it is helpful for the classifier to look at the global "shape" of gene expression patterns over the whole structure, rather than just nearby voxels. 1.530 - 1.531 - 1.532 - 1.533 - 1.534 -there is some small bit of additional information that can be gleaned from knowing the 1.535 - 1.536 -==== Design choices for a supervised learning procedure ==== 1.537 - 1.538 - 1.539 -After all, 1.540 - 1.541 -there is a small correlation between the gene expression levels from distant voxels and 1.542 - 1.543 -Depending on how we intend to use the classifier, we may want to design it so that 1.544 - 1.545 -It is possible for many things to 1.546 - 1.547 -The choice of which data is made part of an instance 1.548 - 1.549 -what we seek is a procedure 1.550 - 1.551 - partition the tissue sample into subregions. 1.552 - 1.553 - each part of the anatomical structure 1.554 - 1.555 - must be One way to rephrase this task is to say that, instead of searching for the location of the subregions, we are looking to partition the tissue sample into subregions. 1.556 - 1.557 - 1.558 -There are various ways to look for marker genes. We will define some terms, and along the way we will describe a few design choices encountered in the process of creating a marker gene finding method, and then we will present four principles that describe which options we have chosen. 1.559 - 1.560 -In developing a procedure for finding marker genes, we are developing a procedure that takes a dataset of experimental observations and produces a result. One can think of the result as merely a list of genes, but really the result is an understanding of a predictive relationship between, on the one hand, the expression levels of genes, and, on the other hand, anatomical subregions. 1.561 - 1.562 -One way to more formally define this understanding is to look at it as a procedure. In this view, the result of the learning procedure is itself a procedure. The result procedure provides a way to use the gene expression profiles of voxels in a tissue sample in order to determine where the subregions are. 1.563 - 1.564 -This result procedure can be used directly, as when an experimenter has a tissue sample and needs to know what subregions are present in it, and, if multiple subregions are present, where they each are. Or it can be used indirectly; imagine that the result procedure tells us that whenever a certain combination of genes are expressed, the local tissue is probably part of a certain subregion. This means that we can then confidentally develop an intervention which is triggered only when that combination of genes are expressed; and to the extent that the result procedure is reliable, we know that the intervention will only be triggered in the target subregion. 1.565 - 1.566 -We said that the result procedure provides "a way to use the gene expression profiles of voxels in a tissue sample" in order to "determine where the subregions are". 1.567 - 1.568 - 1.569 -Does the result procedure get as input all of the gene expression profiles of each voxel in the entire tissue sample, and produce as output all of the subregional boundaries all at once? 1.570 - 1.571 - 1.572 -Or are we given one voxel at a time, 1.573 - 1.574 - 1.575 -In the jargon of the field of machine learning, the result procedure is called a __classifier__. 1.576 - 1.577 - 1.578 -The task of finding genes that mark anatomical areas can be phrased in terms of what the field of machine learning calls a "supervised learning" task. The goal of this task is to learn a function (the "classifier") which 1.579 - 1.580 -If a person knows a combination of genes that mark an area, that implies that the person can be told how strong those genes express in any voxel, and the person can use this information to determine how 1.581 - 1.582 -finding how to infer the areal identity of a voxel if given the gene expression profile of that voxel. 1.583 - 1.584 - 1.585 -For each voxel in the cortex, we want to start with data about the gene expression 1.586 - 1.587 - 1.588 - 1.589 -single voxels, but rather groups of voxels, such that the groups can be placed in some 2-D space. We will call such instances "pixels". 1.590 - 1.591 -We have been speaking as if instances necessarily correspond to single voxels. But it is possible for instances to be groupings of many voxels, in which case each grouping must be assigned the same label (that is, each voxel grouping must stay inside a single anatomical subregion). 1.592 - 1.593 - 1.594 - 1.595 -In some but not all cases, the groups are either rows or columns of voxels. This is the case with the cerebral cortex, in which one may assume that columns of voxels which run perpendicular to the cortical surface all share the same areal identity. In the cortex, we call such an instance a "surface pixel", because such an instance represents the data associated with all voxels underneath a specific patch of the cortical surface. 1.596 +stuff i dunno where to put yet (there is more scattered through grant-oldtext): 1.597 + 1.598 +==== Principle 4: Work in 2-D whenever possible ==== 1.599 + 1.600 +In anatomy, the manifold of interest is usually either defined by a combination of two relevant anatomical axes (todo), or by the surface of the structure (as is the case with the cortex). In the former case, the manifold of interest is a plane, but in the latter case it is curved. If the manifold is curved, there are various methods for mapping the manifold into a plane. 1.601 + 1.602 +The method that we will develop will begin by mapping the data into a 2-D plane. Although the manifold that characterized cortical areas is known to be the cortical surface, it remains to be seen which method of mapping the manifold into a plane is optimal for this application. We will compare mappings which attempt to preserve size (such as the one used by Caret\ref{van_essen_integrated_2001}) with mappings which preserve angle (conformal maps). 1.603 + 1.604 +Although there is much 2-D organization in anatomy, there are also structures whose shape is fundamentally 3-dimensional. If possible, we would like the method we develop to include a statistical test that warns the user if the assumption of 2-D structure seems to be wrong. 1.605 +