Contents

1 Introduction

Single-cell RNA-sequencing (scRNA-seq) technologies have revealed cell heterogeneity. Marker gene expressions are considered as the most intuitive and informative measurements distinguishing different cell types. Although several computational methods have been proposed to do marker gene selection from clusters or cell types, none of them applied marker gene information on bulk RNA-seq experiments to find cell type or cluster enrichment. Here, we present LRcell package, which uses Logistic/Linear Regression to identify the most transcriptionally enriched cell types or clusters when applying cell marker information to a bulk experiment with p-values measurements of differentially expressed genes.

Pre-Loaded Marker genes information: This package offers marker genes information in 1 human PBMC data, 1 human brain region and 9 mouse brain regions.

2 Standard Workflow

2.1 Installation

This is a R Bioconductor package and it can be installed by using BiocManager.

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager") ## this will install the BiocManager package
BiocManager::install("LRcell")
#> Bioconductor version 3.13 (BiocManager 1.30.15), R 4.1.0 (2021-05-18)
#> Warning: package(s) not installed when version(s) same as current; use `force = TRUE` to
#>   re-install: 'LRcell'

To check whether LRcell package is successfully installed:

library(LRcell)
#> Loading required package: ExperimentHub
#> Loading required package: BiocGenerics
#> Loading required package: parallel
#> 
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:parallel':
#> 
#>     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
#>     clusterExport, clusterMap, parApply, parCapply, parLapply,
#>     parLapplyLB, parRapply, parSapply, parSapplyLB
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#> 
#>     Filter, Find, Map, Position, Reduce, anyDuplicated, append,
#>     as.data.frame, basename, cbind, colnames, dirname, do.call,
#>     duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
#>     lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
#>     pmin.int, rank, rbind, rownames, sapply, setdiff, sort, table,
#>     tapply, union, unique, unsplit, which.max, which.min
#> Loading required package: AnnotationHub
#> Loading required package: BiocFileCache
#> Loading required package: dbplyr

2.2 LRcell usage

Once we have LRcell package loaded, we can start using it to analyze the transcriptional engagement of cell types or clusters. LRcell takes both single-cell marker genes list and p-values of bulk DE genes as input to calculate the enrichment of cell-type specific marker genes in the ranked DE genes.

As mentioned above, LRcell provides single-cell marker genes list in 1 human PBMC data, 1 human brain region (Prefrontal Cortex) and 9 mouse brain regions (Frontal Cortex, Cerebellum, Globus Pallidus, Hippocampus, Entopeduncular, Posterior Cortex, Striatum, Substantia Nigra and Thalamus).

  • The human PBMC data comes from volunteers enrolled in an HIV vaccine trial (Hao et al. 2020) at time point of day 0.
  • The human brain data comes from control samples in Major Depressive Disorder studies. (Nagy et al. 2020).
  • The mouse data comes from the whole brain single-cell RNA-seq experiments(Saunders et al. 2018). Another resource for this dataset is DropViz.

The data is stored in another Bioconductor ExperimentHub package named LRcellTypeMarkers. Users can access the data through ExperimentHub by:

## for installing ExperimentHub
# BiocManager::install("ExperimentHub")

## query data
library(ExperimentHub)
eh <- ExperimentHub::ExperimentHub()
#> snapshotDate(): 2021-05-05
eh <- AnnotationHub::query(eh, "LRcellTypeMarkers")  ## query for LRcellTypeMarkers package
eh  ## this will list out EH number to access the calculated gene enrichment scores
#> ExperimentHub with 11 records
#> # snapshotDate(): 2021-05-05
#> # $dataprovider: Emory University
#> # $species: Mus musculus, Homo sapiens
#> # $rdataclass: Matrix
#> # additional mcols(): taxonomyid, genome, description,
#> #   coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
#> #   rdatapath, sourceurl, sourcetype 
#> # retrieve records with, e.g., 'object[["EH4548"]]' 
#> 
#>            title                                
#>   EH4548 | Mouse Frontal Cortex Marker Genes    
#>   EH4549 | Mouse Cerebellum Marker Genes        
#>   EH4550 | Mouse Entopeduncular Marker Genes    
#>   EH4551 | Mouse Globus Pallidus Marker Genes   
#>   EH4552 | Mouse Posterior Cortex Marker Genes  
#>   ...      ...                                  
#>   EH4554 | Mouse Substantia Nigra Marker Genes  
#>   EH4555 | Mouse Thalamus Marker Genes          
#>   EH4556 | Mouse Hippocampus Marker Genes       
#>   EH4557 | Human pre-Frontal Cortex Marker Genes
#>   EH5420 | Human PBMC Marker Genes

