Brain Gene Expression Analysis: a MATLAB toolbox for the analysis of brain-wide gene-expression data

The Allen Brain Atlas project (ABA) generated a genome-scale collection of gene-expression profiles using in-situ hybridization. These profiles were co-registered to the three-dimensional Allen Reference Atlas (ARA) of the adult mouse brain. A set of more than 4,000 such volumetric data are available for the full brain, at a resolution of 200 microns. These data are presented in a voxel-by-gene matrix. The ARA comes with several systems of annotation, hierarchical (40 cortical regions, 209 sub-cortical regions in the whole brain), or non-hierarchical (12 regions in the left hemisphere, with refinement into 94 regions, and cortical layers). The high-dimensional nature of this dataset and the possible connection between anatomy and gene expression pose challenges to data analysis. We developed the Brain Gene Expression Analysis Toolbox, whose functionalities include: determination of marker genes for brain regions, statistical analysis of brain-wide co-expression patterns, and the computation of brain-wide correlation maps with cell-type specific microarray data.


Chapter 2
The voxel-based Allen Atlas of the adult mouse brain 2.1 Presentation of the dataset 2. 1

.1 Co-registered in situ hybridization data
The adult mouse brain is partitioned into V = 49, 742 cubic voxels of side 200 microns, to which in situ hybridization data are registered [2,3] for thousands of genes. For computational purposes, these gene-expression data can be arranged into a voxel-by-gene matrix.
For each voxel v, the expression energy of the gene g is a weighted sum of the greyscale-value intensities I evaluated at the pixels p intersecting the voxel: where M (p) is a Boolean mask that equals 1 if the gene is expressed at pixel p and 0 if it is not.

Data matrices and gene filters Coronal atlas
Some genes in the Allen atlas of the adult mouse brain gave rise to an ISH experiment and coronal sectioning of an entire brain. The resulting data constitute the coronal atlas.
The coronal atlas contains brain-wide data for G all = 4, 104 genes (this is a subset of the genes for which ISH sagittal sectioning took place for a hemisphere after ISH, see next section for documentation on the sagittal atlas). The corresponding voxel-by-gene data matrix has size V = 49, 742 by G all , and is contained in the file ExpEnergy.mat. The list of genes arranged in the same order as the columns of the data matrix are obtained by using the function get genes.m. 1 load ( ' ExpEnergy . mat ' ); 2% the g−th column o f voxel−by−gene matrix E corresponds 3% to the gene geneNamesAll ( g ) 4 geneNamesAll = get_genes ( Ref , ' allNoDup ' , ' allen ' ); 5% Entrez Ids are arranged i n the same order 6% as gene names ( unresolved Entrez i d s are treated as zero ) 7 geneEntrezIdsAll = get_genes ( Ref , ' allNoDup ' , ' entrez ' ); The matrix in ExpEnergyTop75Percent.mat consists of the 3,041 columns of the matrix defined in Equation 2.1 that are best correlated with the corresponding genes in the sagittal atlas. The names and Entrez ids of the genes 1 are obtained as follows: 1 load ( ' E x p E n e r g y t o p 7 5 P e r c e n t . mat ' ); 2% the g−th column o f voxel−by−gene matrix E 3% corresponds to the gene genesAllen ( g ) 4 genesAllen = get_genes ( Ref . Coronal , ' top75corrNoDup ' , ' allen ' ); 5% Entrez Ids are arranged i n the same order It should be noted that the Entrez ids that are not resolved are represented by zeroes in geneEntrezIds, so Entrez ids should be resolved in fine rather than used during computations.
The start-up file mouse_start_up.m loads the Allen Reference Atlas stored in the file refAtlas.mat and the data matrix ExpEnergytop75percent.mat, which can be used through the structure Ref. • Example 1. Systems of brain annotation. The variable cor=Ref.Coronal contains the coronal atlas, and it can be used to check some of the contents of Table 2. 1. a detailed description of the fields in the structure Ref.Coronal.Annotations can be found in the next section. 1 The differences between the two code snippets are the data matrix loaded, and the filter used as teh second argument of the function get_genes. One can check that after executing any of the snippets, the number of gene names equal the number of columns of the data matrix.

