1 Specific aims

2 Massive new datasets obtained with techniques such as in situ hybridization

3 (ISH) and BAC-transgenics allow the expression levels of many genes at many

4 locations to be compared. Our goal is to develop automated methods to relate

5 spatial variation in gene expression to anatomy. We want to find marker genes

6 for specific anatomical regions, and also to draw new anatomical maps based on

7 gene expression patterns. We have three specific aims:

8 (1) develop an algorithm to screen spatial gene expression data for combina-

9 tions of marker genes which selectively target anatomical regions

10 (2) develop an algorithm to suggest new ways of carving up a structure into

11 anatomical subregions, based on spatial patterns in gene expression

12 (3) create a 2-D “flat map” dataset of the mouse cerebral cortex that contains

13 a flattened version of the Allen Mouse Brain Atlas ISH data, as well as

14 the boundaries of cortical anatomical areas. Use this dataset to validate

15 the methods developed in (1) and (2).

16 In addition to validating the usefulness of the algorithms, the application of

17 these methods to cerebral cortex will produce immediate benefits, because there

18 are currently no known genetic markers for many cortical areas. The results

19 of the project will support the development of new ways to selectively target

20 cortical areas, and it will support the development of a method for identifying

21 the cortical areal boundaries present in small tissue samples.

22 All algorithms that we develop will be implemented in an open-source soft-

23 ware toolkit. The toolkit, as well as the machine-readable datasets developed

24 in aim (3), will be published and freely available for others to use.

25 Background and significance

26 Aim 1

27 Machine learning terminology

28 The task of looking for marker genes for anatomical subregions means that one

29 is looking for a set of genes such that, if the expression level of those genes is

30 known, then the locations of the subregions can be inferred.

31 If we define the subregions so that they cover the entire anatomical structure

32 to be divided, then instead of saying that we are using gene expression to find

33 the locations of the subregions, we may say that we are using gene expression to

34 determine to which subregion each voxel within the structure belongs. We call

35 this a classification task, because each voxel is being assigned to a class (namely,

36 its subregion).

37 Therefore, an understanding of the relationship between the combination of

38 their expression levels and the locations of the subregions may be expressed as

39 1

40

41 a function. The input to this function is a voxel, along with the gene expression

42 levels within that voxel; the output is the subregional identity of the target

43 voxel, that is, the subregion to which the target voxel belongs. We call this

44 function a classifier. In general, the input to a classifier is called an instance,

45 and the output is called a label.

46 The object of aim 1 is not to produce a single classifier, but rather to develop

47 an automated method for determining a classifier for any known anatomical

48 structure. Therefore, we seek a procedure by which a gene expression dataset

49 may be analyzed in concert with an anatomical atlas in order to produce a

50 classifier. Such a procedure is a type of a machine learning procedure. The

51 construction of the classifier is called training (also learning), and the initial

52 gene expression dataset used in the construction of the classifier is called training

53 data.

54 In the machine learning literature, this sort of procedure may be thought

55 of as a supervised learning task, defined as a task in whcih the goal is to learn

56 a mapping from instances to labels, and the training data consists of a set of

57 instances (voxels) for which the labels (subregions) are known.

58 Each gene expression level is called a feature, and the selection of which

59 genes to include is called feature selection. Feature selection is one component

60 of the task of learning a classifier. Some methods for learning classifiers start

61 out with a separate feature selection phase, whereas other methods combine

62 feature selection with other aspects of training.

63 One class of feature selection methods assigns some sort of score to each

64 candidate gene. The top-ranked genes are then chosen. Some scoring measures

65 can assign a score to a set of selected genes, not just to a single gene; in this

66 case, a dynamic procedure may be used in which features are added and sub-

67 tracted from the selected set depending on how much they raise the score. Such

68 procedures are called “stepwise” or “greedy”.

69 Although the classifier itself may only look at the gene expression data within

70 each voxel before classifying that voxel, the learning algorithm which constructs

71 the classifier may look over the entire dataset. We can categorize score-based

72 feature selection methods depending on how the score of calculated. Often

73 the score calculation consists of assigning a sub-score to each voxel, and then

74 aggregating these sub-scores into a final score (the aggregation is often a sum or

75 a sum of squares). If only information from nearby voxels is used to calculate a

76 voxel’s sub-score, then we say it is a local scoring method. If only information

77 from the voxel itself is used to calculate a voxel’s sub-score, then we say it is a

78 pointwise scoring method.

79 Key questions when choosing a learning method are: What are the instances?

80 What are the features? How are the features chosen? Here are four principles

81 that outline our answers to these questions.

82 Principle 1: Combinatorial gene expression

83 Above, we defined an “instance” as the combination of a voxel with the “asso-

84 ciated gene expression data”. In our case this refers to the expression level of

85 2

86

87 genes within the voxel, but should we include the expression levels of all genes,

88 or only a few of them?

89 It is too much to hope that every anatomical region of interest will be iden-

90 tified by a single gene. For example, in the cortex, there are some areas which

91 are not clearly delineated by any gene included in the Allen Brain Atlas (ABA)

92 dataset. However, at least some of these areas can be delineated by looking

93 at combinations of genes (an example of an area for which multiple genes are

94 necessary and sufficient is provided in Preliminary Results).

95 Principle 2: Only look at combinations of small numbers of genes

96 When the classifier classifies a voxel, it is only allowed to look at the expression of

97 the genes which have been selected as features. The more data that is available

