单细胞转录组 | 多样本处理与锚定法整合
前言
前几期我们介绍了对单个样本进行处理，本次我们介绍如何处理多个样本以及如何对多样本进行整合矫正。
整合是将多次实验的数据进行整合区分。目的是尽可能地消除测序深度和批次效应的影响，让不同样本均匀地分布在不同的cluster中，使不同的样本之间具有很好的可比性。
本次我们选取单细胞转录组 | GEO数据库介绍及数据下载中的BC21和BC3使用锚定进行多样本整合。
本文框架
1. 安装包
如果已经安装，此步请跳过。
install.packages('Seurat')
install.packages('dplyr')
install.packages('tidyverse')
install.packages('patchwork')
2. 加载包
library(Seurat)
library(dplyr)
library(tidyverse)
library(patchwork)
3. 设置工作路径
setwd("D:/sc-seq/")
请根据自己数据的存放位置自定义路径。
工作路径下存放了需要读取的10×数据文件夹：BC3和BC21。
4. 创建文件的向量
创建需要读取的多样本名向量并命名。
## 将文件夹中的文件名存到"dir_name"中
dir_name <- list.files()
## 查看"dir_name"
dir_name
#[1] "BC21" "BC3"
## 为"dir_name"赋名
names(dir_name) = c('BC21', 'BC3')
## 查看赋名后的"dir_name"
dir_name
# BC21 BC3
#"BC21" "BC3"
5.批量读取数据并创建Seurat对象
CreateSeuratObject函数格式：CreateSeuratObject(counts,project = "SeuratObject",min.cells = 0,min.features = 0)
counts：表达矩阵（原始未标准化的数据，细胞作为列，基因作为行）；
min.cells：指定某基因至少要在多少个细胞中要检测到，低于设定值则丢弃；
min.features：指定某细胞至少有多少个基因表达，低于设定值则丢弃。
## 批量数据处理
scRNAlist <- list()
for(i in 1:length(dir_name)){
counts <- Read10X(data.dir = dir_name[i])
scRNAlist[[i]] <- CreateSeuratObject(counts, min.cells = 3, min.features =300)
}按照步骤4"dir_name"向量中的顺序，[[1]]中为BC21，[[2]]中为BC3。
6. 批量计算线粒体和红细胞比例
PercentageFeatureSet函数格式：PercentageFeatureSet(object,pattern = NULL,……)
object：Seurat对象；
pattern：用于匹配特征的正则表达式模式;
features：已定义的基因集，如果提供了pattern，将忽略该模式匹配。
for(i in 1:length(scRNAlist)){
sc <- scRNAlist[[i]]
# 计算线粒体比例
sc[["mt_percent"]] <- PercentageFeatureSet(sc, pattern = "^MT-")
# 计算红细胞比例
HB_genes <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ")
HB_m <- match(HB_genes, rownames(sc@assays$RNA))
HB_genes <- rownames(sc@assays$RNA)[HB_m]
HB_genes <- HB_genes[!is.na(HB_genes)]
sc[["HB_percent"]] <- PercentageFeatureSet(sc, features=HB_genes)
# 将sc赋值给scRNAlist[[i]]
scRNAlist[[i]] <- sc
# 删除sc
rm(sc)
}
"^MT-":表示匹配所有以"MT-"开头的基因；在gene symbol中人线粒体的基因格式全部都以"MT-"开头（小鼠为"mt-"）
6.1 查看计算后的线粒体和红细胞比例信息
以[[1]]BC21为例，计算后的线粒体和红细胞数据储存在下图红框"meta.data"中。
查看BC21的"meta.data"信息：
scRNAlist[[1]]@meta.data[1:5,1:5]
# orig.ident nCount_RNA nFeature_RNA mt_percent HB_percent
#BC21_AAACCCAAGAGATGCC-1 BC21 2570 954 5.019455 0.00000000
#BC21_AAACCCAAGCCTCTCT-1 BC21 1215 691 6.913580 0.00000000
#BC21_AAACCCAAGCGCTTCG-1 BC21 10485 3509 7.200763 0.00000000
#BC21_AAACCCAGTCACTCGG-1 BC21 15725 3972 3.821940 0.00635930
#BC21_AAACCCAGTGGCTAGA-1 BC21 7569 1616 16.633637 0.01321178
6.2 批量绘制过滤前小提琴图
# 批量绘制小提琴图
violin_before <- list()
for(i in 1:length(scRNAlist)){
violin_before[[i]] <- VlnPlot(scRNAlist[[i]],
features = c("nFeature_RNA", "nCount_RNA", "mt_percent","HB_percent"),
pt.size = 0.01,
ncol = 4)
}
# 合并图片
violin_before_merge <- CombinePlots(plots = violin_before,nrow=length(scRNAlist),legend='none')
# 将图片输出到画板上
violin_before_merge
# 保存图片
ggsave("violin_before_merge.pdf", plot = violin_before_merge, width = 15, height =7)
图片展示
nFeature_RNA：每个细胞检测表达的基因数目；
nCount_RNA：每个细胞测序的UMI count的表达量（即：每个细胞中基因的表达量进行相加，相加结果即为nCount）；
mt_percent：每个细胞线粒体基因表达量占总体基因的比例；
HB_percent：每个细胞红细胞基因表达量占总体基因的比例。
6.3 合并过滤前的样本绘制小提琴图
## 合并数据框
scm1 <- merge(scRNAlist[[1]],y=scRNAlist[[2:length(scRNAlist)]])
## 查看频数
table(scm$orig.ident)
#BC21 BC3
#6342 8684
## 绘图
violin_before_scm <- VlnPlot(scm1,
features = c("nFeature_RNA", "nCount_RNA", "mt_percent","HB_percent"),
pt.size = 0.01,
ncol = 4)
## 将图片输出到画板
violin_before_scm
# 保存图片
ggsave("violin_before_scm.pdf", plot = violin_before_scm, width = 15, height =7)
图片展示
7. 批量过滤细胞
一般默认线粒体含量至少要小于20%，红细胞的数目要至少小于5%；
在这里我们将过滤严格一点，调整为：
nFeature_RNA：每个细胞检测表达的基因数目大于300，小于7000；
nCount_RNA:每个细胞测序的UMI count含量大于1000，且剔除最大的前3%的细胞；
mt_percent:每个细胞的线粒体基因表达量占总体基因的比例小于10%；
HB_percent：每个细胞红细胞基因表达量占总体基因的比例小于3%。
scRNAlist <- lapply(X = scRNAlist, FUN = function(x){
x <- subset(x,
subset = nFeature_RNA > 300 & nFeature_RNA < 7000 &
mt_percent < 10 &
HB_percent < 3 &
nCount_RNA < quantile(nCount_RNA,0.97) &
nCount_RNA > 1000)})
没有固定的阈值标准，要根据自己的数据调整参数不断尝试，才能找到最佳结果。
7.1 批量绘制过滤后小提琴图
# 绘图
violin_after <- list()
for(i in 1:length(scRNAlist)){
violin_after[[i]] <- VlnPlot(scRNAlist[[i]],
features = c("nFeature_RNA", "nCount_RNA", "mt_percent","HB_percent"),
pt.size = 0.01,
ncol = 4)
}
# 合并图片
violin_after_merge <- CombinePlots(plots = violin_before,nrow=length(scRNAlist),legend='none')
# 将图片输出到画板上
violin_after_merge
# 保存图片
ggsave("violin_after_merge.pdf", plot = violin_after_merge, width = 15, height =7)
图片展示
7.2 合并过滤后的样本绘制小提琴图
## 合并数据框
scm2 <- merge(scRNAlist[[1]],scRNAlist[[2:length(scRNAlist)]])
## 查看频数
table(scm2$orig.ident)
#BC21 BC3
#3854 4250
## 绘图
violin_after_scm <- VlnPlot(scm2,
features = c("nFeature_RNA", "nCount_RNA", "mt_percent","HB_percent"),
pt.size = 0.01,
ncol = 4)
## 将图输出到画板
violin_after_scm
## 保存图片
ggsave("violin_after_scm.pdf", plot = violin_after_scm, width = 15, height =7)
图片展示
7.3 合并质控前后图片
# 合并
scplot_merge <- CombinePlots(list(violin_before_scm,violin_after_scm),nrow=2,legend="none")
# 将图输出到画板
scplot_merge
# 保存图片
ggsave("scplot_merge.pdf", plot = scplot_merge, width = 15, height =7)
图片展示
8. 批量数据归一化与筛选高变基因
NormalizeData函数格式：NormalizeData(object,normalization.method = "LogNormalize",scale.factor = 10000,……)
object：过滤后的Seurat对象；
normalization.method：归一化的方法(LogNormalize、CLR、RC)；
scale.factor：设置细胞归一化的比例因子。
整段意思为：对每个细胞的每个基因的表达量除以总表达量，然后乘以比例因子10000（不乘以10000取Log后数据小数点会很多，不好看），然后进行log归一化（LogNormalize目的是让整体的数据服从正态分布）。
FindVariableFeatures函数格式：FindVariableFeatures(object,selection.method = "vst",nfeatures = 2000,……)
object：归一化后的Seurat对象；
selection.method：高亮基因选择方法（vst、mvp、disp），推荐vst更主流；
nfeatures：设置高亮基因数量。
for (i in 1:length(scRNAlist)) {
scRNAlist[[i]] <- NormalizeData(scRNAlist[[i]],normalization.method = "LogNormalize", scale.factor = 10000)
scRNAlist[[i]] <- FindVariableFeatures(scRNAlist[[i]], selection.method = "vst",nfeatures = 3000)
}
官方推荐一般选择2000个高变基因，很多文章也有设置30000的，这个因自己的实验项目决定。
8.1 查看归一化数据与高变基因
以[[1]]BC21为例，归一化后的数据存储在下图红框"data"中，高变基因储存在"var.features"中。
9. 锚定法进行多样本整合
9.1 原理
① CCA降维（图A/B）
CCA将数据降维到同一个低维空间，表达相似的细胞特定的簇聚在一起；
② MNN算法获取"锚点细胞"（图C）
MNN（相互最近邻）算法计算得到两个数据集之间互相"距离"最近的细胞，称之为"锚点细胞"；
③ 过滤不正确的锚点（图D）
一般相同类型和状态的细胞才能构成配对锚点细胞（图C灰色线条），但是图C中"Query"黑色细胞团在"reference"中没有相同类型的细胞却也找到了锚点配对细胞（红色线条），需要将这些不正确的锚点过滤掉；
④ 样本整合（图E）
计算差异向量，用此向量校正这个锚点锚定的细胞子集的基因表达值。校正后的基因表达值即消除了技术偏倚，实现了两个单细胞数据集的整合。
9.2 样本整合与降维聚类
① 查找与过滤错误锚点
FindIntegrationAnchors函数格式：FindIntegrationAnchors(object.list,anchor.features=2000,……)
object.list：Seurat对象；
anchor.features：指定锚定点的基因数。
scRNA_anchors <- FindIntegrationAnchors(object.list = scRNAlist,anchor.features = 2000)
② 整合数据
scRNA1 <- IntegrateData(anchorset = scRNA_anchors)
查看整合后信息
整合后我们可以看到"scRNA1"中"assays"生成了两个层级结构，一个是"RNA"（20513个基因），一个是"integrated"（2000个锚定点/基因）。
同时需要注意的是"active.assay"显示的是"integrated"，也就是现在默认的assay是在"integrated"里面。
Tips：锚定点是用于整合数据的不算用于后续作图，如果后续作图需要将"integrated"改变为RNA。
③ 数据降维聚类
# 将默认"scRNA1"设置成"integrated"
DefaultAssay(scRNA1) <- "integrated"
# "integrated"数据标准化
scRNA1=ScaleData(scRNA1)
# PCA降维：其他降维的基础
scRNA1 <- RunPCA(scRNA1, npcs = 30, verbose = T)
# 聚类
scRNA1 <- FindNeighbors(scRNA1, reduction = "pca", dims = 1:15)
scRNA1 <- FindClusters(scRNA1, resolution = 1)
# umap/tsne降维
scRNA1 <- RunUMAP(scRNA1,dims = 1:15)
scRNA1 <- RunTSNE(scRNA1,dims = 1:15)
# 绘图
umap_integrated1 <- DimPlot(scRNA1, reduction = "umap", group.by = "orig.ident")
umap_integrated2 <- DimPlot(scRNA1, reduction = "umap", label = TRUE)
tsne_integrated1 <- DimPlot(scRNA1, reduction = "tsne", group.by = "orig.ident")
tsne_integrated2 <- DimPlot(scRNA1, reduction = "tsne", label = TRUE)
# 合并图片
umap_tsne_integrated <- CombinePlots(list(tsne_integrated1,tsne_integrated2,umap_integrated1,umap_integrated2),ncol=2)
# 将图片输出到画板
umap_tsne_integrated
# 保存图片
ggsave("umap_tsne_integrated.pdf",umap_tsne_integrated,wi=25,he=15)
图片展示
10. 对"RNA"数据标准化
锚定整合过程中一定需要做两次数据标准化，第一次是对"integrated"标准化，第二次是对"RNA"标准化。
# 将默认"scRNA1"设置成"RNA"
DefaultAssay(scRNA1) <- "RNA"
# 数据标准化
scRNA2 <- ScaleData(scRNA1)
11. 保存数据
save(scRNA1,scRNA2,scRNAlist,file = "scdata1.Rdata")腾讯云开发者
