单细胞转录组 | 多样本处理与Harmony整合
前言
上期推文单细胞转录组 | 多样本处理与锚定法整合介绍了使用锚定法进行多个样本整合，本期我们来介绍另一个多样本整合的主流方法：Harmony。
本文框架
1. Rtools安装
使用harmony需要安装Rtools，如果不安装后续分析会报错。若已经安装此步请跳过。
① 打开网站选择红色箭头指向的"RTools 4.0"，进入详情页；
https://cran.r-project.org/bin/windows/Rtools/
② 选择红色箭头指向的"rtools40-x86_64.exe"下载；
③ 下载安装完成后（安装路径随意），进入Rstudio，在控制台输入；
write('PATH="${RTOOLS40_HOME}\\usr\\bin;${PATH}"', file = "~/.Renviron", append = TRUE)
④ 如果一切正常没有报错，重启Rstudio；
⑤ 测试路径配置是否成功。
Sys.which("make")
只要输出不是空字符串就表明路径配置成功。
可能出现的报错：
In file(file, ifelse(append, "a", "w")) : cannot open file 'D:/????/.Renviron': Invalid argument
解决办法：
① 重启Rstudio后运行getwd()命令，获取此时的工作目录。在工作目录中创建txt文件，将 PATH="
② 再次重启Rstudio；
③ 输入"Sys.which("make")"测试路径配置是否成功。
Sys.which("make")
只要输出不是空字符串就表明路径配置成功。
2. 安装包
如果已经安装，此步请跳过。
install.packages('Seurat')
install.packages('dplyr')
install.packages('tidyverse')
install.packages('patchwork')
install.packages("devtools")
install_github("immunogenomics/harmony")
3. 加载包
library(Seurat)
library(dplyr)
library(tidyverse)
library(patchwork)
library(devtools)
library(harmony)
4. 设置工作路径
setwd("D:/sc-seq/")
请根据自己数据的存放位置自定义路径。
本次示例工作路径下存放了需要读取的10×数据文件夹：BC3和BC21。
5. 创建读取文件的向量
## 将文件夹中的文件名存到"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"
6. 批量读取数据并创建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。
7. 批量计算线粒体和红细胞比例
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-"）。
以[[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
8. 批量过滤细胞
一般默认线粒体含量至少要小于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)})
没有固定的阈值标准，要根据自己的数据调整参数不断尝试，才能找到最佳结果。
9. 合并Seurat对象
## 合并
scRNAlist <- merge(scRNAlist[[1]],y=scRNAlist[[2]])
## 统计细胞数
table(scRNAlist[[]]$orig.ident)
#BC21 BC3
#3854 4250
合并BC21与BC3后的scRNAlist：
10. 筛选高变基因与降维
harmony整合是基于PCA降维结果进行的。
scRNAlist <- NormalizeData(scRNAlist) %>%
FindVariableFeatures(selection.method = "vst",nfeatures = 3000) %>%
ScaleData() %>%
RunPCA(npcs = 30, verbose = T)
归一化后的数据存储在下图红框"data"中，高变基因储存在黄框"var.features"中，PCA降维后的数据储存在蓝框pca中。
11. harmony整合
11.1 原理
Step1：Harmony概率性地将细胞分配给cluster，从而使每个cluster内数据集的多样性最大化;
Step2：Harmony计算每个cluster的所有数据集的全局中心，以及特定数据集的中心;
Step3：在每个cluster中，Harmony基于中心为每个数据集计算校正因子;
Step4：Harmony使用基于Step3的特定于细胞的因子校正每个细胞。由于Harmony使用软聚类，因此可以通过多个因子的线性组合对其A中进行的软聚类分配进行线性校正，来修正每个单细胞。重复步骤A到D，直到收敛为 止。聚类分配和数据集之间的依赖性随着每一轮的减少而减小。
11.2 整合
整合需要指定Seurat对象和metadata中需要整合的变量名。
scRNA_harmony <- RunHarmony(scRNAlist, group.by.vars = "orig.ident")
查看Harmony矫正之后的信息：
scRNA_harmony@reductions[["harmony"]][[1:5,1:5]]
# harmony_1 harmony_2 harmony_3 harmony_4 harmony_5
#BC21_AAACCCAAGAGATGCC-1 3.135328 5.669093 2.477945 1.7110049 9.2685516
#BC21_AAACCCAAGCCTCTCT-1 0.572945 6.311063 2.589623 -0.4654396 5.9814832
#BC21_AAACCCAAGCGCTTCG-1 -32.778030 -3.843416 -10.457380 4.6612394 0.8715071
#BC21_AAACCCAGTCACTCGG-1 9.638793 -10.191707 -6.170584 -4.5852456 1.2507496
#BC21_AAACCCAGTTCCATTT-1 -9.391769 3.392380 -3.588324 -1.8205631 0.2678927
11.3 聚类
后续都是基于Harmony矫正之后的数据，不是基因表达数据和直接的PCA降维数据。
设置reduction = 'harmony'，后续分析是基于Harmony矫正之后的数据。
# 聚类
scRNA_harmony <- FindNeighbors(scRNA_harmony, reduction = "harmony", dims = 1:15) %>% FindClusters(resolution = 1)
# umap/tsne降维
scRNA_harmony <- RunTSNE(scRNA_harmony, reduction = "harmony", dims = 1:15)
scRNA_harmony <- RunUMAP(scRNA_harmony, reduction = "harmony", dims = 1:15)
11.4 绘图
# 绘图
umap_integrated1 <- DimPlot(scRNA_harmony, reduction = "umap", group.by = "orig.ident")
umap_integrated2 <- DimPlot(scRNA_harmony, reduction = "umap", label = TRUE)
tsne_integrated1 <- DimPlot(scRNA_harmony, reduction = "tsne", group.by = "orig.ident")
tsne_integrated2 <- DimPlot(scRNA_harmony, 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)
图片展示
12. 保存数据
save(scRNA_harmony,scRNAlist,file = "scdata2.Rdata")腾讯云开发者