98 to a classifier, the better that it can do. For example, perhaps there are weak

99 correlations over many genes that add up to a strong signal. So, why not include

100 every gene as a feature? The reason is that we wish to employ the classifier in

101 situations in which it is not feasible to gather data about every gene. For

102 example, if we want to use the expression of marker genes as a trigger for some

103 regionally-targeted intervention, then our intervention must contain a molecular

104 mechanism to check the expression level of each marker gene before it triggers.

105 It is currently infeasible to design a molecular trigger that checks the level of

106 more than a handful of genes. Similarly, if the goal is to develop a procedure to

107 do ISH on tissue samples in order to label their anatomy, then it is infeasible

108 to label more than a few genes. Therefore, we must select only a few genes as

109 features.

110 Principle 3: Use geometry in feature selection

111 When doing feature selection with score-based methods, the simplest thing to

112 do would be to score the performance of each voxel by itself and then combine

113 these scores; this is pointwise scoring. A more powerful approach is to also use

114 information about the geometric relations between each voxel and its neighbors;

115 this requires non-pointwise, local scoring methods. See Preliminary Results for

116 evidence of the complementary nature of pointwise and local scoring methods.

117 Principle 4: Work in 2-D whenever possible

118 There are many anatomical structures which are commonly characterized in

119 terms of a two-dimensional manifold. When it is known that the structure that

120 one is looking for is two-dimensional, the results may be improved by allowing

121 the analysis algorithm to take advantage of this prior knowledge. In addition,

122 it is easier for humans to visualize and work with 2-D data.

123 Therefore, when possible, the instances should represent pixels, not voxels.

124 3

125

126 Aim 3

127 Background

128 The cortex is divided into areas and layers. To a first approximation, the par-

129 cellation of the cortex into areas can be drawn as a 2-D map on the surface

130 of the cortex. In the third dimension, the boundaries between the areas con-

131 tinue downwards into the cortical depth, perpendicular to the surface. The layer

132 boundaries run parallel to the surface. One can picture an area of the cortex as

133 a slice of many-layered cake.

134 Although it is known that different cortical areas have distinct roles in both

135 normal functioning and in disease processes, there are no known marker genes

136 for many cortical areas. When it is necessary to divide a tissue sample into

137 cortical areas, this is a manual process that requires a skilled human to combine

138 multiple visual cues and interpret them in the context of their approximate

139 location upon the cortical surface.

140 Even the questions of how many areas should be recognized in cortex, and

141 what their arrangement is, are still not completely settled. A proposed division

142 of the cortex into areas is called a cortical map. In the rodent, the lack of a

143 single agreed-upon map can be seen by contrasting the recent maps given by

144 Swanson?? on the one hand, and Paxinos and Franklin?? on the other. While

145 the maps are certainly very similar in their general arrangement, significant

146 differences remain in the details.

147 Significance

148 The method developed in aim (1) will be applied to each cortical area to find

149 a set of marker genes such that the combinatorial expression pattern of those

150 genes uniquely picks out the target area. Finding marker genes will be useful

151 for drug discovery as well as for experimentation because marker genes can be

152 used to design interventions which selectively target individual cortical areas.

153 The application of the marker gene finding algorithm to the cortex will

154 also support the development of new neuroanatomical methods. In addition to

155 finding markers for each individual cortical areas, we will find a small panel

156 of genes that can find many of the areal boundaries at once. This panel of

157 marker genes will allow the development of an ISH protocol that will allow

158 experimenters to more easily identify which anatomical areas are present in

159 small samples of cortex.

160 The method developed in aim (3) will provide a genoarchitectonic viewpoint

161 that will contribute to the creation of a better map. The development of present-

162 day cortical maps was driven by the application of histological stains. It is

163 conceivable that if a different set of stains had been available which identified

164 a different set of features, then the today’s cortical maps would have come out

165 differently. Since the number of classes of stains is small compared to the number

166 of genes, it is likely that there are many repeated, salient spatial patterns in

167 the gene expression which have not yet been captured by any stain. Therefore,

168 4

169

170 current ideas about cortical anatomy need to incorporate what we can learn

171 from looking at the patterns of gene expression.

172 While we do not here propose to analyze human gene expression data, it is

173 conceivable that the methods we propose to develop could be used to suggest

174 modifications to the human cortical map as well.

175 Related work

176 Preliminary work

177 Justification of principles 1 thur 3

178 Principle 1: Combinatorial gene expression

179 Here we give an example of a cortical area which is not marked by any single

180 gene, but which can be identified combinatorially. according to logistic regres-

181 sion, gene wwc11 is the best fit single gene for predicting whether or not a pixel

182 on the cortical surface belongs to the motor area (area MO). The upper-left

183 picture in Figure shows wwc1’s spatial expression pattern over the cortex. The

184 lower-right boundary of MO is represented reasonably well by this gene, however

185 the gene overshoots the upper-left boundary. This flattened 2-D representation

186 does not show it, but the area corresponding to the overshoot is the medial

187 surface of the cortex. MO is only found on the lateral surface (todo).

188 Gnee mtif22 is shown in figure the upper-right of Fig. . Mtif2 captures MO’s

189 upper-left boundary, but not its lower-right boundary. Mtif2 does not express

190 very much on the medial surface. By adding together the values at each pixel

191 in these two figures, we get the lower-left of Figure . This combination captures

192 area MO much better than any single gene.

