单细胞降维聚类分群的另外一个工具选择Pagoda2

我五年前的单细胞数据处理视频公益课程一直都在b站无限制免费给大家学习，其实就是强调大家需要熟练掌握5个R包，分别是: scater,monocle,Seurat,scran,M3Drop，主要是 需要熟练掌握它们的对象，他们的分析流程也大同小异：

step1: 创建对象
step2: 质量控制
step3: 表达量的标准化和归一化
step4: 去除干扰因素(多个样本整合)
step5: 判断重要的基因
step6: 多种降维算法
step7: 可视化降维结果
step8: 多种聚类算法
step9: 聚类后找每个细胞亚群的标志基因
step10: 继续分类

我都是这样教导学生完成单细胞学习的，基础课程学完后需要完成作业：https://mp.weixin.qq.com/s/lpoHhZqi-_ASUaIfpnX96w

但是最近五年，五大R包的说法已经是不复存在了，目前大家认可度比较好的就剩下了 Seurat，它的单细胞降维聚类分群的一条龙服务大家都超级熟悉了。但是其它包也是可以类似的处理，而且不同包的对比可以加深大家对这个流程的认识，所以我才会强调大家需要熟练掌握5个R包，分别是: scater,monocle,Seurat,scran,M3Drop。

其实也是有很多其它后起之秀，比如现在我们要分享的单细胞降维聚类分群的另外一个工具选择Pagoda2，之所以会接触它，主要是我们跑RNA速率分析的时候，velocyto.R 官方文档指出来的降维聚类分群方法就是它。

而且我简单搜索了一下，发现它也确实有不少引用，但是主要是Smart-Seq2的单细胞转录组数据，比如；

• 2020年七月瑞典卡洛琳斯卡研究所/阿斯利康综合心脏代谢中心Urban Lendahl 和 Christer Betsholtz于Nature communications上共同发表名为《Single-cell analysis uncovers fibroblast heterogeneity and criteria for fibroblast and mural cell identification and discrimination》的单细胞分析文章
• 以及之前在单细胞天地公众号看到的：单细胞水平看骨髓瘤的细胞状态和基因调控 ，选取8个relapsed/refractory MM病人骨髓或血液中的骨髓瘤细胞以及CD45+ 免疫细胞，以及2个健康供体，使用SmartSeq2得到单细胞全长转录组，共 6,955 cells，利用PAGODA2进行分群。

大家可以主要参考：http://pklab.med.harvard.edu/peterk/nano2019.html 文档进行学习。

安装和测试数据的认识

因为pagoda2是成熟的R包，在CRAN可以直接下载，同时安装conos包，因为里面有测试数据。代码如下所示：

install.packages('pagoda2')
install.packages('conos')

如果这两个包你安装有困难，建议花三五天重新学一学R语言哦，还有一个在GitHub的包也需要安装：

remotes::install_github('kharchenkolab/conosPanel')

都是超级简单的安装代码，但是你大概率会报错。

测试数据

以大家熟知的pbmc3k数据集为例。大家先安装这个数据集对应的包，并且对它进行降维聚类分群，参考前面的例子：人人都能学会的单细胞聚类分群注释 ，而且每个亚群找高表达量基因，都存储为Rdata文件。

# 0.安装R包 ----
# devtools::install_github('caleblareau/BuenColors')
# utils::install.packages(pkgs = "ggstatsplot")
# InstallData("pbmc3k") 

library(SeuratData) #加载seurat数据集  
getOption('timeout')
options(timeout=10000)
#InstallData("pbmc3k")  
data("pbmc3k")  
sce <- pbmc3k.final   
library(Seurat)
table(Idents(sce))
p1=DimPlot(sce,label = T)

因为这个数据集已经进行了不同单细胞亚群的标记，所以我们很容易拿到单细胞表达量矩阵和细胞对应的表型信息。

质量控制方面，这个pagoda2只需要一个表达量矩阵即可，代码如下所示：

cm=sce@assays$RNA@counts
str(cm) 
par(mfrow=c(1,2), mar = c(3.5,3.5,2.0,0.5), mgp = c(2,0.65,0), cex = 1.0)
hist(log10(colSums(cm)+1),main='molecules per cell',col='cornsilk',xlab='log10(molecules per cell)')
hist(log10(rowSums(cm)+1),main='molecules per gene',col='cornsilk',xlab='log10(molecules per gene])')