From matrices to volumes 2.2.1 The Allen Reference Atlas (ARA)
The ARA comes in six different versions, described in Table 2. 1. Each of these versions corresponds to an annotation of the three-dimensional grid by digital ids. Each of these digital ids corresponds to a brain region. The correspondence between ids and names of brain regions can be resolved using the fields Ref.Coronal.ids and Ref.Coronal.labels.
• Example 2. From brain regions to digital ids in the ARA. Consider the fine annotation (identifier index 4 in Table 1) and let us work out a threedimensional grid with ones at the voxels corresponding to the caudoputamen in this annotation, and zeroes everywhere else. We can use this grid to compute the number of voxels in the caudoputamen. 1% the three−dimensional g r i d containing the f i n e annotation o f the brain 2 annotFine = get_annotation ( cor , ' fine ' ); 3% the names o f brain r e g i o n s 4 labels = ann . labels { 4 }; 5% the numerical i d s o f the r e g i o n s 6 ids = ann . ids { 4 }; 7% where i s the caudoputamen i n the l i s t ? There may be spaces 8% i n the l a b e l s , s k i p them 9 labelsNoSpace = regexprep ( labels , '\ W ' , ' ' ); 10 ca ud ou pu ta me nI nd ex = find ( strcmp ( labelsNoSpace , ' Caudoputamen ' ) == 1 ); 11 caudouputamenId = ids ( c au dou pu ta me nI nd ex ); 12% put z e r o s at a l l the voxels that do not belong to caudoputamen 13 volCaudoputamen = annotFine ; 14 volCaudoputamen ( volCaudoputamen~= caudouputamenId ) = 0; 15 volCaudoputamen ( volCaudoputamen == caudouputamenId ) = 1; 16% count the voxels i n caudoputamen 17 nu m V o x I n C a u d o p u t a m e n = sum ( sum ( sum ( volCaudoputamen ) ) ); The value of the variable numVoxInCaudoputamen computed in the above code snippet should be 1248.