193 Principle 2: Only look at combinations of small numbers of genes

194 In order to see how well one can do when looking at all genes at once, we ran

195 a support vector machine to classify cortical surface pixels based on their gene

196 expression profiles. We achieved classification accuracy of about 81%3. As noted

197 above, however, a classifier that looks at all the genes at once isn’t practically

198 useful.

199 The requirement to find combinations of only a small number of genes limits

200 us from straightforwardly applying many of the most simple techniques from

201 the field of supervised machine learning. In the parlance of machine learning,

202 our task combines feature selection with supervised learning.

203 __________________________

204 1“WW, C2 and coiled-coil domain containing 1”; EntrezGene ID 211652

205 2“mitochondrial translational initiation factor 2”; EntrezGene ID 76784

206 3Using the Shogun SVM package (todo:cite), with parameters type=GMNPSVM (multi-

207 class b-SVM), kernal = gaussian with sigma = 0.1, c = 10, epsilon = 1e-1 – these are the

208 first parameters we tried, so presumably performance would improve with different choices of

209 parameters. 5-fold cross-validation.

210 5

211

212

213

214 Figure 1: Upper left: wwc1. Upper right: mtif2. Lower left: wwc1 + mtif2

215 (each pixel’s value on the lower left is the sum of the corresponding pixels in

216 the upper row). Within each picture, the vertical axis roughly corresponds to

217 anterior at the top and posterior at the bottom, and the horizontal axis roughly

218 corresponds to medial at the left and lateral at the right. The red outline is

219 the boundary of region MO. Pixels are colored approximately according to the

220 density of expressing cells underneath each pixel, with red meaning a lot of

221 expression and blue meaning little.

222 6

223

224

225

226 Figure 2: The top row shows the three genes which (individually) best predict

227 area AUD, according to logistic regression. The bottom row shows the three

228 genes which (individually) best match area AUD, according to gradient similar-

229 ity. From left to right and top to bottom, the genes are Ssr1, Efcbp1, Aph1a,

230 Ptk7, Aph1a again, and Lepr

231 Principle 3: Use geometry

232 To show that local geometry can provide useful information that cannot be

233 detected via pointwise analyses, consider Fig. . The top row of Fig. displays

234 the 3 genes which most match area AUD, according to a pointwise method4. The

235 bottom row displays the 3 genes which most match AUD according to a method

236 which considers local geometry5 The pointwise method in the top row identifies

237 genes which express more strongly in AUD than outside of it; its weakness is that

238 this includes many areas which don’t have a salient border matching the areal

239 border. The geometric method identifies genes whose salient expression border

240 seems to partially line up with the border of AUD; its weakness is that this

241 includes genes which don’t express over the entire area. Genes which have high

242 rankings using both pointwise and border criteria, such as Aph1a in the example,

243 may be particularly good markers. None of these genes are, individually, a

244 perfect marker for AUD; we deliberately chose a “difficult” area in order to

245 better contrast pointwise with geometric methods.

246 __________________________

247 4For each gene, a logistic regression in which the response variable was whether or not a

248 surface pixel was within area AUD, and the predictor variable was the value of the expression

249 of the gene underneath that pixel. The resulting scores were used to rank the genes in terms

250 of how well they predict area AUD.

251 5For each gene the gradient similarity (see section ??) between (a) a map of the expression

252 of each gene on the cortical surface and (b) the shape of area AUD, was calculated, and this

253 was used to rank the genes.

254 7

255

256 Principle 4: Work in 2-D whenever possible

257 In anatomy, the manifold of interest is usually either defined by a combination

258 of two relevant anatomical axes (todo), or by the surface of the structure (as is

259 the case with the cortex). In the former case, the manifold of interest is a plane,

260 but in the latter case it is curved. If the manifold is curved, there are various

261 methods for mapping the manifold into a plane.

262 The method that we will develop will begin by mapping the data into a

263 2-D plane. Although the manifold that characterized cortical areas is known

264 to be the cortical surface, it remains to be seen which method of mapping the

265 manifold into a plane is optimal for this application. We will compare mappings

266 which attempt to preserve size (such as the one used by Caret??) with mappings

267 which preserve angle (conformal maps).

268 Although there is much 2-D organization in anatomy, there are also struc-

269 tures whose shape is fundamentally 3-dimensional. If possible, we would like

270 the method we develop to include a statistical test that warns the user if the

271 assumption of 2-D structure seems to be wrong.

272 ——

273 Massive new datasets obtained with techniques such as in situ hybridization

274 (ISH) and BAC-transgenics allow the expression levels of many genes at many

275 locations to be compared. This can be used to find marker genes for specific

276 anatomical structures, as well as to draw new anatomical maps. Our goal is

277 to develop automated methods to relate spatial variation in gene expression to

278 anatomy. We have five specific aims:

279 (1) develop an algorithm to screen spatial gene expression data for combi-

280 nations of marker genes which selectively target individual anatomical

281 structures

282 (2) develop an algorithm to screen spatial gene expression data for combina-

283 tions of marker genes which can be used to delineate most of the bound-

284 aries between a number of anatomical structures at once

285 (3) develop an algorithm to suggest new ways of dividing a structure up into

286 anatomical subregions, based on spatial patterns in gene expression

287 (4) create a flat (2-D) map of the mouse cerebral cortex that contains a flat-

288 tened version of the Allen Mouse Brain Atlas ISH dataset, as well as the

289 boundaries of anatomical areas within the cortex. For each cortical layer,

