1.3: Plot the Brain

Next we would like to see if neighboring brain areas share SGSE estimates. To do this we will use CreateBrain to merge all of the microdissected tissues returned from testEnrich so that each general brain area present within the Allen Brain Atlas reference map is supported by the respective set of underlying tissues.

As mentioned above, note that the general brain areas within the Allen Brain Atlas reference map and the microdissected tissues present within the Allen Brain Atlas transcription datasets do not have a 1:1 relationship. Therefore some areas are not supported by transcriptional information in the reference map while other areas are supported by more than one tissue.

CreateBrain takes the results from SpatialEnrichment and testEnrich as input. There are two additional arguments to consider.

CreateBrain arguments:

• slice: Image from within either the developing or adult brain set to plot. The choices for both the developing and adult brain are 1-10.
• pcut: p-value threshold for considering a microdissected tissue as being enriched or depleted

Since we initially found the LHAa interesting, let’s see what slice of the developing brain contains general brain areas composed of the LHAa.

tis_in_region(composite, "LHAa")
#>   slice areas
#> 1     6   THM
#> 2     7   THM

Slices 6 and 7 both have relevant information for LHAa in the thalamus (THM). Let’s use slice 6 for downstream analysis.

composite <- CreateBrain(composite, res, slice = 5, pcut = 0.05)

Once the microdissected tissues have been merged into the corresponding general brain areas, the image can be plotted.

PlotBrain(composite)

From this image it looks like the thalamus is indeed enriched for the vth gene set.Interestingly surrounding areas are also enriched for this set. By referencing the Allen Brain Atlas we can see that those are general brain areas like the Putamen (Pu) and Caudate Nucleus (Ca).

Note that you can either reference the Allen Brain Atlas directly or use available_areanames to identify the list of available general brain area names to query.

available_areanames(composite, slice = 6)
#>  [1] "BNST" "CP"   "Ca"   "FWM"  "HTH"  "IZ"   "LGE"  "MGE"  "MZ"   "Pu"  
#> [11] "R"    "SG"   "SN"   "SP"   "SZ"   "THM"  "VA"   "VZ"   "alCP" "peCP"
#> [21] "peMZ"

Let’s look more closely at the general brain area of the Putamen (Pu). Using tis_set we can work backwards and identify which microdissected tissues support the final color code in a given area of interest.

tis_set(composite, area.name = "Pu", slice = 6)
#>   abrev               structure_name
#> 1   Pld laterodorsal part of putamen
#> 2   Pmv medioventral part of putamen
#> 3   PuC               caudal putamen
#> 4   PuR              rostral putamen

And now we can go back and observe the gene overlap in a microdissected tissue within the Putamen. Here we focus in on the medioventral part of putamen (Pmv)

vth_pld_overlap <- GetGenes(vth, composite, tissue_abrev = "Pmv")
length(vth_pld_overlap)
#> [1] 154
##vth_go <- enrichGO(gene = vth_pld_overlap,
##                   OrgDb = org.Hs.eg.db,
##                   keytype = 'SYMBOL',
##                   pvalueCutoff = 0.05, 
##                   qvalueCutoff = 0.05)

##dotplot(vth_go, showCategory=30)

The GO term “hormone activity” has showed up in both the comparison with LHAa and Pmv. Perhaps there is a shared pathway across these tissues. Let’s use whichtissues to see if the genes with hormone activity in the Pmv are also present in the LHAa.

#grab the genes associated with hormone activity
#vth_go2 <- data.frame(vth_go)
#vth_match <- vth_go2$Description == "hormone activity"
#vth_pmv_hormone <- vth_go2[vth_match, "geneID"]
#vth_pmv_hormone <- unlist(strsplit(vth_pmv_hormone,"/"))


#identify the presence of these genes across all tissues
#vth_pmv_hormone_tis <- whichtissues(vth_pmv_hormone, refset = "developing")
#vth_pmv_hormone_tis[,1:10]


#which genes are present (1) in the LHAa
#all <- vth_pmv_hormone_tis[,"LHAa"]
#vth_pmv_hormone[ all == 1]

Indeed all of the genes involved in hormone activity in the Pmv are also expressed in the LHAa.