## get mouse brain Frontal Cortex enriched genes
enriched.g <- eh[["EH4548"]]
#> see ?LRcellTypeMarkers and browseVignettes('LRcellTypeMarkers') for documentation
#> loading from cache
marker.g <- get_markergenes(enriched.g, method="LR", topn=100)

Users are also encouraged to process a read count matrix with cell annotation information into a gene enrichment scores matrix.

enriched.g <- LRcell_gene_enriched_scores(expr, annot, power=1, parallel=TRUE, n.cores=4)

Here, expr is a read count matrix with rows as genes and columns as cells. annot is a named-vector with names as cell names (which is in accordance with the column names of expr) and values as annotated cell types. power is a hyper-parameter controlling how much penalty for the proportion of cells expressing certain gene. parallel and n.cores are two parameters for executing function in parallel to accelerate the calculation.

2.2.1 Directly indicate species and brain region in LRcell

Compared to processing data yourself, a much easier way is to indicate species and brain region or tissue. In this way, marker genes are extracted from ExperimentHub accordingly. For example, we can use mouse Frontal Cortex marker genes to do LRcell analysis on the example bulk experiment(Swarup et al. 2019). (The example contains 23, 420 genes along with p-values calculated from DESeq2. Data is processed from a mouse Alzheimer’s disease model (GEO: GSE90693), which is 6 months after treatment in Frontal Cortex brain region.)

# load example bulk data
data("example_gene_pvals")
head(example_gene_pvals, n=5)
#>         Xkr4          Rp1        Sox17       Mrpl15       Lypla1 
#> 1.742186e-05 4.103134e-02 9.697389e-02 1.206500e-02 6.609558e-01

Here, we use Linear Regression:

res <- LRcell(gene.p = example_gene_pvals,
              marker.g = NULL,
              species = "mouse",
              region = "FC",
              method = "LiR")
#> snapshotDate(): 2021-05-05
#> see ?LRcellTypeMarkers and browseVignettes('LRcellTypeMarkers') for documentation
#> loading from cache
FC_res <- res$FC
# exclude leading genes for a better view
sub_FC_res <- subset(FC_res, select=-lead_genes)
head(sub_FC_res)
#>                        ID genes_num          coef   p.value       FDR
#> 1  FC_1-1.Interneuron_CGE        98  6.065974e-04 0.2365456 0.4088820
#> 2 FC_1-10.Interneuron_CGE        92 -7.992953e-05 0.9154947 0.9515587
#> 3 FC_1-11.Interneuron_CGE        95  1.183490e-06 0.9987899 0.9987899
#> 4  FC_1-2.Interneuron_CGE        98  4.642689e-04 0.3418545 0.5114907
#> 5  FC_1-3.Interneuron_CGE        85 -2.765110e-04 0.6647116 0.8273606
#> 6  FC_1-4.Interneuron_CGE        98  5.427872e-06 0.9950063 0.9987899
#>         cell_type
#> 1 Interneuron_CGE
#> 2 Interneuron_CGE
#> 3 Interneuron_CGE
#> 4 Interneuron_CGE
#> 5 Interneuron_CGE
#> 6 Interneuron_CGE

Plot out the result:

plot_manhattan_enrich(FC_res, sig.cutoff = .05, label.topn = 5)