290 a layer-specific flat dataset will be created. A single combined flat dataset

291 will be created which averages information from all of the layers. These

292 datasets will be made available in both MATLAB and Caret formats.

293 (5) validate the methods developed in (1), (2) and (3) by applying them to

294 the cerebral cortex datasets created in (4)

295 All algorithms that we develop will be implemented in an open-source soft-

296 ware toolkit. The toolkit, as well as the machine-readable datasets developed in

297 8

298

299 aim (4) and any other intermediate dataset we produce, will be published and

300 freely available for others to use.

301 In addition to developing generally useful methods, the application of these

302 methods to cerebral cortex will produce immediate benefits that are only one

303 step removed from clinical application, while also supporting the development

304 of new neuroanatomical techniques. The method developed in aim (1) will be

305 applied to each cortical area to find a set of marker genes. Currently, despite

306 the distinct roles of different cortical areas in both normal functioning and

307 disease processes, there are no known marker genes for many cortical areas.

308 Finding marker genes will be immediately useful for drug discovery as well as for

309 experimentation because once marker genes for an area are known, interventions

310 can be designed which selectively target that area.

311 The method developed in aim (2) will be used to find a small panel of genes

312 that can find most of the boundaries between areas in the cortex. Today, finding

313 cortical areal boundaries in a tissue sample is a manual process that requires a

314 skilled human to combine multiple visual cues over a large area of the cortical

315 surface. A panel of marker genes will allow the development of an ISH protocol

316 that will allow experimenters to more easily identify which anatomical areas are

317 present in small samples of cortex.

318 For each cortical layer, a layer-specific flat dataset will be created. A single

319 combined flat dataset will be created which averages information from all of

320 the layers. These datasets will be made available in both MATLAB and Caret

321 formats.

322 —-

323 New techniques allow the expression levels of many genes at many locations

324 to be compared. It is thought that even neighboring anatomical structures have

325 different gene expression profiles. We propose to develop automated methods

326 to relate the spatial variation in gene expression to anatomy. We will develop

327 two kinds of techniques:

328 (a) techniques to screen for combinations of marker genes which selectively

329 target anatomical structures

330 (b) techniques to suggest new ways of dividing a structure up into anatomical

331 subregions, based on the shapes of contours in the gene expression

332 The first kind of technique will be helpful for finding marker genes associated

333 with known anatomical features. The second kind of technique will be helpful in

334 creating new anatomical maps, maps which reflect differences in gene expression

335 the same way that existing maps reflect differences in histology.

336 We intend to develop our techniques using the adult mouse cerebral cortex

337 as a testbed. The Allen Brain Atlas has collected a dataset containing the

338 expression level of about 4000 genes* over a set of over 150000 voxels, with a

339 spatial resolution of approximately 200 microns[?].

340 We expect to discover sets of marker genes that pick out specific cortical

341 areas. This will allow the development of drugs and other interventions that

342 selectively target individual cortical areas. Therefore our research will lead

343 9

344

345 to application in drug discovery, in the development of other targeted clinical

346 interventions, and in the development of new experimental techniques.

347 The best way to divide up rodent cortex into areas has not been completely

348 determined, as can be seen by the differences in the recent maps given by Swan-

349 son on the one hand, and Paxinos and Franklin on the other. It is likely that our

350 study, by showing which areal divisions naturally follow from gene expression

351 data, as opposed to traditional histological data, will contribute to the creation

352 of a better map. While we do not here propose to analyze human gene expres-

353 sion data, it is conceivable that the methods we propose to develop could be

354 used to suggest modifications to the human cortical map as well.

355 In the following, we will only be talking about coronal data.

356 The Allen Brain Atlas provides “Smoothed Energy Volumes”, which are

357 One type of artifact in the Allen Brain Atlas data is what we call a “slice

358 artifact”. We have noticed two types of slice artifacts in the dataset. The first

359 type, a “missing slice artifact”, occurs when the ISH procedure on a slice did

360 not come out well. In this case, the Allen Brain investigators excluded the slice

361 at issue from the dataset. This means that no gene expression information is

362 available for that gene for the region of space covered by that slice. This results

363 in an expression level of zero being assigned to voxels covered by the slice. This

364 is partially but not completely ameliorated by the smoothing that is applied to

365 create the Smoothed Energy Volumes. The usual end result is that a region of

366 space which is shaped and oriented like a coronal slice is marked as having less

367 gene expression than surrounding regions.

368 The second type of slice artifact is caused by the fact that all of the slices

369 have a consistent orientation. Since there may be artifacts (such as how well

370 the ISH worked) which are constant within each slice but which vary between

371 different slices, the result is that ceteris paribus, when one compares the genetic

372 data of a voxel to another voxel within the same coronal plane, one would expect

373 to find more similarity than if one compared a voxel to another voxel displaced

374 along the rostrocaudal axis.

375 We are enthusiastic about the sharing of methods, data, and results, and

376 at the conclusion of the project, we will make all of our data and computer

377 source code publically available. Our goal is that replicating our results, or

378 applying the methods we develop to other targets, will be quick and easy for

379 other investigators. In order to aid in understanding and replicating our results,

380 we intend to include a software program which, when run, will take as input

381 the Allen Brain Atlas raw data, and produce as output all numbers and charts

382 found in publications resulting from the project.

383 To aid in the replication of our results, we will include a script which takes