From gene-expression vectors to volumes
At a resolution of 200 microns, the Allen Reference Atlas (ARA) is embedded in a three-dimensional grid of size 67 × 41 × 58. Out of the V tot = 67 × 41 × 58 in the grid, V = 49, 742 are in in the brain according to the ARA. The brain voxels can be mapped to a three-dimensional grid using the function make_volume_from_labels.m and a specified voxel filter, corresponding to one of the versions of the ARA (as per Table 2.1).
• Example 3. Whole-brain filter. Let us take the whole first column of the data matrix and map it to a three-dimensional grid, using the whole-brain filter, corresponding to the 'standard' annotation: 1 wholeBrainFilter = Ref . Coronal . Annotations . Filter { 1 }; 2 display ( wholeBrainFilter ); 3 brainFilter = get_voxel_filter ( cor , wholeBrainFilter ); 4% a column vector with 49 ,742 elements 5 col1 = E ( : , 1 ); 6 display ( size ( col1 ) ) 7% map t h i s column vector to a volume 8 vol1 = m a k e _ v o l u m e _ f r o m _ l a b e l s ( col1 , brainFilter ); 9 display ( size ( vol1 ) ); Some of the annotations do not extend to the whole brain, as can be seen from Table  2. 1. when comparing columns of the matrix of gene-expression energies to regions in the ARA, it is important to restrict the matrix to the rows that correspond to annotated voxels. For each system of annotation, the field Ref.Coronal.Annotations.Filter is the list of voxels in a 67 × 41 × 58 grid (not a list of row indices in the matrix of expression energies), that are annotated. The example below shows how to recover the list of row indices in the gene-expression matrix corresponding to a given filter.
• Example 4. More filters. For instance, one can check that the following two snippets produce the same matrix EFiltered, consisting of the rows of the full geneexpression matrix corresponding to voxels in the fine annotation: 1 % Compute the s e t o f rows using a volume o f 2 %i n d i c e s and the f i l t e r 3 cor = Ref . Coronal ; 4 identifierIndex = 4; 5 brainFilter = get_voxel_filter ( cor , ' brainVox ' ); 6 numVox = numel ( brainFilter ); 7 % l a b e l the voxels by i n t e g e r s 8 indsBrainVoxels = 1 : numVox ; 9 % arrange the i n t e g e r s i n a volume 10 indsVol = m a k e _ v o l u m e _ f r o m _ l a b e l s ( indsBrainVoxels , brainFilter ); 11 filter = get_voxel_filter ( cor , ann . filter { identifierIndex } ); 12 % r e s t r i c t the volume to the voxels that are i n the f i l t e r 13 indsFiltered = indsVol ( filter ); 14 EFiltered = E ( indsFiltered , : );}} 1 %apply the f i l t e r to each column o f the data matrix 2 cor = Ref . Coronal ; 3 identifierIndex = 4; 4 brainFilter = get_voxel_filter ( cor , ' brainVox ' ); 5 filter = get_voxel_filter ( cor , ann . filter { identifierIndex } ); 6 numGenes = size ( E , 2 ); 7 % r e s t r i c t each colum o f the data matrix to voxels 8 % that are i n the f i l t e r 9 for gg = 1 : numGenes Note also that the voxel filter can be recomputed from the three-dimensional annotation and the brain-wide filter: 1% the three−dimensional g r i d containing the f i n e annotation o f the brain 2 annotFine = get_annotation ( cor , ' fine ' ); 3 labels = ann . labels { 4 }; 4 ids = ann . ids { 4 }; 5 labelsNoSpace = regexprep ( labels , '\ W ' , ' ' ); 6 ca ud ou pu ta me nI nd ex = find ( strcmp ( labelsNoSpace , ' Caudoputamen ' ) == 1 ); 7 caudouputamenId = ids ( c au dou pu ta me nI nd ex ); 8 volCaudoputamen = annotFine ; 9 volCaudoputamen ( volCaudoputamen~= caudouputamenId ) = 0; 10 volCaudoputamen ( volCaudoputamen == caudouputamenId ) = 1; 11 p l o t _ i n t e n s i t y _ p r o j e c t i o n s ( volCaudoputamen ); The fineAnatomy filter. It is manifest from the projection of the characteristic function of the caudoputamen in the fine annotation ( Figure 2.2) that the fine annotation contains only the left hemisphere of the brain. We can confirm this by applying the voxel filter of the fine annotation (called fineAnatomy, which is the value of Ref.Coronal.Annotations.filter{4}) to the gene-expression vector of Gabra6 . The following snippet should reproduce Figure 2.3, which is equivalent to Figure 5.1, with all the voxels that are not in the fineAnatomy filter filled with zeros.  1% take the brain−wide expression o f Gabra6 2 indexGabra6 = find ( strcmp ( genesAllen , ' Gabra6 ' ) == 1 ); 3 columnGabra6 = E ( : , indexGabra6 ); 4 volGabra6 = m a k e _ v o l u m e _ f r o m _ l a b e l s ( columnGabra6 , brainFilter ); 5% c o n s i d e r the voxel f i l t e r o f the ' f i n e ' annotation 6 fineFilter = get_voxel_filter ( cor , ' fineAnatomy ' ); 7% take the expression values at the voxels i n the f i l t e r , and arrange them 8% i n a column vector ( one can check that 9% dataGabra6Filtered has the same number o f 10% elements as f i n e F i l t e r 11 colG abra6F iltere d = volGabra6 ( fineFilter ); 12% map these values to a three−dimensional grid , using f i n e F i l t e r  • Example 8. Brain-wide versus left hemisphere. As can be seen on Table  2.1, the standard annotation contains all the voxels in the brain (49,742 of them). We can compute the characteristic function of the caudoputamen in this annotation, and check  7% where i s the caudoputamen i n the l i s t ? 8 labelsNoSpace = regexprep ( labels , '\ W ' , ' ' ) ; 9 ca ud ou pu ta me nI nd ex = find ( strcmp ( labelsNoSpace , ' Caudoputamen ' ) == 1 ); 10 caudouputamenId = ids ( c au dou pu ta me nI nd ex ); 11% put z e r o s at a l l the voxels that do not belong to caudoputamen 12 volCaudoputamen = annotStandard ; 13 volCaudoputamen ( volCaudoputamen~= caudouputamenId ) = 0; 14 volCaudoputamen ( volCaudoputamen == caudouputamenId ) = 1; 15% count the voxels i n caudoputamen 16 n u m V o x I n C a u d o p u t a m e n S t a n d a r d = sum ( sum ( sum ( volCaudoputamen ) ) ); 17 display ( n u m V o x I n C a u d o p u t a m e n S t a n d a r d ); 18% reproduce Figure

Sections
The function flip_through_sections.m allows to go through the sections of a (67×41×58) volume, of a kind specified in the options. It pauses between sections. The duration of the pause is one second by default, it can be adjusted using the field secondsOfPause of the options. If the value of secondsOfPause is negative, the user will have to press a key to display the next section. • Example 9. Sections of the average across all genes in the dataset. The following code allows to visualize the coronal, sagittal and axial sections of the average expression E across all genes in the data matrix, defined in Equation 2.2:

Relation between the various annotations in the ARA
The field Ref.Coronal.Annotations.parentIdx in the reference data structure gives the indices (not the numerical id) of the parents of the regions in a given annotation, arranged in the same order as these regions. If a region does not have a parent in the considered system of annotation, the corresponding entry in parentIdx is zero. NB: the list of regions in the 'cortex' annotation is not closed under the operation of taking the parent of a region in a hierarchy. The field parentSymbols has to be used instead of parentIdx to work out the hierarchy.
• Example 10. Hierarchical and non-hierarchical annotations. Let us check that the big12 and fine annotations are non-hierarchical, and investigate parents and descendants of the caudoputamen and of the cerebellum in the standard annotation.
Let us work out which regions (if any) correspond to an empty set of voxels in the standard annotation. Only two of the 7 regions corresponding to an empty set of voxels are therefore leaves of the hierarchical tree. We can conclude that they are too small to be represented by voxels at a spatial resolution of 200 microns. The other 5 regions have descendants in the hierarchy, but they can only be resolved by taking the reunion of the voxels belonging to their descendants (and descendants thereof). One can note a curiousity about Striatum dorsal region, which only has one descendant, caudoputamen (which corresponds to the index uu=1 in the above loop). The inclusion of caudoputamen in the dorsal region of the striatum is therefore trivial, and the two labels Striatum dorsal region and Caudoputamen can be treated as synonyms of each other.
The fine annotation compared to the big12 annotation. Let us show that the 'fine annotation is a refinement of the big12 annotation. It is manifest from Table 2.1 that the fine annotation covers fewer voxels than, but we can check that 1) all these voxels are also in the big12 annotation, and 2) each region in the fine annotation intersects only one region in the big12 annotation. The toolbox contains these code snippets as two functions that work out the refinement of the big12 annotation, and the organisation of the fine annotation into larger regions of the brain. One can check that the following reproduces the above results: For instance, the 16-th region in the fine is Nucleus Accumbens. Check it: ( 16 ) ); 2% check that t h i s region i s a the r i g h t index i n annotationFineToBig12 3 display ( a n n o t a t i o n F i n e T o B i g 1 2 . labelsFine ( nAccIndex ) ); 4% In which region o f Big12 i s i t included ? 5 display ( a n n o t a t i o n F i n e T o B i g 1 2 . l a b e l s P a r e n t B i g A t l a s ( nAccIndex ) ); 6% Apart from nucleus accumbens which are the other subregions o f the Genes versus neuroanatomy 3.1 Localization scores 3. 1

.1 Localization scores of a single gene in the ARA
Let us define [4] the localization score λ ω (g) of a gene g in a region ω as the ratio of the squared L 2 -norm of the expression energy of gene g in region ω to the squared L 2 -norm of the expression energy of gene g in the set Ω of voxels that are annotated in the version of the ARA that contains ω.
It can be computed from a voxel-by-gene matrix (with one score per column for a fixed region ω), using the function localization_from_id.m, given the numerical id of the region ω and the numerical identifier corresponding to an annotation containing ω.
• Example 13. One gene, one region. Let us compute the localization score of Pak7 in the cerebral cortex, as defined by the big12 annotation of the left hemisphere.

Localization scores of sets of genes in the ARA
where G is the number of genes in our dataset (G = 3, 041 by default when using the start-up file mouse_start_up.m).
Let us define the localization score in the brain region ω of a weighted set of genes encoded by Equation 3.2 as • Example 16. Let us compute the best generalized localization scores in the big12 annotation.    The fitting score φ ω (g) of a gene g in a region ω is defined where E norm g is the L 2 -normalized g-th column E g of the matrix of gene-expression energies: the symbol Ω denotes the set of voxels in a given system of annotation that contains region ω. The definition is the only decreasing affine function of the squared L 2 norm of the difference between the normalized gene-expression vector of gene g and the characteristic function χ ω of the region ω: where the denominator in Equation 3.6 ensure the L 2 -normalization V v=1 χ ω (v) 2 = 1.
Given a gene, a system of annotation chosen among the ones in Table 2.1, and the numerical id of a region in this system of annotation, the function fitting_from_id.m computes the fitting score of the gene to the corresponding region.
• Example 17. One gene, one region in a given annotation. Let us compute the fitting score of Pak7 in ω = Cerebral cortex, in the big12 annotation (hence identifierIndex = 5, see Table 2 The function fitting_from_id.m can also used to compute the fitting score of several genes. If the second argument of the function is a voxel-by-gene matrix with p columns, the function returns an array of p fitting scores arranged in the same order as the columns. Given a version of the ARA (specified by the index identifierIndex), a region-bygene matrix of all the fitting scores of all genes corresponding to the columns of the data matrix. Note that the columns of this region-by-gene score matrix do not sum to a constant (the squares of the entries of each column sum to the square of the fraction of the gene-expression that projects onto the set of voxels in the annotation, which is at most 1).
• Example 18. Fitting scores of all genes in all the regions in the big12 annotation.

Fitting scores of sets of genes in the ARA
Like the localization score, the fitting score can be generalized to linear combinations of sets of genes: where E norm α is the L 2 -normalized gene-expression vector corresponding to the coefficients (α 1 , . . . , α G ) which can be minimized wrt the weights of the genes using Matlab code implementing an interior-point method by Koh [5]:             Genes versus genes: co-expression networks 4

.1 Co-expression networks of genes in the Allen Atlas
The co-expression of two genes in the Allen Atlas is defined as the cosine similarity between their gene-expression vectors in voxel space. Given a voxel-by-gene matrix containing the brain-wide expression energies (as in Equation 2.1), the corresponding gene-by-gene matrix of co-expressions of the full set of genes, or coExpr full , is the symmetric matrix with entries equal to the co-expression of pairs of genes, as in Equation 4. 1. (4.1) The matrix coExpr full can be computed as follows in Matlab: The function co_expression_matrix.m takes a voxel-by-gene matrix as an argument and returns the gene-by-gene co-expression matrix defined by Equation 4.1, which equals the matrix coExpressionFull defined in the above code snippet if the full voxel-by-gene matrix of gene expression energies is used as an argument. Other versions of the Allen Atlas than the brain-wide standard annotations can be specified in the options to restrict the voxels to one of the annotations described in Table 2. 1.
• Example 20. Distribution of co-expression coefficients. The diagonal elements of the co-expression matrix equal 1 by construction, and the co-expression matrix is symmetric. Hence, the distibution of co-expression coefficients in the atlas is given by the upper diagonal coeefficients of the co-expression matrix, which can be extracted using the function upper_diagonal_coeffs.m, as in the code snipped below, which plots the distribution of brain-wide co-expression coefficients (Figure 4.1)) and compares it to the one of co-expression coefficients in the left hemisphere (as defined in the big12 annotation, Figure 4 4.2 Monte Carlo analysis of brain-wide co-expression networks 4

.2.1 Special sets of genes versus full atlas
The gene-by-gene matrix coExpr full defines a universe in which we would like to study co-expression networks of special sets of genes, in a probabilistic way.
Given a set of G special genes of interest, corresponding to the column indices (g 1 , . . . , g special ) in the data matrix, their co-expression matrix coExpr special is obtained by extracting the submatrix of coExpr full corresponding to these indices (see Equation 4.2.4).

Cumulative distribution function of co-expression coefficients in sets of genes drawn from the Allen Brain Atlas
Having observed (Figure 4.1) that the distribution of pairwise co-expression coefficients of genes in the whole coronal atlas is roughly linear in a large domain of co-expression, we can study the cumulative distribution function of the co-expression coefficients in the special set, and compare it to the one resulting from random sets of genes (with the same number of genes as the special set, in order to eliminate the sample-size bias).
These cumulative distribution functions are evaluated in the following way. Again let G special denote the size of the matrix coExpr special , i.e. the number of genes from which coExpr special was computed. Consider the set of entries above the diagonal of coExpr special above the diagonal (which are the meaningful quantities in coExpr special ): The elements of this set are numbers between 0 and 1. For every number between 0 and 1, the cumulative distribution function (c.d.f.) of C special , denoted by cdf C is defined as the fraction of the elements of C special that are smaller than this number: where C special = G special (G special − 1)/2.
For any set of genes, cdf special is a growing function cdf special (0) = 0 and cdf special (1) = 1. For highly co-expressed genes, the growth of cdf special is concentrated at high values of the argument (in the limit where all the genes in the special set have the same brain-wide expression vector, all the entries of the co-expression matrix go to 1 and the cumulative distribution function converges to a Dirac measure supported at 1). To compare the function cdf special to what could be expected by chance, let us draw R random sets of G special genes from the Atlas, compute their co-expression network by extracting the corresponding entries from the full co-expression matrix of the atlas (coExpr full ). This induces a family of R growing functions cdf i , 1 ≤ ı ≤ R on the interval [0, 1] From this family of functions, we can estimate a mean cumulative distribution function cdf of the co-expression of sets of G special genes drawn from the Allen Atlas, by taking the mean of the values of cdf i across the random draws: The functions cdf special , cdf and cdf dev are fields of the output of the function cumul_co_expr.m whose usage is illustrated in the example below.
Example 21. Cumulative distribution function of co-expression of a special set of genes. Consider the set of 288 genes from the NicSNP database, whose position in the data matrix is encoded in the Matlab file nicotineGenes.mat. NB: the code snippet below can be applied to any special set of genes upon changing the variable indsSmall. It should reproduce

Thresholding the co-expression matrix
The co-expression matrix coExpr special corresponding to a special subset of the genes in the Allen Atlas (Equation 4. 2.4) is symmetric, like coExpr full , and its entries are in the interval [0, 1]. It can be mapped to a weighted graph in the following way (see [6] for details). The vertices of the graph are the genes, and the edges are as follows: -genes g and g are linked by an edge if their co-expression is strictly positive.
-If an edge exists, it has weight coExpr special gg .
Let us define the following thresholding procedure on co-expression graphs: given a threshold ρ between 0 and 1, put to zero all the entries of coExpr special that are lower than this coefficient. The underlying graph is obtained by taking the graph corresponding to coExpr special , and cutting all the links with weight below ρ.
The more-co-expressed a set of genes is, the larger the connected components of the thresholded graphs uderlying coExpr special ρ will be, for any value of the threshold ρ. For instance we can study the average size of connected components of thresholded co-expression matrices and the size of the largest connected component as a function of the threshold ρ : 12) where N ρ (k) is the number of connected components with size k. The connected components are worked out using the implementation of Tarjan's algorithm [17] in the Matlab function graphconncomp.m.

Statistics of sizes of connected components
For any quantity worked out from the special co-expression matrix co-expression matrix coExpr special defined in Equation , we can simulate its probability distribution by repeatedly drawing random random sets of genes from the atlas, and recomputing the same quantity on for this set.
Let us use the above-defined thresholding procedure to study a set of G special = 288 genes obtained by intersecting the NicSNP database [16] with the set of G = 3, 041 genes given by get_genes( Ref.Coronal, 'top75CorrNoDup', 'allen'). We would like to acertain whether this set of genes is more co-expressed than expected by chance for a set of this size taken from genesAllen. At each value of a regular grid the threshold ρ between zero and 1, the function co_expression_island_bootstrap.m computes the maximal size and average size of connected components of the thresholded co-expression graph, and draws R random sets of genes of size G special from the atlas. This induces a distribution of R partitions of sets of G special genes into connected components, obtained by applying Tarjan's algorithm to each of the R sets of genes.
Example 22. Consider the set of 288 genes from the NicSNP database, whose position in the data matrix is encoded in the Matlab file nicotineGenes.mat.

Neuroanatomical properties of connected components Relation between fitting scores and co-expression
If the co-expression of pairs of genes is defined as the cosine of the angle between their expression vectors in voxel space, as in Equation 4.1, the fitting score of a gene g to a region ω of the brain equals the co-expression of gene g and a (hypothetical) gene whose expression profile would be proportional to the characteristic function of region ω: where the second equality comes from the normalization of E mathrmnorm g and χ ω in voxel space.

Fitting scores of sums of gene-expression vectors
At a given level of the threshold on co-expression, a set of genes is partitioned into connected components induced by the graph underlying the thresholed matrix defined in Equation .
For a set of genes of fixed size numGenes, the function fitting_distribution_in_atlas.m estimates the distribution of fitting scores of the sum of sets of genes of a given size, extracted from the Allen Atlas.
In these functions, the genes are not weighted by coefficients to be optimized, they are simply summed over, so that the set of gene-expression vectors is different from the one explored in marker genes. The fitting score φ sum ω of the sum of genes depends only on a list of K distinct K genes {g 1 , . . . , g K }, and a region ω in the Allen Reference Atlas: where E norm ({g 1 , . . . , g K }) is the sum of gene-expression vectors of the genes {g 1 , . . . , g K }, normalized in the L 2 sense: .  One can search the results by P -value of fitting scores, and/or size of connected components, and/or brain region: Example 25. Neuroanatomy of nicotine-related genes in the big12 annotation (continued). 1% show a l l the connected components , at any value o f 2% the threshold , whose f i t t i n g s c o r e to any region i s estimated to be i n Gene-based and cell-based expression data The computational techniques exposed in this chapter allow to reproduce the results of [7,11,12,8,9,10] (see also [13,14,15] for analyses of the Allen Atlas data in terms of cell types).

Cell-type-specific microarray data
The file G_t_means.txt contains a matrix of microarray data. The rows correspond to the cell-type-specific samples coming from several studies and analyzed in [18]. The columns correspond to genes (the matrix is a type-by-gene matrix). The following code snippet (included in the file cell_type_start_up.m) defines two matrices with the same numbers of columns, corresponding to the intersection of the sets of genes in the coronal Allen Brain Atlas and in the microarray data:

Brain-wide correlation profiles
The function cell_types_correls.m computes the voxel-by-type correlation matrix between the coronal Allen Brain Atlas and a set of cell types, as in equation 5. 1. 2) Example 26. Brain-wide correlation profiles between cell-type-specific data and the Allen Atlas. To compute the correlations for all the available cell types, the type indices specified in the variables span the whole range of row indices of the matrix C:  The results are saved in the file cellTypesCorrelations.mat. It has V = 49, 742 rows and T = 64 columns (the matrix cellTypesCorrelations is a voxel-by-type matrix, consistently with the definition of Corr(v, t) in Equation 5.1), so that each column can be mapped to a volume and visualized in exactly the same way as a column of the data matrix E: Example 27. Visualization of brain-wide correlation profiles. Let us plot maximal-intensity projections of the correlations between the 20-th row of C and all the rows of the Allen Atlas.

