单细胞多样本整合之R语言版scVI
上个帖子简单介绍了scVI和scANVI，以及其python环境部署，并尝试运行了一个示例数据，详见：
利用reticulate可以在R语言中运行scVI，本期推文对此做一简单介绍。
scVI官网教程在https://docs.scvi-tools.org/en/stable/tutorials/notebooks/harmonization.html，github在https://github.com/scverse/scvi-tools scVI属于比较高级的整合算法（常见于各大高分文章），但仍需强调的是，scVI/scANVI整合分析比较吃资源，如果有GPU的话速度会快5-10倍；如果仅用CPU运行scVI，还需谨慎。因为CPU模式的scVI运行速度挺慢的，特别是针对大样本单细胞数据的运行。
一. 环境部署
虽然使用R语言运行，但仍需部署一个scvi的Python环境，这部分与推文【单细胞多样本整合之scVI和scANVI】的环境部署部分一样：
mamba create -n scvi python=3.9
conda activate scvi
#安装软件
pip install scvi-tools anndata numpy scanpy scib certifi scib-metrics pymde scvi-colab
#Please be sure to install a version of PyTorch that is compatible with your GPU (if applicable).
还需要在Python里运行：
from scvi_colab import install
install()
如果你的服务器有GPU的话，可以安装下面的包进行GPU加速：
mamba install pytorch torchvision torchaudio pytorch-cuda=11.7 -c pytorch -c nvidia -y
mamba install jax jaxlib -c conda-forge -y
如果下述R包没有安装的，需要在R语言中安装：
# install.packages("Seurat")
# install.packages("reticulate")
# install.packages("cowplot")
# install.packages("devtools")
# devtools::install_github("satijalab/seurat-data")
# SeuratData::InstallData("pbmc3k")
# install.packages("https://seurat.nygenome.org/src/contrib/ifnb.SeuratData_3.0.0.tar.gz", repos = NULL, type = "source")
# SeuratData::InstallData("ifnb")
# devtools::install_github("cellgeni/sceasy")
二. 运行示例数据
Step1. 读入数据
示例数据采用经典的ifnb数据集。
## R包加载
library(Seurat)
library(stringr)
library(ggplot2)
library(patchwork)
library(reticulate)
library(sceasy)
library(readr)
library(cowplot)
library(SeuratData)
sc <- import("scanpy", convert = FALSE)
scvi <- import("scvi", convert = FALSE)
### 1.读入数据
data("ifnb")
Step2. 标准化+高变基因选取
### 2. 标准分析流程
seurat.data <- ifnb %>% NormalizeData(normalization.method = "LogNormalize",
scale.factor = 10000) %>%
FindVariableFeatures(selection.method = "vst", nfeatures = 2000)
top2000 <- head(VariableFeatures(seurat.data), 2000)
seurat.data <- seurat.data[top2000]
Step3. 将Seurat数据转化为AnnData对象
关于Seurat数据转化为AnnData对象（R语言单细胞与Python数据互转），很久之前我介绍过一个R/Python包diopy，详见：
这里想和大家分享另一款R包sceasy，加上outFile参数即可输出h5ad的数据，衔接Python非常方便：
### 3.将seurat转为AnnData
adata <- convertFormat(seurat.data,
#outFile = "./python.data/NatCellBiol.Chen.2021.Batch1.h5ad",
from="seurat",
to="anndata",
main_layer="counts",
drop_single_values=FALSE,
#outFile='filename.h5ad'
)
print(adata) # Note generally in Python, dataset conventions are obs x var
AnnData object with n_obs × n_vars = 13997 × 2000
obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'stim', 'seurat_annotations', 'percent.mt'
var: 'vst.mean', 'vst.variance', 'vst.variance.expected', 'vst.variance.standardized', 'vst.variable'
Step4. scVI整合
### 3.整合
# run setup_anndata, use column stim for batch
scvi$model$SCVI$setup_anndata(adata, batch_key = 'stim')
# create the model
model = scvi$model$SCVI(adata)
# train the model
model$train()
# to specify the number of epochs when training:
# model$train(max_epochs = as.integer(400))
# get the latent represenation
latent = model$get_latent_representation()
# put it back in our original Seurat object
latent <- as.matrix(latent)
rownames(latent) = colnames(ifnb)
ifnb[["scvi"]] <- CreateDimReducObject(embeddings = latent,
key = "scvi_",
assay = DefaultAssay(ifnb))
ifnb <- RunUMAP(ifnb, dims = 1:10, reduction = "scvi", n.components = 2)
这里可以学习一下如何使用CreateDimReducObject函数哦。
UMAP可视化：
p1 <- DimPlot(ifnb, reduction = "umap", group.by = "stim", pt.size=2)
p1
image-20230806160515248
FeaturePlot(ifnb, features = c("SELL", "CREM", "CD8A", "GNLY", "CD79A", "FCGR3A",
"CCL2", "PPBP"), min.cutoff = "q9")
image-20230806160549246
FeaturePlot(ifnb, features = c("GNLY", "IFI6"), split.by = "stim", max.cutoff = 3,
cols = c("grey", "red"))
image-20230806160617758
三. 番外
总体感受：R语言版本的scVI还是Python版的运行的舒服一些，因此我建议能使用Python的尽量还是使用Python。另外，以后有志于专注单细胞数据分析的，特别是大样本单细胞数据分析，我觉得还是有必要学习一点Python。
最后，上期推文【单细胞多样本整合之scVI和scANVI 】我使用CPU跑了好几个小时。前几天我借到一台GPU服务器，然后运行同样的示例数据只需要10分钟：
希望有一天能够服务器自由，GPU自由~
- END -
腾讯云开发者