384 as input the dataset in aim (3) and provides as output all of the tables in figures

385 in our publications .

386 We also expect to weigh in on the debate about how to best partition rodent

387 cortex

388 be useful for drug discovery as well

389 * Another 16000 genes are available, but they do not cover the entire cerebral

390 cortex with high spatial resolution.

391 10

392

393 User-definable ROIs Combinatorial gene expression Negative as well as pos-

394 itive signal Use geometry Search for local boundaries if necessary Flatmapped

395 Specific aims

396 Develop algorithms that find genetic markers for anatomical regions

397 1. Develop scoring measures for evaluating how good individual genes are at

398 marking areas: we will compare pointwise, geometric, and information-

399 theoretic measures.

400 2. Develop a procedure to find single marker genes for anatomical regions: for

401 each cortical area, by using or combining the scoring measures developed,

402 we will rank the genes by their ability to delineate each area.

403 3. Extend the procedure to handle difficult areas by using combinatorial cod-

404 ing: for areas that cannot be identified by any single gene, identify them

405 with a handful of genes. We will consider both (a) algorithms that incre-

406 mentally/greedily combine single gene markers into sets, such as forward

407 stepwise regression and decision trees, and also (b) supervised learning

408 techniques which use soft constraints to minimize the number of features,

409 such as sparse support vector machines.

410 4. Extend the procedure to handle difficult areas by combining or redrawing

411 the boundaries: An area may be difficult to identify because the bound-

412 aries are misdrawn, or because it does not “really” exist as a single area,

413 at least on the genetic level. We will develop extensions to our procedure

414 which (a) detect when a difficult area could be fit if its boundary were

415 redrawn slightly, and (b) detect when a difficult area could be combined

416 with adjacent areas to create a larger area which can be fit.

417 Apply these algorithms to the cortex

418 1. Create open source format conversion tools: we will create tools to bulk

419 download the ABA dataset and to convert between SEV, NIFTI and MAT-

420 LAB formats.

421 2. Flatmap the ABA cortex data: map the ABA data onto a plane and draw

422 the cortical area boundaries onto it.

423 3. Find layer boundaries: cluster similar voxels together in order to auto-

424 matically find the cortical layer boundaries.

425 4. Run the procedures that we developed on the cortex: we will present, for

426 each area, a short list of markers to identify that area; and we will also

427 present lists of “panels” of genes that can be used to delineate many areas

428 at once.

429 11

430

431 Develop algorithms to suggest a division of a structure into anatom-

432 ical parts

433 1. Explore dimensionality reduction algorithms applied to pixels: including

434 TODO

435 2. Explore dimensionality reduction algorithms applied to genes: including

436 TODO

437 3. Explore clustering algorithms applied to pixels: including TODO

438 4. Explore clustering algorithms applied to genes: including gene shaving,

439 TODO

440 5. Develop an algorithm to use dimensionality reduction and/or hierarchial

441 clustering to create anatomical maps

442 6. Run this algorithm on the cortex: present a hierarchial, genoarchitectonic

443 map of the cortex

444 gradient similarity is calculated as: ∑

445 pixels cos(abs(∠∇1 - ∠∇2)) ⋅|∇1|+|∇2|

446 2 ⋅

447 pixel_value1+pixel_value2

448 2

449 (todo) Technically, we say that an anatomical structure has a fundamen-

450 tally 2-D organization when there exists a commonly used, generic, anatomical

451 structure-preserving map from 3-D space to a 2-D manifold.

452 Related work:

453 The Allen Brain Institute has developed an interactive web interface called

454 AGEA which allows an investigator to (1) calculate lists of genes which are se-

455 lectively overexpressed in certain anatomical regions (ABA calls this the “Gene

456 Finder” function) (2) to visualize the correlation between the genetic profiles of

457 voxels in the dataset, and (3) to visualize a hierarchial clustering of voxels in

458 the dataset [?]. AGEA is an impressive and useful tool, however, it does not

459 solve the same problems that we propose to solve with this project.

460 First we describe AGEA’s “Gene Finder”, and then compare it to our pro-

461 posed method for finding marker genes. AGEA’s Gene Finder first asks the

462 investigator to select a single “seed voxel” of interest. It then uses a clustering

463 method, combined with built-in knowledge of major anatomical structures, to

464 select two sets of voxels; an “ROI” and a “comparator region”*. The seed voxel

465 is always contained within the ROI, and the ROI is always contained within the

466 comparator region. The comparator region is similar but not identical to the

467 set of voxels making up the major anatomical region containing the ROI. Gene

468 Finder then looks for genes which can distinguish the ROI from the comparator

469 region. Specifically, it finds genes for which the ratio (expression energy in the

470 ROI) / (expression energy in the comparator region) is high.

471 Informally, the Gene Finder first infers an ROI based on clustering the seed

472 voxel with other voxels. Then, the Gene Finder finds genes which overexpress

473 in the ROI as compared to other voxels in the major anatomical region.

474 There are three major differences between our approach and Gene Finder.

475 12

476

477 First, Gene Finder focuses on individual genes and individual ROIs in isola-

478 tion. This is great for regions which can be picked out from all other regions by a

479 single gene, but not all of them can (todo). There are at least two ways this can

480 miss out on useful genes. First, a gene might express in part of a region, but not

481 throughout the whole region, but there may be another gene which expresses

482 in the rest of the region*. Second, a gene might express in a region, but not in

483 any of its neighbors, but it might express also in other non-neighboring regions.