Using the above mentioned functions, brainImageR has provided us with a easy way to query SGSE with respect to the postmortem human brain. We calculated the significance of the enrichment of gene overlaps and visualized these enrichments with images of the human brain. Now let’s move on to predicting developmental time from a dataset.

SECTION 2: DEVELOPMENTAL TIME POINT PREDICTION

2.0: Preparing the data

An additional functionality of brainImageR is predicting the developmental timepoint of the sample with reference to the postmortem human brain. This analysis takes a normalized expression matrix as input. The columns of this matrix should be samples, and the rows should be gene names in the Gene Symbol format. An example is provided here in dat.

data(dat)
dim(dat)
#> [1] 10192     6
head(rownames(dat))
#> [1] "A2M"   "NAT2"  "ACAT1" "ADA"   "ADRB2" "ADRB3"
head(colnames(dat))
#> [1] "NPC.1"     "NPC.2"     "NPC.3"     "NPC.4"     "Neurons.5" "Neurons.6"

state <- do.call("rbind", strsplit(colnames(dat), ".", fixed = TRUE))[,1]

2.1: Predict developmental time

In the simplest case, one can predict developmental time with respect to the full Allen Brain Atlas dataset by providing predict_time. It has been shown across multiple studies that neurons derived in-vitro primarily recapitulate prenatal development. Therefore we’ve provided additional flexibility in predict_time to obtain higher resolution within this prenatal window. This is done by restricting the analysis to samples within the Allen Brain Atlas that come from a given class of data. You can see which datasets that we’ve precomputed using availabledatasets.

#availabledatasets()

Let’s start with the default settings which will use all of the available samples to perform the temoral prediction.

time <- predict_time(dat,  minage = 8, maxage = 40)

PlotPred will show the predictions overlaid on top of predictions against the Allen Brain Atlas reference set (aba) to show how the model performs on the starting dataset. This plot enables the user to have an intuitive understanding of how the model is performing

PlotPred(time)

As we can see, using all samples, the prediction is linear across all timepoints. As with most in vitro studies of neurons, our samples are clustering with the prenatal samples. Although our model is good at predicting time over all ages, there is the issue that the prenatal samples are predicted within a small window of each other. To provide higher resolution within this time period lets restrict the analysis to use only prenatal samples in the prediction.

If we select “prenatal” this will restrict the predictions to only those samples that are < 40 weeks post-conception.

#time <- predict_time(dat,dataset = "prenatal")
#PlotPred(time)

Now we can see that the model has higher resolution within the prenatal timepoints. Note that this model should not be used for postnatal samples, as (conversely to the default model) it has a reduced ability to resolve differences between postnatal timepoints.

Now that we have a model that is useful for examining samples with respect to prenatal development, lets compare the NPCs to the Neurons.

#time2 <- data.frame(pred_age = time@pred_age, state)
#time2$state <- factor(time2$state, c("NPC","Neurons"))

#ggplot(time2, aes(state,pred_age))+
#  geom_boxplot()+
#  ylab("Predicted time, weeks post-conception")+
#  labs(title = "predicted developmental time,\ndataset = prenatal")

As expected the neural progenitor cells are predicted to be younger in developmental time than the neurons. In the manuscript associated with this package there are additional examples of predicting developmental time on post-mortem and in vitro-derived neurons showing the utility of predict_time across datasets and types.

There are additional options in predict_time that allow the user to further refine their analysis.

predict_time arguments:

• dat: A normalized expression matrix where columns are samples, and rows are genes defined by Gene Symbol.
• dataset: If desired, the user can work from a pre-computed model. Use availabledatasets to see what options are available. default = “all”.
• genelist: If you would like to restrict the predictions to a certain subset of genes you can retain your full dataset in the argument dat and input the restricted gene list here.
• minage: Minimum weeks post-conception for running the prediction.
• maxage: Maximum weeks post-conception for running the prediction.
• tissue: Tissue of interest to restrict the dataset to.
• minrsq: (0-1). minrsq is correlated with the predictive strength of the model. Reduce the value to allow weaker models to be calculated.

Using these additional options we have more flexibility in our predictions. For example, let’s say that we have generated neurons that we expect to be functionally similar to the adult amygdala and therefore would prefer a model that most closely represents that profile. We can set minage to 40 and tissue to AMY and recalculate our predictions.

#time <- predict_time(dat,minage = 40, tissue = "AMY")