LRcell: Differential cell type change analysis using Logistic/linear Regression. (original) (raw)
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.21 (BiocManager 1.30.25), R 4.5.0 RC (2025-04-04 r88126)
#> Warning: package(s) not installed when version(s) same as or greater than current; use
#> `force = TRUE` to re-install: 'LRcell'
#> Old packages: 'cobalt'To check whether LRcell package is successfully installed:
library(LRcell)
#> Loading required package: ExperimentHub
#> Loading required package: BiocGenerics
#> Loading required package: generics
#>
#> Attaching package: 'generics'
#> The following objects are masked from 'package:base':
#>
#> as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
#> setequal, union
#>
#> Attaching package: 'BiocGenerics'
#> 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, aperm, append,
#> as.data.frame, basename, cbind, colnames, dirname, do.call,
#> duplicated, eval, evalq, get, grep, grepl, is.unsorted, lapply,
#> mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#> rank, rbind, rownames, sapply, saveRDS, table, tapply, unique,
#> unsplit, which.max, which.min
#> Loading required package: AnnotationHub
#> Loading required package: BiocFileCache
#> Loading required package: dbplyrLRcell 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()
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 15 records
#> # snapshotDate(): 2025-04-11
#> # $dataprovider: Emory University
#> # $species: Mus musculus, Homo sapiens
#> # $rdataclass: Matrix, list
#> # 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
#> ... ...
#> EH5420 | Human PBMC Marker Genes
#> EH6727 | MSigDB C8 MANNO MIDBRAIN
#> EH6728 | MSigDB C8 ZHENG CORD BLOOD
#> EH6729 | MSigDB C8 FAN OVARY
#> EH6730 | MSigDB C8 RUBENSTEIN SKELETAL MUSCLE
## 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.
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 fromDESeq2. 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-01Here, we use Linear Regression:
res <- LRcell(gene.p = example_gene_pvals,
marker.g = NULL,
species = "mouse",
region = "FC",
method = "LiR")
#> 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_CGEPlot 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.
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
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_CGEPlot out the result:
plot_manhattan_enrich(res, sig.cutoff = .05, label.topn = 5)We can clearly find that Microglia (FC_11-1.Microglia) is the most transcriptionally enriched when using Logistic Regression analysis. The result is clear and sound because proliferation and activation of Microglia in the brain is a known feature of Alzheimer’s Disease (AD).
LRcell is used to give a hint on which cell types are potentially involved in diseases when a paired scRNA-seq and bulk RNA-seq experiments are unavailable. Thus, LRcell can be applied to the considerable amount of bulk DEG experiments to shed light on biological discoveries.
Calculate gene enrichment scores from expression dataframe
The gene enrichment score calculation is based on algorithm proposed in (Marques et al. 2016), which takes both enrichment of genes in certain cell types and fraction of cells expressing the gene into consideration. If interested, please take a look at the Supplementary Materials in the original paper.
LRcell_gene_enriched_scores() function takes the gene-by-cell expression matrix and a cell-type annotation as input, which means cell type assignments should be done first. The columns of the expression matrix should be in accordance with cell names in the annotation vector.
Take a randomly-generated data as example.
# print out the generated expression matrix
print(sim.expr)
#> cell1 cell2 cell3 cell4 cell5 cell6 cell7 cell8 cell9 cell10
#> gene1 3 0 2 8 10 6 1 0 0 2
#> gene2 7 5 8 1 0 5 2 3 2 1
#> gene3 8 10 6 7 8 9 5 8 6 8
# print out the cell-type annotation
print(sim.annot)
#> cell1 cell2 cell3 cell4 cell5 cell6
#> "celltype1" "celltype1" "celltype1" "celltype2" "celltype2" "celltype2"
#> cell7 cell8 cell9 cell10
#> "celltype3" "celltype3" "celltype3" "celltype3"This toy example contains 3 genes and 10 cells. As you can tell from the matrix, gene1 is a marker gene of celltype2; gene2 is a marker gene of celltype1; gene3 is a house keeping gene.
# generating the enrichment score
enriched_res <- LRcell_gene_enriched_scores(expr = sim.expr,
annot = sim.annot, parallel = FALSE)
#> Generate enrichment score for each gene..
#>
|
| | 0%
|
|======================= | 33%
|
|=============================================== | 67%
|
|======================================================================| 100%
#> total gene number: 3
enriched_res
#> celltype1 celltype2 celltype3
#> gene1 0.3472222 2.5000000 0.1171875
#> gene2 1.9607843 0.3921569 0.5882353
#> gene3 1.0666667 1.0666667 0.9000000According to the result, it is a gene-by-celltype dataframe indicating the enrichment score of genes in different cell type. Since gene2 is a marker gene of celltype1, thus it has the highest enrichment score.
When choosing marker genes using top 1 as threshold:
marker_res <- get_markergenes(enriched.g = enriched_res,
method = "LR", topn=1)
marker_res
#> $celltype1
#> [1] "gene2"
#>
#> $celltype2
#> [1] "gene1"
#>
#> $celltype3
#> [1] "gene3"