484 To take advantage of these types of genes, we propose to find combinations of

485 genes which, together, can identify the boundaries of all subregions within the

486 containing region.

487 Second, Gene Finder uses a pointwise metric, namely expression energy ratio,

488 to decide whether a gene is good for picking out a region. We have found better

489 results by using metrics which take into account not just single voxels, but also

490 the local geometry of neighboring voxels, such as the local gradient (todo). In

491 addition, we have found that often the absence of gene expression can be used

492 as a marker, which will not be caught by Gene Finder’s expression energy ratio

493 (todo).

494 Third, Gene Finder chooses the ROI based only on the seed voxel. This

495 often does not permit the user to query the ROI that they are interested in. For

496 example, in all of our tests of Gene Finder in cortex, the ROIs chosen tend to

497 be cortical layers, rather than cortical areas.

498 In summary, when Gene Finder picks the ROI that you want, and when this

499 ROI can be easily picked out from neighboring regions by single genes which

500 selectively overexpress in the ROI compared to the entire major anatomical re-

501 gion, Gene Finder will work. However, Gene Finder will not pick cortical areas

502 as ROIs, and even if it could, many cortical areas cannot be uniquely picked out

503 by the overexpression of any single gene. By contrast, we will target cortical

504 areas, we will explore a variety of metrics which can complement the shortcom-

505 ings of expression energy ratio, and we will use the combinatorial expression of

506 genes to pick out cortical areas even when no individual gene will do.

507 * The terms “ROI” and “comparator region” are our own; the ABI calls

508 them the “local region” and the “larger anatomical context”. The ABI uses the

509 term “specificity comparator” to mean the major anatomic region containing

510 the ROI, which is not exactly identical to the comparator region.

511 ** In this case, the union of the area of expression of the two genes would

512 suffice; one could also imagine that there could be situations in which the in-

513 tersection of multiple genes would be needed, or a combination of unions and

514 intersections.

515 Now we describe AGEA’s hierarchial clustering, and compare it to our pro-

516 posal. The goal of AGEA’s hierarchial clustering is to generate a binary tree of

517 clusters, where a cluster is a collection of voxels. AGEA begins by computing

518 the Pearson correlation between each pair of voxels. They then employ a recur-

519 sive divisive (top-down) hierarchial clustering procedure on the voxels, which

520 means that they start with all of the voxels, and then they divide them into clus-

521 ters, and then within each cluster, they divide that cluster into smaller clusters,

522 etc***. At each step, the collection of voxels is partitioned into two smaller

523 13

524

525 clusters in a way that maximizes the following quantity: average correlation

526 between all possible pairs of voxels containing one voxel from each cluster.

527 There are three major differences between our approach and AGEA’s hier-

528 archial clustering. First, AGEA’s clustering method separates cortical layers

529 before it separates cortical areas.

530 following procedure is used for the purpose of dividing a collection of voxels

531 into smaller clusters: partition the voxels into two sets, such that the following

532 quantity is maximized:

533 *** depending on which level of the tree is being created, the voxels are

534 subsampled in order to save time

535 does not allow the user to input anything other than a seed voxel; this means

536 that for each seed voxel, there is only one

537 The role of the “local region” is to serve as a region of interest for which

538 marker genes are desired; the role of the “larger anatomical context” is to be

539 the structure

540 There are two kinds of differences between AGEA and our project; differ-

541 ences that relate to the treatment of the cortex, and differences in the type of

542 generalizable methods being developed. As relates

543 indicate an ROI

544 explore simple correlation-based relationships between voxels, genes, and

545 clusters of voxels.

546 There have not yet been any studies which describe the results of applying

547 AGEA to the cerebral cortex; however, we suspect that the AGEA metrics are

548 not optimal for the task of relating genes to cortical areas. A voxel’s gene

549 expression profile depends upon both its cortical area and its cortical layer,

550 however, AGEA has no mechanism to distinguish these two. As a result, voxels

551 in the same layer but different areas are often clustered together by AGEA. As

552 part of the project, we will compare the performance of our techniques against

553 AGEA’s.

554 —

555 The Allen Brain Institute has developed interactive tools called AGEA which

556 allow an investigator to explore simple correlation-based relationships between

557 voxels, genes, and clusters of voxels. There have not yet been any studies

558 which describe the results of applying AGEA to the cerebral cortex; however,

559 we suspect that the AGEA metrics are not optimal for the task of relating

560 genes to cortical areas. A voxel’s gene expression profile depends upon both

561 its cortical area and its cortical layer, however, AGEA has no mechanism to

562 distinguish these two. As a result, voxels in the same layer but different areas

563 are often clustered together by AGEA. As part of the project, we will compare

564 the performance of our techniques against AGEA’s.

565 Another difference between our techniques and AGEA’s is that AGEA allows

566 the user to enter only a voxel location, and then to either explore the rest of

567 the brain’s relationship to that particular voxel, or explore a partitioning of

568 the brain based on pairwise voxel correlation. If the user is interested not in a

569 single voxel, but rather an entire anatomical structure, AGEA will only succeed

570 to the extent that the selected voxel is a typical representative of the structure.

571 14

572

573 As discussed in the previous paragraph, this poses problems for structures like

574 cortical areas, which (because of their division into cortical layers) do not have

575 a single “typical representative”.

576 By contrast, in our system, the user will start by selecting, not a single voxel,