Brain-wide density estimates of cell types
To decompose the gene-expression at every voxel in the Allen Atlas into its cell-typespecific components, let us introduce the positive quantity ρ t (v) denoting the contribution of cell-type t at voxel v, and propose the following linear model: Both sides are estimators of the number of mRNAs for gene g at voxel v. The residual term in Equation 5.4 reflects the fact that T = 64 cell types are not enough to sample the whole diversity of cell types in the mouse brain, as well as noise in the measurements, reproducibility issues, the non-linearity of the relations between numbers of mRNAs, expression energies and microarray data.
To find the best fit of the model, we have to minimize the residual term by solving the following problem : (ρ t (v)) 1≤t≤T,1≤v≤V = argmin φ∈R + (T,V ) E E,C (φ), (5.5) where (as defined in Equation 5.1), one can rank the regions in a given annotation of the Allen Reference Atlas by computing the average correlation across in each of the regions: For each cell type t, each region r in the ARA is ranked according to its average correlation. Call its rank χ(r). In particular, the correlation profile yields a top region ρ χ (t) in the ARA, for which χ(ρ χ (t)) = 1. This is the region whose voxels are most highly correlated to cell type t on average.
Ranking of regions by estimated density of a cell type. Consider a system of annotation of the ARA, taken from Table 2.1 Big12 is the default option in the Matlab implementation). For each cell-type index t and each region V r (where r is an integer index taking values in [1 ..R], where R is the number of regions in the chosen system of annotation), we can compute the contribution of the voxels of the region to the brainwide density profile of cell type t: Given a cell type t, the regions in the ARA is ranked according to its average correlation. Call its rank φ(r). In particular, the correlation profile yields a top region ρ φ (t) in the ARA, for which φ(ρ φ (t)) = 1. This is the region whose voxels bring the largest total contribution to the vector ρ t .
These two rankings are implemented in the function classify_pattern.m (see example below). The values of the fraction of total density and average correlation for a given cell type can be plotted (in the big12 version of the ARA), using 1 the functions plot_correlations_for_big12.m and plot_densities_for_big12.m. Given a region in the ARA (our best guess according to one of the above criteria), we need to pick a section that intersects this region. The function cell_type_vol_prepare.m implements the choice of the section. It chooses the section (which can be sagittal, coronal or axial, specified by the field options.sectionStyle of the options) that intersect the desired region (specified by the fields options.identifierIndex and options.regionIndexForSection) along the largest number of voxels, unless the field options.customIndex equals 1. Then it works out the section (of the required style) at a position specified by option.desiredIndex.