hist(log10(rowSums(counts)+1),main='Molecules per gene',xlab='molecules (log10)',col='cornsilk')
abline(v=1,lty=2,col=2)

counts <- counts[rowSums(counts)>=10,]
dim(counts) 
rownames(counts) <- make.unique(rownames(counts))

经过了上面的操作，我们就熟悉了这个pagoda2包的表达量矩阵输入信息，可以构建一个对象啦：

r <- Pagoda2$new(counts,log.scale=TRUE, n.cores=2) 
> class(r)
[1] "Pagoda2" "R6" 

接下来，所有的分析都是基于这个对象，变量名字有点普通，就是，后续分析就是对 r进行各种各样的处理。

说句实话，这个对象还是超级大难以理解，主要是不习惯这样的R6架构。

降维聚类分群

差不多也是一条路，寻找并且选择高变基因，PCA降维，KNN聚类分群，代码如下所示：

# 一切的输入数据，都是 counts 这样的 纯粹的表达量矩阵哦  
r <- Pagoda2$new(counts,log.scale=TRUE, n.cores=2) 
r$adjustVariance(plot=T,gam.k=10) 
r$calculatePcaReduction(nPcs=50,n.odgenes=3e3) 
r$makeKnnGraph(k=40,type='PCA',center=T,distance='cosine') 
r$getKnnClusters(method=infomap.community,type='PCA')

M <- 30; 
r$getEmbedding(type='PCA',embeddingType = 'largeVis', 
               M=M,perplexity=30,gamma=1/M,alpha=1)
r$plotEmbedding(type='PCA',show.legend=F,
                shuffle.colors=F,
              alpha=0.1 ) + ggtitle('clusters (largeVis)')

最简单的降维聚类分群结果如下所示：

降维聚类分群结果

当然了，我们大概率并不会直接使用PCA这样的原始降维形式，会喜欢tSNE等等，这个pagoda2也是可以很简单的操作，代码如下所示：

r$getEmbedding(type='PCA',embeddingType='tSNE',perplexity=50,verbose=F,n.cores=30)
r$plotEmbedding(type='PCA',embeddingType='tSNE',show.legend=F,
                min.group.size=1,
                shuffle.colors=F,mark.cluster.cex=1,
                alpha=0.1,main='clusters (tSNE)')

下面是tSNE的效果：

tSNE的效果

一些可视化

我们前面提到了，使用seurat标准流程有五个标准图表，详见以前我们做的投票：可视化单细胞亚群的标记基因的5个方法，下面的5个基础函数相信大家都是已经烂熟于心了：

• VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
• FeaturePlot(pbmc, features = c("MS4A1", "CD79A"))
• RidgePlot(pbmc, features = c("MS4A1", "CD79A"), ncol = 1)
• DotPlot(pbmc, features = unique(features)) + RotatedAxis()
• DoHeatmap(subset(pbmc, downsample = 100), features = features, size = 3)

同样的，这个pagoda2流程也是可以的，首先是检查不同的基因表达量，类似于FeaturePlot，代码：

gene <-"HBB" 
r$plotEmbedding(type='PCA',embeddingType='tSNE',
                colors=r$counts[,gene],shuffle.colors=F,
              alpha=0.1,main=gene)

可以看到HBB这个基因在前面的6和12的表达量比较高，如下所示：

因为这些是pbmc3k的示范数据，所以每个单细胞其实亚群生物学名字是已知的，很容易看。而且我们前面介绍的每个步骤的每个函数都有很大参数可以选择，大家可以自行阅读文档去提高自己：

• https://mran.microsoft.com/snapshot/2021-08-08/web/packages/pagoda2/vignettes/pagoda2.walkthrough.html
• https://cran.r-project.org/web/packages/pagoda2/pagoda2.pdf
• https://github.com/kharchenkolab/pagoda2/blob/main/doc/pagoda2.walkthrough.md

多个单细胞数据集的整合

这个时候需要借助 前面提到的另外一个R包了，名字是con，我们后续再分享哈。

写在文末

我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 jmzeng1314@163.com

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你