According to the result, when using enrichment scores as a predictor variable in Linear Regression, one cluster of Astrocytes (FC_8-2.Astrocytes) is the mostly enriched along with Microglia (FC_11-1.Microglia). (Note that although cluster FC_11-4 was annotated as unknown in the publication, according to our research, FC_11-* belong to Microglia/Macrophage cell types.) Recent publications have shown that Astrocytes are involved in Alzheimer’s Disease.

2.2.2 Marker gene download and do LRcell analysis

Marker gene list downloading example (mouse, Frontal Cortex, Logistic Regression):

library(ExperimentHub)
eh <- ExperimentHub::ExperimentHub()  ## use ExperimentHub to download data
#> snapshotDate(): 2021-05-05
eh <- query(eh, "LRcellTypeMarkers")
enriched_genes <- eh[["EH4548"]]  # use title ID which indicates FC region
#> see ?LRcellTypeMarkers and browseVignettes('LRcellTypeMarkers') for documentation
#> loading from cache
# get marker genes for LRcell in logistic regression
FC_marker_genes <- get_markergenes(enriched_genes, method="LR", topn=100)

# to have a glance of the marker gene list
head(lapply(FC_marker_genes, head))
#> $`FC_1-1.Interneuron_CGE`
#> [1] "Npy"    "Pde11a" "Kit"    "Gad2"   "Fam46a" "Pnoc"  
#> 
#> $`FC_1-10.Interneuron_CGE`
#> [1] "Krt73"   "Npas1"   "Yjefn3"  "Dlx6os1" "Igf1"    "Sln"    
#> 
#> $`FC_1-11.Interneuron_CGE`
#> [1] "Htr3a"     "Tnfaip8l3" "Npy2r"     "Adarb2"    "Npas1"     "Crabp1"   
#> 
#> $`FC_1-2.Interneuron_CGE`
#> [1] "Ndnf"  "Kit"   "Reln"  "Npy"   "Cplx3" "Gad1" 
#> 
#> $`FC_1-3.Interneuron_CGE`
#> [1] "Calb2"  "Crh"    "Nr2f2"  "Dlx1"   "Penk"   "Pcdh18"
#> 
#> $`FC_1-4.Interneuron_CGE`
#> [1] "Htr3a"    "Tac2"     "Vip"      "Crh"      "Crispld2" "Adarb2"

Then, we can run LRcell analysis by using LRcellCore() function using Logistic Regression.

res <- LRcellCore(gene.p = example_gene_pvals,
           marker.g = FC_marker_genes,
           method = "LR", min.size = 5, 
           sig.cutoff = 0.05)
#> Warning: glm.fit: algorithm did not converge
## curate cell types
res$cell_type <- unlist(lapply(strsplit(res$ID, '\\.'), '[', 2))
head(subset(res, select=-lead_genes))
#>                        ID genes_num        coef odds_ratio     p.value
#> 1  FC_1-1.Interneuron_CGE        98 0.011908189   1.076812 0.006035417
#> 2 FC_1-10.Interneuron_CGE        92 0.006805227   1.043199 0.358384832
#> 3 FC_1-11.Interneuron_CGE        95 0.006160366   1.039027 0.441406086
#> 4  FC_1-2.Interneuron_CGE        98 0.009983947   1.064012 0.039942444
#> 5  FC_1-3.Interneuron_CGE        85 0.004671113   1.029455 0.651009137
#> 6  FC_1-4.Interneuron_CGE        98 0.006406322   1.040616 0.400647777
#>          FDR       cell_type
#> 1 0.01357969 Interneuron_CGE
#> 2 0.41470245 Interneuron_CGE
#> 3 0.48316072 Interneuron_CGE
#> 4 0.05882433 Interneuron_CGE
#> 5 0.66749038 Interneuron_CGE
#> 6 0.45707704 Interneuron_CGE

Plot out the result:

plot_manhattan_enrich(res, sig.cutoff = .05, label.topn = 5)