577 but rather, an anatomical superstructure to be divided into pieces (for example,

578 the cerebral cortex). We expect that our methods will take into account not

579 just pairwise statistics between voxels, but also large-scale geometric features

580 (for example, the rapidity of change in gene expression as regional boundaries

581 are crossed) which optimize the discriminability of regions within the selected

582 superstructure.

583 —–

584 screen for combinations of marker genes which selectively target anatom-

585 ical structures pick delineate the boundaries between neighboring anatomical

586 structures. (b) techniques to screen for marker genes which pick out anatomical

587 structures of interest

588 , techniques which: (a) screen for marker genes , and (b) suggest new

589 anatomical maps based on

590 whose expression partitions the region of interest into its anatomical sub-

591 structures, and (b) use the natural contours of gene expression to suggest new

592 ways of dividing an organ into

593 The Allen Brain Atlas

594 –

595 to: brooksl@mail.nih.gov

596 Hi, I’m writing to confirm the applicability of a potential research project to

597 the challenge grant topic ”New computational and statistical methods for the

598 analysis of large data sets from next-generation sequencing technologies”.

599 We want to develop methods for the analysis of gene expression datasets that

600 can be used to uncover the relationships between gene expression and anatomical

601 regions. Specifically, we want to develop techniques to (a) given a set of known

602 anatomical areas, identify genetic markers for each of these areas, and (b) given

603 an anatomical structure whose substructure is unknown, suggest a map, that

604 is, a division of the space into anatomical sub-structures, that represents the

605 boundaries inherent in the gene expression data.

606 We propose to develop our techniques on the Allen Brain Atlas mouse brain

607 gene expression dataset by finding genetic markers for anatomical areas within

608 the cerebral cortex. The Allen Brain Atlas contains a registered 3-D map of

609 gene expression data with 200-micron voxel resolution which was created from

610 in situ hybridization data. The dataset contains about 4000 genes which are

611 available at this resolution across the entire cerebral cortex.

612 Despite the distinct roles of different cortical areas in both normal function-

613 ing and disease processes, there are no known marker genes for many cortical

614 areas. This project will be immediately useful for both drug discovery and clini-

615 cal research because once the markers are known, interventions can be designed

616 which selectively target specific cortical areas.

617 This techniques we develop will be useful because they will be applicable to

618 the analysis of other anatomical areas, both in terms of finding marker genes

619 15

620

621 for known areas, and in terms of suggesting new anatomical subdivisions that

622 are based upon the gene expression data.

623 —-

624 It is likely that our study, by showing which areal divisions naturally fol-

625 low from gene expression data, as opposed to traditional histological data, will

626 contribute to the creation of

627 there are clear genetic or chemical markers known for only a few cortical

628 areas. This makes it difficult to target drugs to specific

629 As part of aims (1) and (5), we will discover sets of marker genes that pick

630 out specific cortical areas. This will allow the development of drugs and other

631 interventions that selectively target individual cortical areas. As part of aims

632 (2) and (5), we will also discover small panels of marker genes that can be used

633 to delineate most of the cortical areal map.

634 With aims (2) and (4), we

635 There are five principals

636 In addition to validating the usefulness of the algorithms, the application of

637 these methods to cerebral cortex will produce immediate benefits that are only

638 one step removed from clinical application.

639 todo: remember to check gensat, etc for validation (mention bias/variance)

640 Why it is useful to apply these methods to cortex

641 There is still room for debate as to exactly how the cortex should be parcellated

642 into areas.

643 The best way to divide up rodent cortex into areas has not been completely

644 determined,

645 not yet been accounted for in

646 that the expression of some genes will contain novel spatial patterns which

647 are not account

648 that a genoarchitectonic map

649 This principle is only applicable to aim 1 (marker genes). For aim 2 (partition

650 a structure in into anatomical subregions), we plan to work with many genes at

651 once.

652 tood: aim 2 b+s?

653 Principle 5: Interoperate with existing tools

654 In order for our software to be as useful as possible for our users, it will be

655 able to import and export data to standard formats so that users can use our

656 software in tandem with other software tools created by other teams. We will

657 support the following formats: NIFTI (Neuroimaging Informatics Technology

658 Initiative), SEV (Allen Brain Institute Smoothed Energy Volume), and MAT-

659 LAB. This ensures that our users will not have to exclusively rely on our tools

660 when analyzing data. For example, users will be able to use the data visualiza-

661 tion and analysis capabilities of MATLAB and Caret alongside our software.

662 16

663

664 To our knowledge, there is no currently available software to convert between

665 these formats, so we will also provide a format conversion tool. This may be

666 useful even for groups that don’t use any of our other software.

667 todo: is “marker gene” even a phrase that we should use at all?

668 note for aim 1 apps: combo of genes is for voxel, not within any single cell

669 , as when genetic markers allow the development of selective interventions;

670 the reason that one can be confident that the intervention is selective is that it

671 is only turned on when a certain combination of genes is turned on and off. The

672 result procedure is what assures us that when that combination is present, the

673 local tissue is probably part of a certain subregion.

674 The basic idea is that we want to find a procedure by

675 The task of finding genes that mark anatomical areas can be phrased in

676 terms of what the field of machine learning calls a “supervised learning” task.

677 The goal of this task is to learn a function (the “classifier”) which

678 If a person knows a combination of genes that mark an area, that implies

679 that the person can be told how strong those genes express in any voxel, and

680 the person can use this information to determine how

