CCA 和 Harmony在整合pbmc3k和pbmc5k的效果比较
挖掘到了一个段子手学徒，总是给我惊喜。把枯燥无味的知识点讲解的让人捧腹大笑！
单细胞的多组对照设计（例如正常组与给药组）可以为细胞类型水平比较提供以往Bulk RNA-seq分析所不能达到的精度。对此一般有两种进阶分析思路：
（1）DE(Differential expression)--两组样本的同一细胞类型的基因表达差异分析； （2）DA(Differential abundance)--两组样本的同一细胞类型的丰度差异分析
参考：教程第14章 http://bioconductor.org/books/release/OSCA/overview.html
这样的多分组样品往往是需要整合， 目前比较流行的是seurat包自带的CCA 和 Harmony两个方法， 所以安排学徒在整合pbmc3k和pbmc5k的效果比较。
前情提要
“随着分析流程的逐渐完善，上游分析的套路也愈发的固定，各方用户表示用的安心，回车键按的放心！”
以上是学徒办公室前线发来的报到！
但是根据曾老师的指导，情况在这时发生了新的变化：
“尊重不同数据的实际情况，选择合适的、能解决实际问题的整合算法！”
以下是详细报道：
最近的单细胞数据练习中，多次出现CCA拉胯的现象，比如在数据集 GSE149165 的分析练习中，出现了一个奇怪的分叉，分组差异这么明显吗？这合理吗？
然后经过对原文的深入理解，找到了问题的关键，文章明确说到了Harmony，并且参照原文的分析结果，组间差异没有这么明显，然后又重新做了一下，这个分叉就没了，这样的差异在分群的时候都没有被单独分为一个亚群，对后续的注释以及深入理解分析结果的影响不会超过细胞类型之间的差异，勉强属于可接受范围，但是新的风暴已经来临—— GSE162115
这个是已经注释的结果，可以看到 CAF 亚群的数量在 G1/G2 之间的差异十分明显，不同于前一个数据集明确要求使用Harmony整合方式，该数据集的文章并没有明确，而且由于该文章对于一些模糊分群的骚操作（删除Mark-genes不显著的分群），所以也无法对照原文结果来确定这个差异到底是组间差异，还是整合算法不合适。
CCA和Harmony之间的差异到底相差多少？
两者之间的差距是否会影响分析结果？
黑猫白猫
为了减少文献中来自预先实验设计造成的批次影响，这个测试我们选择10X官方提供的，PBMC的3k和5k单细胞数据，其相关的分析教程和结果在互联网中随手可得，话不多说，这就开始：
1、众生平等的过滤
过滤环节基本就是千年王八的屁股——老规定了，回车键敲起来！
#########################################################################################
####################### 文件读取 #####################################################
# 这个 input 文件夹存放了 pbmc3k和pbmc5k 两个文件夹，各种里面都是10x的3个文件哦
samples=list.files('input/')
samples
sceList = lapply(samples,function(pro){
#pro=samples[1]
folder=file.path('input/',pro)
print(pro)
print(folder)
print(list.files(folder))
sce=CreateSeuratObject(counts = Read10X(folder),
project = pro )
return(sce)
})
names(sceList)
names(sceList) = samples
sceList
sce <- merge(sceList[[1]], y= sceList[ -1 ] ,
add.cell.ids=samples)
sce.all <- sce
########################################################################################
####################### 数据质控 ####################################################
dir.create("./1-QC")
setwd("./1-QC")
#计算线粒体基因比例v
mito_genes=rownames(sce.all)[grep("^MT-", rownames(sce.all))]
mito_genes #13个线粒体基因
sce.all=PercentageFeatureSet(sce.all, "^MT-", col.name = "percent_mito")
#计算核糖体基因比例
ribo_genes=rownames(sce.all)[grep("^Rp[sl]", rownames(sce.all),ignore.case = T)]
ribo_genes
sce.all=PercentageFeatureSet(sce.all, "^RP[SL]", col.name = "percent_ribo")
fivenum(sce.all@meta.data$percent_ribo)
#计算红血细胞基因比例
rownames(sce.all)[grep("^HB[^(p)]", rownames(sce.all))]
sce.all=PercentageFeatureSet(sce.all, "^HB[^(p)]", col.name = "percent_hb")
#可视化细胞的上述比例情况
feats <- c("nFeature_RNA", "nCount_RNA", "percent_mito", "percent_ribo", "percent_hb")
feats <- c("nFeature_RNA", "nCount_RNA")
p1=VlnPlot(sce.all, group.by = "orig.ident", features = feats, pt.size = 0.01, ncol = 2) +
NoLegend()
p1
library(ggplot2)
ggsave(filename="Vlnplot1.png",plot=p1)
feats <- c("percent_mito", "percent_ribo", "percent_hb")
p2=VlnPlot(sce.all, group.by = "orig.ident", features = feats, pt.size = 0.01, ncol = 3, same.y.lims=T) +
scale_y_continuous(breaks=seq(0, 100, 5)) +
NoLegend()
p2
ggsave(filename="Vlnplot2.png",plot=p2)
p3=FeatureScatter(sce.all, "nCount_RNA", "nFeature_RNA", group.by = "orig.ident", pt.size = 0.5)
p3
ggsave(filename="Scatterplot.png",plot=p3)
#根据上述指标，过滤低质量细胞/基因
#过滤指标1:最少表达基因数的细胞&最少表达细胞数的基因
selected_c <- WhichCells(sce.all, expression = nFeature_RNA > 600)
selected_f <- rownames(sce.all)[Matrix::rowSums(sce.all@assays$RNA@counts > 0 ) > 3]
sce.all.filt <- subset(sce.all, features = selected_f, cells = selected_c)
dim(sce.all)
dim(sce.all.filt)
table(sce.all@meta.data$orig.ident)
table(sce.all.filt@meta.data$orig.ident)
# 可以看到，主要是过滤了基因，其次才是细胞
C=sce.all.filt@assays$RNA@counts
dim(C)
C=Matrix::t(Matrix::t(C)/Matrix::colSums(C)) * 100
most_expressed <- order(apply(C, 1, median), decreasing = T)[50:1]
pdf("TOP50_most_expressed_gene.pdf",width=14)
boxplot(as.matrix(Matrix::t(C[most_expressed, ])),
cex = 0.1, las = 1,
xlab = "% total count per cell",
col = (scales::hue_pal())(50)[50:1],
horizontal = TRUE)
dev.off()
rm(C)
#过滤指标2:线粒体/核糖体基因比例(根据上面的violin图)
selected_mito <- WhichCells(sce.all.filt, expression = percent_mito < 15)
selected_ribo <- WhichCells(sce.all.filt, expression = percent_ribo > 3)
selected_hb <- WhichCells(sce.all.filt, expression = percent_hb < 0.1)
length(selected_hb)
length(selected_ribo)
length(selected_mito)
sce.all.filt <- subset(sce.all.filt, cells = selected_mito)
sce.all.filt <- subset(sce.all.filt, cells = selected_ribo)
sce.all.filt <- subset(sce.all.filt, cells = selected_hb)
dim(sce.all.filt)
table(sce.all.filt$orig.ident)
#可视化过滤后的情况
feats <- c("nFeature_RNA", "nCount_RNA")
p1_filtered=VlnPlot(sce.all.filt, group.by = "orig.ident", features = feats, pt.size = 0.1, ncol = 2) +
NoLegend()
ggsave(filename="Vlnplot1_filtered.png",plot=p1_filtered,width = 8)
feats <- c("percent_mito", "percent_ribo", "percent_hb")
p2_filtered=VlnPlot(sce.all.filt, group.by = "orig.ident", features = feats, pt.size = 0.1, ncol = 3) +
NoLegend()
ggsave(filename="Vlnplot2_filtered.png",plot=p2_filtered, width = 8)
#细胞周期评分
sce.all.filt = NormalizeData(sce.all.filt)
s.genes=Seurat::cc.genes.updated.2019$s.genes
g2m.genes=Seurat::cc.genes.updated.2019$g2m.genes
sce.all.filt=CellCycleScoring(object = sce.all.filt,
s.features = s.genes,
g2m.features = g2m.genes,
set.ident = TRUE)
p4=VlnPlot(sce.all.filt, features = c("S.Score", "G2M.Score"), group.by = "orig.ident",
ncol = 2, pt.size = 0.1)
ggsave(filename="Vlnplot4_cycle.png",plot=p4)
sce.all.filt@meta.data %>% ggplot(aes(S.Score,G2M.Score))+geom_point(aes(color=Phase))+
theme_minimal()
ggsave(filename="cycle_details.png" )
# S.Score较高的为S期，G2M.Score较高的为G2M期，都比较低的为G1期
#检测doublets
sce.all.filt = FindVariableFeatures(sce.all.filt)
sce.all.filt = ScaleData(sce.all.filt,
vars.to.regress = c("nFeature_RNA", "percent_mito"))
sce.all.filt = RunPCA(sce.all.filt, npcs = 20)
sce.all.filt = RunTSNE(sce.all.filt, npcs = 20)
sce.all.filt = RunUMAP(sce.all.filt, dims = 1:10)
sce.all.filt
sce.all.filt <- FindNeighbors(sce.all.filt, dims = 1:15)
sce.all.filt <- FindClusters(sce.all.filt, resolution = 0.8)
table(sce.all.filt@meta.data$RNA_snn_res.0.8)
save(sce.all.filt,file = 'sce.all.filt.Rdata')
setwd('../')
2、CCA
开始抓耗子！CCA也是比较常见的环节了：
load(file = '../../1-QC/sce.all.filt.Rdata')
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
sce.all=sce.all.filt
#拆分为 个seurat子对象
sce.all.list <- SplitObject(sce.all, split.by = "orig.ident")
##执行CCA算法整合3k和5k的数据
for (i in 1:length(sce.all.list)) {
print(i)
sce.all.list[[i]] <- NormalizeData(sce.all.list[[i]], verbose = FALSE)
sce.all.list[[i]] <- FindVariableFeatures(sce.all.list[[i]], selection.method = "vst",
nfeatures = 2000, verbose = FALSE)
}
alldata.anchors <- FindIntegrationAnchors(object.list = sce.all.list, dims = 1:30,
reduction = "cca")
sce.all.int <- IntegrateData(anchorset = alldata.anchors, dims = 1:30, new.assay.name = "CCA")
##执行将为算法
sce.all.int=ScaleData(sce.all.int)
sce.all.int=RunPCA(sce.all.int, npcs = 30)
sce.all.int=RunTSNE(sce.all.int, dims = 1:30)
sce.all.int=RunUMAP(sce.all.int, dims = 1:30)
names(sce.all.int@reductions)
names(sce.all@reductions)
colnames(sce.all@meta.data)
#比较整合前后的降维结果
p1.compare=plot_grid(ncol = 3,
DimPlot(sce.all, reduction = "pca", group.by = "orig.ident")+NoAxes()+ggtitle("PCA raw_data"),
DimPlot(sce.all, reduction = "tsne", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE raw_data"),
DimPlot(sce.all, reduction = "umap", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP raw_data"),
DimPlot(sce.all.int, reduction = "pca", group.by = "orig.ident")+NoAxes()+ggtitle("PCA integrated"),
DimPlot(sce.all.int, reduction = "tsne", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE integrated"),
DimPlot(sce.all.int, reduction = "umap", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP integrated")
)
p1.compare
ggsave(plot=p1.compare,filename="Before&After_int.pdf", width = 15, height = 8)
###### step4:分cluster######
sce.all=FindNeighbors(sce.all, dims = 1:30, k.param = 60, prune.SNN = 1/15)
#设置不同的分辨率，观察分群效果(选择哪一个？)
for (res in c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1,1.2,1.5,2)) {
sce.all=FindClusters(sce.all, graph.name = "CCA_snn", resolution = res, algorithm = 1)}
apply(sce.all@meta.data[,grep("CCA_snn_res",colnames(sce.all@meta.data))],2,table)
p2_tree=clustree(sce.all@meta.data, prefix = "CCA_snn_res.")
ggsave(plot=p2_tree, filename="Tree_diff_resolution.pdf", width = 8, height = 10)
#接下来分析，按照分辨率为0.8进行
sel.clust = "CCA_snn_res.0.8"
sce.all <- SetIdent(sce.all, value = sel.clust)
table(sce.all@active.ident)
saveRDS(sce.all, "sce.all_int.rds")
这个代码段主要是执行了CCA的整合算法和降维算法，这里有两个比较关键的地方，就是一个整合前后的比较图：
从这里可以看出整合环节是十分必要的，而且效果不错。此外就是针对不同的 resolution 阈值进行的分群聚类树，一般都需要根据这个结果来选择合适的分群阈值，并且这个结果会受到数据本身情况的直接影响：
这个图可以看出，三个相对独立的分群之间井水不犯河水，这种情况就紧着高的用吧，本次测试的阈值是 0.8。
然后就是注释：
th=theme(axis.text.x = element_text(angle = 45,
vjust = 0.5, hjust=0.5))
sce.all <- SetIdent(sce.all, value = sce.all$CCA_snn_res.0.8)
#这里的Mark比一般分析的少，主要是已知细胞类型不那么多
DotPlot(sce.all, features = c("MS4A1", "GNLY", "CD3E",
"CD14", "FCER1A", "FCGR3A",
'C1QA', 'C1QB', # mac
'S100A9', 'S100A8', 'MMP19',# monocyte
'LAMP3', 'IDO1','IDO2',## DC3
'CD1E','CD1C', # DC2
"LYZ", "PPBP", "CD8A"), assay = "RNA") + coord_flip() +th
ggsave("anno_marker_gene_check.png", width = 10, height = 8)
DimPlot(sce.all, reduction = "umap", pt.size = 1.0, label = T)+NoLegend()
ggsave(filename = "UMAP_res_0.8.png")
# 需要自行看图，定细胞亚群：
celltype=data.frame(ClusterID=0:10,
celltype='Tcells')
celltype[celltype$ClusterID %in% c( 4,5 ),2]='Bcells'
#celltype[celltype$ClusterID %in% c( 1,8,10),2]='myeloid'
celltype[celltype$ClusterID %in% c( 1 ),2]='mono'
celltype[celltype$ClusterID %in% c( 9 ),2]='mac'
celltype[celltype$ClusterID %in% c( 10 ),2]='DC'
#把注释结果导入seurat对象
head(celltype)
celltype
table(celltype$celltype)
sce.all@meta.data$celltype = "NA"
for(i in 1:nrow(celltype)){
sce.all@meta.data[which(sce.all@meta.data$CCA_snn_res.0.8 == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}
table(sce.all@meta.data$celltype)
看一下最后的比对结果，可以说是岁月静好：
3、Harmony
为什么总觉得这个单词眼熟呢？（小声BB）
然后是在最近两次分析练习中的黑马选手上台，它的车很快，是基于主成分分析，而且用的是迭代矫正过弯（围观自：科学网—“harmony”整合不同平台的单细胞数据之旅 - 陈同的博文 (sciencenet.cn)），我没有看清车牌，只看到发表在 Cell 上，如果你们有人认识它，麻烦跟它说一下，今天晚上，CCA在秋名山等它去批次效应！
load(file = '../../1-QC/sce.all.filt.Rdata')
table(sce.all.filt$orig.ident)
library(harmony)
#添加一个肘图参数，可以更好的确定pca阈值
sce.all.int <- RunHarmony(sce.all.filt,c( "orig.ident" ), plot_convergence = T)
ggsave(filename = "Harmony_plot_convergence.png", width = 10, height = 8)
harmony_embeddings <- Embeddings(sce.all.int, 'harmony')
#注意一下，这里执行降维算法的时候，是基于harmony，不指定的话，默认是PCA
sce.all.int=RunTSNE(sce.all.int,reduction = "harmony", dims = 1:20)
sce.all.int=RunUMAP(sce.all.int,reduction = "harmony",dims = 1:20)
names(sce.all.int@reductions)
names(sce.all@reductions)
DimPlot(sce.all.int , reduction = "umap",group.by ='orig.ident' )
ggsave('umap_by_orig.ident_after_harmony.png')
#比对整合结果
sce.all=sce.all.filt
p1.compare=plot_grid(ncol = 3,
DimPlot(sce.all, reduction = "pca", group.by = "orig.ident")+NoAxes()+ggtitle("PCA raw_data"),
DimPlot(sce.all, reduction = "tsne", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE raw_data"),
DimPlot(sce.all, reduction = "umap", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP raw_data"),
DimPlot(sce.all.int, reduction = "pca", group.by = "orig.ident")+NoAxes()+ggtitle("PCA integrated"),
DimPlot(sce.all.int, reduction = "tsne", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE integrated"),
DimPlot(sce.all.int, reduction = "umap", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP integrated")
)
p1.compare
这里也先看一下与原始数据的比较结果：
没毛病老铁！Harmony还有个好处就是，可以在降维过程中直接绘制一个肘图来确定后续分析阈值，只要在RunHarmony()里加一个plot_convergence = T参数就好了，这里可以看一下这次分析的结果，这个肘图相对PCA的结果来说，断层很明显，不存在选择困难
看来只能进一步的深入分析来决胜负了，注释环节走一个：
sce.all=sce.all.int
#这里也是要指定harmony的结果的，还要指定dims的阈值
sce.all=FindNeighbors(sce.all, reduction = "harmony",
dims = 1:20 )
colnames(sce.all@meta.data)
sce.all@graphs$RNA_snn
#设置不同的分辨率，观察分群效果(选择哪一个？)
for (res in c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1)) {
sce.all=FindClusters(sce.all,
resolution = res,
algorithm = 1)}
apply(sce.all@meta.data[,grep("RNA_snn_res",
colnames(sce.all@meta.data))],2,table)
p2_tree=clustree(sce.all@meta.data, prefix = "RNA_snn_res.")
ggsave(plot=p2_tree, filename="Tree_diff_resolution.pdf",width = 7, height = 10)
p2_tree
table(sce.all@meta.data$RNA_snn_res.0.3)
在注释之前，差异从这个代码段开始也逐渐拉大，主要的差异是resolution的阈值上：
可以看到从 res=0.3 开始，亚群间的细胞交换逐渐混乱，所以后续的分析选 0.3。相比CCA的整合结果，较为混乱一些，根据这样的结果来注释看看
th=theme(axis.text.x = element_text(angle = 45,
vjust = 0.5, hjust=0.5))
DotPlot(sce.all, features = c("MS4A1", "GNLY", "CD3E",
"CD14", "FCER1A", "FCGR3A",
'C1QA', 'C1QB', # mac
'S100A9', 'S100A8', 'MMP19',# monocyte
'LAMP3', 'IDO1','IDO2',## DC3
'CD1E','CD1C', # DC2
"LYZ", "PPBP", "CD8A",
'CD203c','CD63', 'CD123',"CD68"), assay = "RNA") + coord_flip() +th
ggsave("anno_marker_gene_check.png", width = 10, height = 8)
DimPlot(sce.all, reduction = "umap", pt.size = 1.0, label = T)+NoLegend()
ggsave(filename = "UMAP_res_0.3.png")
# "DC", "Platelet"
# 需要自行看图，定细胞亚群：
celltype=data.frame(ClusterID=0:9,
celltype='Tcells')
celltype[celltype$ClusterID %in% c( 4,5,9 ),2]='Bcells'
celltype[celltype$ClusterID %in% c( 2 ),2]='monocyte'
celltype[celltype$ClusterID %in% c( 6 ),2]='mac'
celltype[celltype$ClusterID %in% c( 7 ),2]='DC'
celltype[celltype$ClusterID %in% c( 8 ),2]='Platelet'
head(celltype)
celltype
table(celltype$celltype)
sce.all@meta.data$celltype = "NA"
for(i in 1:nrow(celltype)){
sce.all@meta.data[which(sce.all@meta.data$RNA_snn_res.0.3 == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}
table(sce.all@meta.data$celltype)
这里的根据相同的 Mark-genes 进行了注释，并且 res=0.3 所得到的分群数量（10）和之前CCA的0.8差的不多（n=11）所以相对的结果也差不多，但是这里多出来了一个血小板细胞的小分群
抛砖引玉
黑猫白猫抓到耗子就是好猫，那么CCA和Harmony两者明显都抓到耗子了，这个耗子大不大、毛色如何等标准实际还是有得比，这里的分析结果其实还是局限在了抓没抓住的问题上。
从左边的 Mark-genes 的气泡图里甚至分不出差别来，而且注释的结果也是差不多的，那就抽丝剥茧，详细的看看这个结果的差异：
1、数量
单细胞分析，既然已经是单细胞的维度，那么细胞数量本身就已经是一个关键指标了，所以先把各个细胞类型的细胞数量比较一下：
这个图直观的展示了两者的注释结果之间的差异，可以看到还是存在一些离散的数值，但都是个位数了
2、表达量
通过AverageExpression()获取对应的细胞类型的表达量，然后绘制表达量热 图。可以看到两次分析之间各细胞类型都是强相关性（x轴为CCA，y轴为Harmony）。
3、总结
两相比较之下，大概可以得出以下结论：
(1) CCA 和 Harmony 对于批次效应，都有相差不大的整合效力
(2) Harmony 对于 resolution 的取值范围相对 CCA 存在差异
(3) Harmony 对于小分群较为敏感
实际上单细胞分析中的批次效应来源是多样化的，在复杂的实验方案中，不仅仅是测序批次，还有组间、分析平台间等多种批次效应，而在各种情况，不同的整合算法是否各有优势，又或者相差不大，这些都是需要深入讨论的。此外，对于生物学差异和批次效应的薛定谔效应也需要深厚的背景知识来处理。本文只是抛砖引玉的对目前遇到的两种算法进行了比较，希望随着学习的深入能够更新这个系列，或者有大佬能来个降维打击，一步到位~[手动狗头]
预知后事如何，且听下回分解！