Chapter 6
Clustering 6.1 Kullback-Leibler distance from gene-expression to a given probability distribution One can use the Kullback-Leibler (KL) divergence to compare the brain-wide expression of a gene to a given probability distribution over the brain.
This probability distribution can be for instance the normalized average E avg full across all genes in the coronal atlas for instance KL avg (g) : whereẼ avg full andẼ g are probability densities over the brain obtained by normalization: . (6.2)

Construction of a bipartite graph from a voxelby-gene matrix
The Allen Atlas can be mapped to a weighted bipartite graph in the following way: • the first set of vertices consists of voxels, numbered from 1 to V , • the second set of vertices consists of genes, numbered from 1 to G, • each of the edges connects one voxel, say v to one gene, say g and has a weight given by the expression energy of E(v, g) of the gene at the voxel.
We looked for partitions of this weighted bipartite graph into subgraphs such that the weights of the internal edges of the subgraphs are strong compared to the weights of the edges between the subgraphs. This is the isoperimetric problem addressed by the algorithm of [19] (the graph need not be bipartite to apply this algorithm, but since we started with a bipartite graph, each of the subgraphs, or biclusters, returned by the algorithm, is bipartite, and therefore corresponds to a set of voxels and a set of genes).
Given a weighted graph, the algorithm cuts some of the links, thus partitioning the graph into a subset S and its complementaryS, such that the sum of weights in the set of cut edges is minimized relative to the total weight of internal edges in S. The sum of weights in the set of cut edges is analogous to a boundary term, while the total weight of internal edges is analogous to a volume term. In that sense the problem is an isoperimetric optimization problem, and the optimal set S minimizes the isoperimetric ratio ρ over all the possible subgraphs: where the quantity W ij is the weight of the link between vertex i and vertex j. Once S has been worked out, the algorithm can be applied separately to S and its complemen-taryS. This recursive application goes on until the isoperimetric ratio reaches a stopping ratio, representing the highest allowed isoperimetric ratio. This value is a parameter of the algorithm. Rising it results in a higher number of clusters, as it rises the number of acceptable cuts. The implementation of the recursive partition of the bipartite graph in the present toolbox is due to Grady and Schwartz [19].