681 finding how to infer the areal identity of a voxel if given the gene expression

682 profile of that voxel.

683 For each voxel in the cortex, we want to start with data about the gene

684 expression

685 There are various ways to look for marker genes. We will define some terms,

686 and along the way we will describe a few design choices encountered in the

687 process of creating a marker gene finding method, and then we will present four

688 principles that describe which options we have chosen.

689 In developing a procedure for finding marker genes, we are developing a

690 procedure that takes a dataset of experimental observations and produces a

691 result. One can think of the result as merely a list of genes, but really the result

692 is an understanding of a predictive relationship between, on the one hand, the

693 expression levels of genes, and, on the other hand, anatomical subregions.

694 One way to more formally define this understanding is to look at it as a

695 procedure. In this view, the result of the learning procedure is itself a procedure.

696 The result procedure provides a way to use the gene expression profiles of voxels

697 in a tissue sample in order to determine where the subregions are.

698 This result procedure can be used directly, as when an experimenter has

699 a tissue sample and needs to know what subregions are present in it, and,

700 if multiple subregions are present, where they each are. Or it can be used

701 indirectly; imagine that the result procedure tells us that whenever a certain

702 combination of genes are expressed, the local tissue is probably part of a certain

703 subregion. This means that we can then confidentally develop an intervention

704 which is triggered only when that combination of genes are expressed; and to

705 the extent that the result procedure is reliable, we know that the intervention

706 will only be triggered in the target subregion.

707 We said that the result procedure provides “a way to use the gene expression

708 profiles of voxels in a tissue sample” in order to “determine where the subregions

709 are”.

710 17

711

712 Does the result procedure get as input all of the gene expression profiles

713 of each voxel in the entire tissue sample, and produce as output all of the

714 subregional boundaries all at once?

715 it is helpful for the classifier to look at the global “shape” of gene expression

716 patterns over the whole structure, rather than just nearby voxels.

717 there is some small bit of additional information that can be gleaned from

718 knowing the

719 Design choices for a supervised learning procedure

720 After all,

721 there is a small correlation between the gene expression levels from distant

722 voxels and

723 Depending on how we intend to use the classifier, we may want to design it

724 so that

725 It is possible for many things to

726 The choice of which data is made part of an instance

727 what we seek is a procedure

728 partition the tissue sample into subregions.

729 each part of the anatomical structure

730 must be One way to rephrase this task is to say that, instead of searching

731 for the location of the subregions, we are looking to partition the tissue sample

732 into subregions.

733 There are various ways to look for marker genes. We will define some terms,

734 and along the way we will describe a few design choices encountered in the

735 process of creating a marker gene finding method, and then we will present four

736 principles that describe which options we have chosen.

737 In developing a procedure for finding marker genes, we are developing a

738 procedure that takes a dataset of experimental observations and produces a

739 result. One can think of the result as merely a list of genes, but really the result

740 is an understanding of a predictive relationship between, on the one hand, the

741 expression levels of genes, and, on the other hand, anatomical subregions.

742 One way to more formally define this understanding is to look at it as a

743 procedure. In this view, the result of the learning procedure is itself a procedure.

744 The result procedure provides a way to use the gene expression profiles of voxels

745 in a tissue sample in order to determine where the subregions are.

746 This result procedure can be used directly, as when an experimenter has

747 a tissue sample and needs to know what subregions are present in it, and,

748 if multiple subregions are present, where they each are. Or it can be used

749 indirectly; imagine that the result procedure tells us that whenever a certain

750 combination of genes are expressed, the local tissue is probably part of a certain

751 subregion. This means that we can then confidentally develop an intervention

752 which is triggered only when that combination of genes are expressed; and to

753 the extent that the result procedure is reliable, we know that the intervention

754 will only be triggered in the target subregion.

755 18

756

757 We said that the result procedure provides “a way to use the gene expression

758 profiles of voxels in a tissue sample” in order to “determine where the subregions

759 are”.

760 Does the result procedure get as input all of the gene expression profiles

761 of each voxel in the entire tissue sample, and produce as output all of the

762 subregional boundaries all at once?

763 Or are we given one voxel at a time,

764 In the jargon of the field of machine learning, the result procedure is called

765 a classifier.

766 The task of finding genes that mark anatomical areas can be phrased in

767 terms of what the field of machine learning calls a “supervised learning” task.

768 The goal of this task is to learn a function (the “classifier”) which

769 If a person knows a combination of genes that mark an area, that implies

770 that the person can be told how strong those genes express in any voxel, and

771 the person can use this information to determine how

772 finding how to infer the areal identity of a voxel if given the gene expression

773 profile of that voxel.

774 For each voxel in the cortex, we want to start with data about the gene

775 expression

776 single voxels, but rather groups of voxels, such that the groups can be placed

777 in some 2-D space. We will call such instances “pixels”.

778 We have been speaking as if instances necessarily correspond to single voxels.

779 But it is possible for instances to be groupings of many voxels, in which case

780 each grouping must be assigned the same label (that is, each voxel grouping

781 must stay inside a single anatomical subregion).

782 In some but not all cases, the groups are either rows or columns of voxels.

783 This is the case with the cerebral cortex, in which one may assume that columns

784 of voxels which run perpendicular to the cortical surface all share the same areal

785 identity. In the cortex, we call such an instance a “surface pixel”, because such

786 an instance represents the data associated with all voxels underneath a specific

787 patch of the cortical surface.

788 19

789

790