The Allen Institute for Brain Science provides an RNA sequencing (RNA-Seq) data resource for studying transcriptional mechanisms involved in human brain development known as BrainSpan. This resource serves as an additional control in research involving RNA sequencing of human brain tissue. BrainSABER facilitates comparisons of user data with the various developmental stages and brain structures found in the BrainSpan atlas. It extends the SummarizedExperiment class into a self-validating container for user data and produces similarity matrices containing samples of the user’s data analyzed against each sample in the BrainSpan database. These matrices are presented as dynamic heat maps, and include both cosine and Euclidean similarity methods.
BrainSABER requires Bioconductor 3.11 or higher. For more information on Bioconductor, please see their website at https://bioconductor.org. To install Bioconductor and BrainSABER, run the following commands in R or RStudio:
# for R version >= 4.0
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("BrainSABER")After installing, attach the BrainSABER package with the following command:
library(BrainSABER)To leverage the Allen Institute for Brain Science Brain Atlas RNA sequence data, BrainSABER includes a function, ‘buildAIBSARNA’, that will download the data, format it into a RangedSummarizedExperiment, and use biomaRt to add RefSeq identifiers. This function also stores the downloaded zipped file using the BiocFileCache package in a local cached directory, so after the initial run the process should be faster and won’t require an active internet connection. For more information on how BiocFileCache stores data use the link or see the package vignette.
# Extract dataset from online source or from local BiocCache
AIBSARNA <- buildAIBSARNA()Once AIBSARNA has been build, this function does not need to be run again for subsequent work unless the AIBSARNA variable is removed or a new R session is started. If desired, AIBSARNA can be saved as an R data object and reloaded for future work without the need to rebuild it.
The BrainSABER workflow can be used to compare the Allen Institute data to a CellScabbard object. The CellScabbard class extends the SummarizedExperiment class, and in addition requires a sample of the Allen Institute’s data.
To build a CellScabbard object requires the following inputs:
exprsData, a numeric matrix of expression values, where rows are genes, and columns are cells/samples
phenoData, a data frame or DataFrame where rows are cells/samples, and columns are cell/sample attributes (such as cell type, culture condition, day captured, etc.)
featureData, a data frame or DataFrame where rows are features (e.g. genes), and columns are gene attributes, such as gene identifiers, gc content, etc.
AIBSARNA, a SummarizedExperiment or RangedSummarizedExperiment object containing sample AIBSARNA data.
The expression value matrix must have the same number of columns as the phenoData has rows, and it must have the same number of rows as the featureData data frame has rows. Row names of the phenoData object should match the column names of the expression matrix. Row names of the featureData object should match row names of the expression matrix.
For this vignette, we will construct a toy CellScabbard object by extracting 50 samples and 200 genes from AIBSARNA, using the code below:
# Obtain the sample indexes to use for subsetting (not random)
sample_idx <- 1:50 * 10 - 1
# Set the RNG seed for repeatable results
set.seed(8)
# Get the total number of genes available
totalGenes <- nrow(AIBSARNA)
# Sample the indices of 200 random genes
gene_idx <- sample.int(totalGenes, 200)
# Subset AIBSARNA
toy_exprs <- assay(AIBSARNA)[gene_idx, sample_idx]
toy_fd <- rowData(AIBSARNA)[gene_idx, ]
toy_pd <- colData(AIBSARNA)[sample_idx, ]
# Create toy CellScabbard
toySet <- CellScabbard(exprsData = toy_exprs, phenoData = toy_pd,
featureData = toy_fd, AIBSARNA = AIBSARNA,
autoTrim = TRUE)Before we can compare our data to the Allen Institute data, we must remove any genes that are not present in the Allen Institute data. To accomplish that, we need to explore different combinations of gene identifiers for our data and the Allen Institute data to find which combination yields the most matched genes. We can use BrainSABER’s “getExternalVector()” function to extract trimmed vectors from our toySet, and compare lengths to find which one contains the most matches.
# Try comparing different identifiers
length(getExternalVector(toySet, index = 1, AIBSARNA = AIBSARNA,
dataSetId = "gene_id", AIBSARNAid = "gene_id"))## [1] 200length(getExternalVector(toySet, index = 1, AIBSARNA = AIBSARNA,
dataSetId = "ensembl_gene_id",
AIBSARNAid = "ensembl_gene_id"))## [1] 200length(getExternalVector(toySet, index = 1, AIBSARNA = AIBSARNA,
dataSetId = "gene_symbol",
AIBSARNAid = "gene_symbol"))## [1] 200length(getExternalVector(toySet, index = 1, AIBSARNA = AIBSARNA,
dataSetId = "entrez_id", AIBSARNAid = "entrez_id"))## [1] 200length(getExternalVector(toySet, index = 1, AIBSARNA = AIBSARNA,
dataSetId = "refseq_ids", AIBSARNAid = "refseq_ids"))## [1] 200Because the Allen Institute data was originally sequenced using Ensembl identifiers, the “ensembl_id” column of AIBSARNA is the only column with no missing values, making it the ideal choice if the user data also contains ensembl identifiers. If not, the Allen Institute data contains other identifiers. More information is available at their website, http://brainspan.org.
Before we can compare our data to the Allen Institute data, we must remove any genes that are not present in the Allen Institute data. BrainSABER includes a function to trim down user data, which returns a new SummarizedExperiment object. Although this process is done automatically when the user sets the autoTrim argument of CellScabbard() to true, the user can also manually trim their data by specifying the subsetting columns. We will use the “ensembl_gene_id” columns to trim by, as they showed the most matches in the previous step.
trimmed_toySet <- getTrimmedExternalSet(dataSet = toySet,
dataSetId = "ensembl_gene_id",
AIBSARNA = AIBSARNA,
AIBSARNAid = "ensembl_gene_id")We must also filter AIBSARNA to obtain only the genes present in our data set, using BrainSABER’s “getRelevantGenes()” function. Although the CellScabbard constructor does this automatically and stores the results under the ‘relevantGenes’ slot, we will use the “ensembl_gene_id” columns to demonstrate the process manually.
trimmed_AIBSARNA <- getRelevantGenes(data = toySet,
dataSetId = "ensembl_gene_id",
AIBSARNA = AIBSARNA,
AIBSARNAid = "ensembl_gene_id")
# Or extract the results directly from our toySet
autotrim_AIBSARNA <- relevantGenes(toySet)Note that both the autotrimming and autofiltering processes in ‘CellScabbard()’ use the highest-matching column identifier.
Once we have filtered the data to contain only genes present in both the user data and the Allen Institute data, we can obtain a data frame containing similarity scores comparing each sample in the user data to each sample in the Allen Institute data. This can be done using the manually trimmed user and AIBSARNA data, or any trimmed CellScabbard set, which already contains the necessary info in its relevantGenes slot. The getSimScores() function currently supports scoring based on either euclidean distance or cosine similarity, and users should store results in the ‘similarityScores’ slot if using a CellScabbard set. The scores range from 0 to 1, where 1 is a perfect match, and 0 is completely dissimilar.
# Using manually filtered data sets
euc_sim <- getSimScores(data = trimmed_toySet,
relevantGenes = trimmed_AIBSARNA,
similarity_method = "euclidean")
cos_sim <- getSimScores(data = trimmed_toySet,
relevantGenes = trimmed_AIBSARNA,
similarity_method = "cosine")
# Or using the auto-trimmed toySet
auto_euc_sim <- getSimScores(data = toySet, similarity_method = "euclidean")
auto_cos_sim <- getSimScores(data = toySet, similarity_method = "cosine")We can use the similarity scores generated in the previous step to generate a matrix or data frame for each sample in our data to identify which ages and brain structures from the Allen Institute data that sample most closely matches. The getSimMatrix function will return single list which contains the matrices for all samples, and the getSimDataFrame function will return a single list containing data frames for all samples, sorted by the highest score. If passing in a CellScabbard, the functions will automatically extract the relevant genes and similarity scores from the data set, but the user will have to store them there first.
# Using manually filtered data scores
euc_mats <- getSimMatrix(sim_score = euc_sim, relevantGenes = trimmed_AIBSARNA)
euc_df <- getSimDataFrame(sim_score = euc_sim,
relevantGenes = trimmed_AIBSARNA,
similarity_method = "euclidean")
cos_mats <- getSimMatrix(sim_score = cos_sim, relevantGenes = trimmed_AIBSARNA)
cos_df <- getSimDataFrame(sim_score = cos_sim,
relevantGenes = trimmed_AIBSARNA,
similarity_method = "cosine")
# Or using the auto-trimmed data scores
# first store the data in the toySet, then call the similarity functions
similarityScores(toySet) <- auto_euc_sim
auto_euc_mats <- getSimMatrix(data = toySet)
auto_euc_df <- getSimDataFrame(data = toySet, similarity_method = "euclidean")
# to determine cosine similarity, reset the similarityScores data and then call similarity functions
similarityScores(toySet) <- auto_cos_sim
auto_cos_mats <- getSimMatrix(data = toySet)
auto_cos_df <- getSimDataFrame(data = toySet, similarity_method = "cosine")
# Store results of euclidean testing in the toySet
similarityMatrices(toySet) <- auto_euc_mats
similarityDFs(toySet) <- auto_euc_dfTo access the similarity matrices and data frames from the toySet, use the ‘similarityMatrices()’ and ‘similarityDFs()’ accessor methods.
Using R’s head() function, we can see the best matches for a sample in our list of data frames:
head(cos_df[[1]])## age structure_acronym cosine_similarity
## 211 24 pcw VFC 0.003296334
## 194 21 pcw OFC 0.003083213
## 220 35 pcw CBC 0.002850345
## 203 24 pcw S1C 0.002796607
## 186 21 pcw MFC 0.002363169
## 208 24 pcw OFC 0.002333901The matrices can be fed into a heatmapping function to produce graphics displaying the similarity for a sample. we will use the heatmaply package to generate an interactive heatmap:
library(heatmaply)## Loading required package: plotly## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:IRanges':
##
## slice## The following object is masked from 'package:S4Vectors':
##
## rename## The following object is masked from 'package:biomaRt':
##
## select## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layout## Loading required package: viridis## Loading required package: viridisLite##
## ======================
## Welcome to heatmaply version 1.4.0
##
## Type citation('heatmaply') for how to cite the package.
## Type ?heatmaply for the main documentation.
##
## The github page is: https://github.com/talgalili/heatmaply/
## Please submit your suggestions and bug-reports at: https://github.com/talgalili/heatmaply/issues
## You may ask questions at stackoverflow, use the r and heatmaply tags:
## https://stackoverflow.com/questions/tagged/heatmaply
## ======================##
## Attaching package: 'heatmaply'## The following object is masked from 'package:BiocGenerics':
##
## normalizeheatmaply(euc_mats[[1]])In addition to similarity scoring, BrainSABER can also be used to create a list of matrices, one for each sample, indicating whether genes are up-reglated, down-regulated, or normal as compared to the Allen Institute data. The getUNDmatrix() function compares the two data sets, either directly with the “discrete” option, or using the log2 fold-change via “log2fc”, and returns either a numerical matrix or a character matrix where ‘-1’ or ‘D’ indicates down-regulation, ‘1’ or ‘U’ indicates up-regulation, and ‘0’ or ‘N’ indicates normal. The rows of these matrices will be named with the rownames of the user’s trimmed gene expression matrix, and so correspond to the genes of the user data, while the columns of the matrices will be named for the columns in the trimmed Allen Institute data gene expression matrix, and thus correspond to each sample in the Allen Institute data
und_num <- getUNDmatrix(dataSet = trimmed_toySet,
relevantGenes = trimmed_AIBSARNA,
method = "log2fc",
matrix_type = "num")
und_num[[1]][1:10, 1:10]## Sample_1 Sample_2 Sample_3 Sample_4 Sample_5 Sample_6
## Gene_13620 0.7805937 1.001816 2.509564 1.620752 1.423441 2.625176
## Gene_42722 -Inf -Inf -Inf -Inf -Inf -Inf
## Gene_47116 -Inf -Inf -Inf -Inf -Inf -Inf
## Gene_42238 -5.0511187 -5.893638 -3.034275 -3.457771 -4.815736 -3.339929
## Gene_5853 -Inf -Inf -Inf -Inf -Inf -Inf
## Gene_35715 -Inf -Inf -Inf -Inf -Inf -Inf
## Gene_17332 -4.4062581 -4.222253 -4.523927 -4.661644 -4.058217 -4.120087
## Gene_34152 -1.4152096 -1.099079 -1.582914 -1.206143 -1.314688 -1.437514
## Gene_26734 -Inf -Inf -Inf -Inf -Inf -Inf
## Gene_14089 -Inf -Inf -Inf -Inf -Inf -Inf
## Sample_7 Sample_8 Sample_9 Sample_10
## Gene_13620 1.680324 2.982117 1.348810 1.266450
## Gene_42722 -Inf -Inf -Inf -Inf
## Gene_47116 -Inf -Inf -Inf -Inf
## Gene_42238 -4.791561 -1.607487 -4.014670 -3.350593
## Gene_5853 -Inf -Inf -Inf -Inf
## Gene_35715 -Inf -Inf -Inf -Inf
## Gene_17332 -4.356516 -4.499664 -3.905202 -3.419984
## Gene_34152 -1.162350 -1.290172 -1.361178 -1.663125
## Gene_26734 -Inf -Inf -Inf -Inf
## Gene_14089 -Inf -Inf -Inf -Infund_char <- getUNDmatrix(dataSet = trimmed_toySet,
relevantGenes = trimmed_AIBSARNA,
method = "log2fc",
matrix_type = "char")
und_char[[1]][1:10, 1:10]## Sample_1 Sample_2 Sample_3 Sample_4 Sample_5 Sample_6
## Gene_13620 0.7805937 1.001816 2.509564 1.620752 1.423441 2.625176
## Gene_42722 -Inf -Inf -Inf -Inf -Inf -Inf
## Gene_47116 -Inf -Inf -Inf -Inf -Inf -Inf
## Gene_42238 -5.0511187 -5.893638 -3.034275 -3.457771 -4.815736 -3.339929
## Gene_5853 -Inf -Inf -Inf -Inf -Inf -Inf
## Gene_35715 -Inf -Inf -Inf -Inf -Inf -Inf
## Gene_17332 -4.4062581 -4.222253 -4.523927 -4.661644 -4.058217 -4.120087
## Gene_34152 -1.4152096 -1.099079 -1.582914 -1.206143 -1.314688 -1.437514
## Gene_26734 -Inf -Inf -Inf -Inf -Inf -Inf
## Gene_14089 -Inf -Inf -Inf -Inf -Inf -Inf
## Sample_7 Sample_8 Sample_9 Sample_10
## Gene_13620 1.680324 2.982117 1.348810 1.266450
## Gene_42722 -Inf -Inf -Inf -Inf
## Gene_47116 -Inf -Inf -Inf -Inf
## Gene_42238 -4.791561 -1.607487 -4.014670 -3.350593
## Gene_5853 -Inf -Inf -Inf -Inf
## Gene_35715 -Inf -Inf -Inf -Inf
## Gene_17332 -4.356516 -4.499664 -3.905202 -3.419984
## Gene_34152 -1.162350 -1.290172 -1.361178 -1.663125
## Gene_26734 -Inf -Inf -Inf -Inf
## Gene_14089 -Inf -Inf -Inf -Inf# Or using the auto-trimmed toySet
auto_und_num <- getUNDmatrix(dataSet = toySet, method = "log2fc", matrix_type = "num")
auto_und_char <- getUNDmatrix(dataSet = toySet, method = "log2fc", matrix_type = "char")sessionInfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] heatmaply_1.4.0 viridis_0.6.2
## [3] viridisLite_0.4.1 plotly_4.10.0
## [5] ggplot2_3.3.6 BrainSABER_1.8.0
## [7] SummarizedExperiment_1.28.0 Biobase_2.58.0
## [9] GenomicRanges_1.50.0 GenomeInfoDb_1.34.0
## [11] IRanges_2.32.0 S4Vectors_0.36.0
## [13] BiocGenerics_0.44.0 MatrixGenerics_1.10.0
## [15] matrixStats_0.62.0 biomaRt_2.54.0
## [17] BiocStyle_2.26.0
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-3 ellipsis_0.3.2 lsa_0.73.3
## [4] XVector_0.38.0 farver_2.1.1 SnowballC_0.7.0
## [7] bit64_4.0.5 AnnotationDbi_1.60.0 fansi_1.0.3
## [10] xml2_1.3.3 codetools_0.2-18 cachem_1.0.6
## [13] knitr_1.40 jsonlite_1.8.3 dbplyr_2.2.1
## [16] png_0.1-7 shiny_1.7.3 BiocManager_1.30.19
## [19] compiler_4.2.1 httr_1.4.4 assertthat_0.2.1
## [22] Matrix_1.5-1 fastmap_1.1.0 lazyeval_0.2.2
## [25] cli_3.4.1 later_1.3.0 htmltools_0.5.3
## [28] prettyunits_1.1.1 tools_4.2.1 gtable_0.3.1
## [31] glue_1.6.2 GenomeInfoDbData_1.2.9 reshape2_1.4.4
## [34] dplyr_1.0.10 rappdirs_0.3.3 Rcpp_1.0.9
## [37] jquerylib_0.1.4 vctrs_0.5.0 Biostrings_2.66.0
## [40] crosstalk_1.2.0 iterators_1.0.14 xfun_0.34
## [43] stringr_1.4.1 mime_0.12 lifecycle_1.0.3
## [46] XML_3.99-0.12 dendextend_1.16.0 ca_0.71.1
## [49] zlibbioc_1.44.0 scales_1.2.1 TSP_1.2-1
## [52] hms_1.1.2 promises_1.2.0.1 RColorBrewer_1.1-3
## [55] yaml_2.3.6 curl_4.3.3 memoise_2.0.1
## [58] gridExtra_2.3 sass_0.4.2 stringi_1.7.8
## [61] RSQLite_2.2.18 foreach_1.5.2 seriation_1.4.0
## [64] filelock_1.0.2 rlang_1.0.6 pkgconfig_2.0.3
## [67] bitops_1.0-7 evaluate_0.17 lattice_0.20-45
## [70] purrr_0.3.5 htmlwidgets_1.5.4 labeling_0.4.2
## [73] bit_4.0.4 tidyselect_1.2.0 plyr_1.8.7
## [76] magrittr_2.0.3 bookdown_0.29 R6_2.5.1
## [79] generics_0.1.3 DelayedArray_0.24.0 DBI_1.1.3
## [82] pillar_1.8.1 withr_2.5.0 KEGGREST_1.38.0
## [85] RCurl_1.98-1.9 tibble_3.1.8 crayon_1.5.2
## [88] utf8_1.2.2 BiocFileCache_2.6.0 rmarkdown_2.17
## [91] progress_1.2.2 grid_4.2.1 data.table_1.14.4
## [94] blob_1.2.3 digest_0.6.30 webshot_0.5.4
## [97] xtable_1.8-4 tidyr_1.2.1 httpuv_1.6.6
## [100] munsell_0.5.0 registry_0.5-1 bslib_0.4.0