scGEM is a nested tree-structured generative model that identifies subtype-specific and -shared gene coexpressing modules (GEMs) in scRNAseq data.
v0.1.4: enable setting species
in initTree
function.
Currently, scGEM works on MAC/Linux. For faster implementation on M1 or M2 MAC, you need to switch your built-in BLAS library to Apple's vecLib. (https://mpopov.com/blog/2021/10/10/even-faster-matrix-math-in-r-on-macos-with-m1/).
devtools::install_github("hansolo-bioinfo/scGEM")
?scGEM::initTree
?scGEM::minibatchInfer
Note: if the Error in Lazy database
pops up, try .rs.restartR()
. For more information, please see (r-lib/devtools#1980).
To test the scGEM performance, we will use processed peripheral blood mononuclear cells (PBMC) 3K dataset as our toy case.
SeuratData::InstallData("pbmc3k")
library(Seurat)
library(SeuratData)
library(magrittr)
library(scGEM)
pbmc3k.final <- LoadData("pbmc3k", type = "pbmc3k.final")
set.seed(111)
X <- pbmc3k.final@assays$RNA@counts
# you can use raw count/log2/or other normalized data for the input
# just make it non-negative and no NAs in the sparse matrix
X@x <- log2(X@x + 1) # log-transform the count data
# for binary learning
# X@x <- ifelse(X@x > 0, 1, 0)
# assuming create a 5-4-3 tree structure
scGEM_init <- initTree(
X,
num_topics = c(5, 4, 3),
blacklist_genes = "^MT|^RP|^HSP|B2M|MALAT1",
rm.housekeeping_genes = TRUE
)
For one-batch training, please set the batch_size
equals to the sample size:
# one-batch
scGEM_trained <- minibatchInfer(
X,
scGEM_init,
model_paras = list(b0 = 0.01,
g1 = 5,
g2 = 1,
g3 = 1/3,
g4 = 2/3),
max_est = 50,
subtree_size = c(1, 20),
batch_size = ncol(X),
n_epoch = 50,
learning_rate = 0.01
)
For mini-batch training (recommended), please do:
# mini-batch
scGEM_trained_mb <- minibatchInfer(
X,
scGEM_init,
model_paras = list(b0 = 0.01,
g1 = 5,
g2 = 1,
g3 = 1/3,
g4 = 2/3),
max_est = 50,
subtree_size = c(1, 20),
batch_size = 1000,
n_epoch = 50,
learning_rate = 0.01
)
Compare the total likelihood during each epoch in one-batch and mini-batch modes.
# total likelihood
plot(scGEM_trained$likelihood, col = "blue", type = "b",
ylab = "Likelihood", xlab = "Epoch")
points(scGEM_trained_mb$likelihood, col = "red", type = "b")
legend("bottomright", legend = c("One-batch", "Mini-batch"),
lty = 1, col = c("blue", "red"))
![Screenshot 2023-08-18 at 10 57 25](https://private-user-images.githubusercontent.com/65295899/261658586-334c3181-cb3b-4e41-a492-43ad2ef57312.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MjM2MTk2MTYsIm5iZiI6MTcyMzYxOTMxNiwicGF0aCI6Ii82NTI5NTg5OS8yNjE2NTg1ODYtMzM0YzMxODEtY2IzYi00ZTQxLWE0OTItNDNhZDJlZjU3MzEyLnBuZz9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNDA4MTQlMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjQwODE0VDA3MDgzNlomWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPWZkMzU2ZGYxMDdkYTM1Y2Y4YjYzNTlkODNkNjhjOWUwNzFmODRkMTc0MjZhNGZkMGYzN2M5NjQyZTMzMjUxOTImWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0JmFjdG9yX2lkPTAma2V5X2lkPTAmcmVwb19pZD0wIn0._xtcrkk67J2K4QRyaSg3TXbDfKqCq9aOLALlxr40YUw)
To get the raw values:
gene_over_gem <- scGEM_trained_mb$centroids # 10640 x 85
gem_over_cell <- scGEM_trained_mb$count_matrix # 85 x 2638
tree_relation <- scGEM_trained_mb$tree # 86 x 2
gene_name <- scGEM_trained_mb$gene
cell_name <- scGEM_trained_mb$cell
gem_name <- paste0("GEM", 1:85)
rownames(gene_over_gem) <- gene_name
colnames(gene_over_gem) <- gem_name
rownames(gem_over_cell) <- gem_name
colnames(gem_over_cell) <- cell_name
To get a relative proportion like in other topic models:
gene_over_gem <- apply(gene_over_gem, 2, function(x) x/sum(x))
gem_over_cell <- apply(gem_over_cell, 2, function(x) x/sum(x))
Plot gene distribution for GEM 1
barplot(sort(gene_over_gem[, 1], decreasing = T), ylab = "weight", xlab = "gene", main = "GEM 1")
![Screenshot 2023-08-18 at 12 28 38](https://private-user-images.githubusercontent.com/65295899/261661957-a191f3a2-ac96-4fbd-a91c-a5d9ace66aaa.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MjM2MTk2MTYsIm5iZiI6MTcyMzYxOTMxNiwicGF0aCI6Ii82NTI5NTg5OS8yNjE2NjE5NTctYTE5MWYzYTItYWM5Ni00ZmJkLWE5MWMtYTVkOWFjZTY2YWFhLnBuZz9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNDA4MTQlMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjQwODE0VDA3MDgzNlomWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPTdlOGZhM2UzYjMxZDA2OGJmNWEyMDc0NzliZjk4NWE4NTJiZTIyMjU2MzA2ZTFhNmE2OGIxZjNlYjQzNzc4ZjYmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0JmFjdG9yX2lkPTAma2V5X2lkPTAmcmVwb19pZD0wIn0.Gx0XS80AbOKf0ygEJTI0v_XM978miiVR_WSmg5rMNSQ)
Plot GEM distribution for the 100th cell
barplot(gem_over_cell[, 100], ylab = "weight", xlab = "GEM",
main = "100th cell: AAGATTACCGCCTT")
![Screenshot 2023-08-18 at 12 31 24](https://private-user-images.githubusercontent.com/65295899/261662359-adced07b-67a8-47bd-90db-46492ebfad6e.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MjM2MTk2MTYsIm5iZiI6MTcyMzYxOTMxNiwicGF0aCI6Ii82NTI5NTg5OS8yNjE2NjIzNTktYWRjZWQwN2ItNjdhOC00N2JkLTkwZGItNDY0OTJlYmZhZDZlLnBuZz9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNDA4MTQlMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjQwODE0VDA3MDgzNlomWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPTEwYTc1MDE5ZTBlZTlhMDQwOWQxZDBjYTVmNWU0OTQzNGRjYWYxODdjYzQ3NzAxOTk3MWNmYjYxM2U4OTk4MDEmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0JmFjdG9yX2lkPTAma2V5X2lkPTAmcmVwb19pZD0wIn0.FXN5XtLvhv9zL20wY5bb-zFY_u15ALmLqEDRxUTbYHg)