降维后进行细胞聚类分群

 今天是生信星球陪你的第461天

   大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~

   就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~

   这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！

豆豆写于19.9.28
目标：利用去除技术噪音的PCA结果进行细胞聚类，推断细胞亚群。

来自Bioconductor放出的单细胞数据分析教程
https://bioconductor.org/packages/release/workflows/vignettes/simpleSingleCell/inst/doc/reads.html

前面去除技术噪音中使用了denoisePCA函数，最后保留了24个主成分，存储在reducedDim(sce, "PCA")中，其中行是细胞，列是主成分（从PC1到PC24）。

思路就是：先提取PCA结果中细胞与主成分PC1、PC2的坐标，然后根据细胞在空间上的欧式距离进行层次聚类（对行进行操作，也就是这里的细胞），并且使用Ward标准去降低每个cluster内的总体方差，保证了每个cluster中的细胞都含有相似的基因表达模式

# 三步走：取PCA结果=》计算细胞欧式距离=》降低cluster内差异
pcs <- reducedDim(sce, "PCA")
my.dist <- dist(pcs)
my.tree <- hclust(my.dist, method="ward.D2")

接着利用dynamic tree cut的方法得到cluster的数量，在复杂的聚类数中相比cutree更有优势 (Langfelder, Zhang, and Horvath 2008) 。但是函数计算的毕竟是统计角度的cluster，我们还需要人工干预帮忙校正一下。比如设置：cutHeight、minClusterSize 。

library(dynamicTreeCut)
my.clusters <- unname(cutreeDynamic(my.tree, distM=as.matrix(my.dist), 
    minClusterSize=10, verbose=0))
# minClusterSize默认值是20，这里设置小一点是避免一些偏远的小群体错误地凑在一起

然后，利用得到的细胞分群结果与已知的分类因子进行结合，首先看得cluster与批次之间的关系：

table(my.clusters, sce$Plate)
##            
## my.clusters 20160113 20160325
##           1       41       39
##           2       19       20
##           3       16       11
##           4       10       14
##           5        5        8

发现其中每个cluster基本都包含两个批次，并且数量都差不多，表明这里的聚类结果不受批次效应影响。

再看cluster与生物差异之间的关系：

table(my.clusters, sce$Oncogene)
##            
## my.clusters induced control
##           1      80       0
##           2       0      39
##           3       0      27
##           4       0      24
##           5      13       0

发现其中的cluster在处理和对照之间存在很大差异，更加验证了是由于生物因素导致的分群

得到cluster后，用tsne结果进行映射

sce$cluster <- factor(my.clusters)
plotTSNE(sce, colour_by="cluster") + fontsize

图27

检查分群结果

一般来说，分群结果分的越开，说明效果越好。如果感觉用肉眼观察tsne图不容易判断的话，还有一种方法可以帮助我们判断，学名叫”silhouette width“（轮廓图）。先上代码和图：

library(cluster)
clust.col <- scater:::.get_palette("tableau10medium") # hidden scater colours
sil <- silhouette(my.clusters, dist = my.dist)
sil.cols <- clust.col[ifelse(sil[,3] > 0, sil[,1], sil[,2])]
sil.cols <- sil.cols[order(-sil[,1], sil[,3])]
plot(sil, main = paste(length(unique(my.clusters)), "clusters"), 
    border=sil.cols, col=sil.cols, do.col.sort=FALSE) 

图28

• 这个图是对上面得到的各个cluster中的细胞做的barplot，每个cluster是一种颜色。

• 如果一个细胞的横坐标width值为正并且它越大，表示相比于其他clusters的细胞，这个细胞和它所在cluster中的其他细胞更接近；如果width为负，就表示这个cluster的这个细胞和其他cluster的细胞更接近，也就是说分群效果不好

• 利用这个图，我们就能调整之前的参数，来调整分群效果，不过也不需要太过纠结完美的分群结果

• 从图中可以看到，大多数细胞的width数值比较小，表明cluster之间的分离效果一般。很有可能是”过度分群“，比如原来都属于treat组的cluster又继续分成了小cluster，那么它们之间分离效果肯定要差。

作者建议

• 聚类也可以用quickCluster函数，关于它的利弊：more robust to noise and normalization errors, but is also less sensitive to subtle changes in the expression profiles

---

点击底部的“阅读原文”，获得更好的阅读体验哦😻

初学生信，很荣幸带你迈出第一步。

我们是生信星球，一个不拽术语、通俗易懂的生信知识平台。由于是2018年新号，竟然没有留言功能。需要帮助或提出意见请后台留言、联系微信或发送邮件到[e-mail]，每一条都会看到的哦~