NC单细胞文章复现(三)：Clustering

2021-07-02 18:30:31 浏览数 (1)

分享是一种态度

我们上次基于各种marker对1189个细胞进行分类，然而，仅基于marker对细胞进行分类可能是不精确的，特别是考虑到scRNA-seq数据的high dropout rate 。因此，在进行t-SNE降维之前，作者又进一步将细胞进行分类。

首先，因为对不同病人样本的效应进行回归分析，以确保得到的聚类结果不是患者样本之间的异质性导致的。此外，选择了在细胞间分散度高的基因表达子集(平均表达量>0.1)，以增加簇的稳健性。首先计算经验离散值（empirical dispersion value ）去拟合离散度-均值（dispersion parameter ）的关系，为每个基因选择离散度参数（dispersion-mean relationship ）。

代码语言：javascript复制
# 将得到的1189个细胞分类结果更新到sceset_final中
colData(sceset_final)$cell_types_markers <- cell_types_simple
pd_norm <- colData(sceset_final)
#使用monocle包对细胞类型进行无监督聚类方法，monocle_unsup_clust_plots是包装好的一个函数，
HSMM_clustering_ct <- monocle_unsup_clust_plots(sceset_obj = sceset_final, mat_to_cluster = mat_norm,                                                       anno_colors = anno_colors, 
                                                name_in_phenodata = "cluster_allregr_disp", disp_extra = 1,                                                 save_plots = 0,
                                                path_plots = NULL, type_pats = "allpats", regress_pat = 1)


HSMM_clustering_ct$Cluster <- HSMM_clustering_ct$cluster_allregr_disp
table(HSMM_clustering_ct$Cluster)
> table(HSMM_clustering_ct$Cluster)

  1   2   3   4   5   6   7 
 64 141  95  74 202 559  54

我们看到，此时可分为7簇，跟文中的9簇不一样，作者这么解释：由于reduceDimension and clusterCells函数的更新，因此结果会有所不同。

代码语言：javascript复制
# due to changes in Monocle's functions (reduceDimension and clusterCells), the resulting clustering is slightly different from
# the original clustering from the paper. for reproducibility, we read in the original cluster assignment

为了更好地进行下面的学习，我们就拿已经处理好的original clustering跑下去看看。

代码语言：javascript复制
#读取已经处理的的original clustering
original_clustering <- readRDS(file="data/original_clustering.RDS")
HSMM_clustering_ct$Cluster <- original_clustering
table(HSMM_clustering_ct$Cluster)
> table(HSMM_clustering_ct$Cluster)

  1   2   3   4   5   6   7   8   9 
 52  18 156 197 144 119 129 328  46

这时候就有9簇细胞，接下来需要对每1簇细胞再进行细分。

代码语言：javascript复制
# 将得到的9簇细胞分类结果更新到sceset_final中
colData(sceset_final) <- cbind(colData(sceset_final), pData(HSMM_clustering_ct)[,c(96:104)])
colData(sceset_final)$cell_types_cl_all <- colData(sceset_final)$cell_types_markers
pd_norm <- colData(sceset_final)

#首先在每个簇中确定unknown and undecided cell types 
cells_to_assign <- list()
cells_to_reassign <- list()
mean_exprs <- list()
mean_reassign_exprs <- list()

# cluster1: 只有52 epithelial cells
cluster_here <- 1
table(pd_norm$cell_types_markers[which(pd_norm$Cluster == cluster_here)])
mean_exprs[[cluster_here]] <- NULL
mean_reassign_exprs[[cluster_here]] <- NULL
cells_to_assign[[cluster_here]] <- NULL
cells_to_reassign[[cluster_here]] <- NULL

#cluster2：确定clsters2的大致细胞类型：3个epithelial 细胞, 13个macrophage细胞, 1个undecided细胞 and 1个unknown细胞
cluster_here <- 2
table(pd_norm$cell_types_markers[which(pd_norm$Cluster == cluster_here)])

> table(pd_norm$cell_types_markers[which(pd_norm$Cluster == cluster_here)])
   epithelial macrophage  undecided    unknown 
         3         13          1          1 
#计算undecided and unknown cells的免疫标志物的平均表达量
to_assign_cluster <- c("undecided", "unknown")
mean_exprs[[cluster_here]] <- compute_mean_expr_types(types = c("epithelial", "immune", "other"), mat_expr = mat_norm, 
                                                      cells_pos = intersect(which(pd_norm$Cluster == cluster_here), which(pd_norm$cell_types_markers %in% to_assign_cluster)),
                                                      epithelial_markers = epithelial_markers, immune_markers = immune_markers, other_markers = other_markers)
mean_exprs[[cluster_here]]

#对于undecided 和 unknown cell，如果免疫标志物的平均表达量最高，指定为巨噬细胞。
cells_to_assign[[cluster_here]] <- intersect(which(pd_norm$Cluster == cluster_here), which(pd_norm$cell_types_markers %in% to_assign_cluster))[which(apply(mean_exprs[[cluster_here]], 1, function(x){if (x[2] >= max(x)) return(1) else {return(0)}}) == 1)]
pd_norm$cell_types_cl_all[cells_to_assign[[cluster_here]]] <- "macrophage"
table(pd_norm$cell_types_cl_all[which(pd_norm$Cluster == cluster_here)])

#确定clusters2中3个上皮细胞（epithelial cells）免疫标志物的平均表达水平
mean_reassign_exprs[[cluster_here]] <- compute_mean_expr_types(types = c("epithelial", "immune", "other"), mat_expr = mat_norm, 
                                                               cells_pos = intersect(which(pd_norm$Cluster == cluster_here), which(pd_norm$cell_types_markers %in% "epithelial")),
                                                               epithelial_markers = epithelial_markers, immune_markers = immune_markers, other_markers = other_markers)

#类似地，clusters2有3个上皮细胞（epithelial cells）免疫标志物的平均表达最强，也被指定为巨噬细胞。
mean_reassign_exprs[[cluster_here]] <- compute_mean_expr_types(types = c("epithelial", "immune", "other"),                                                                  mat_expr = mat_norm, 
                                                               cells_pos = intersect(which(pd_norm$Cluster                                      == cluster_here), which(pd_norm$cell_types_markers %in% "epithelial")),
                                                               epithelial_markers = epithelial_markers,                                        immune_markers = immune_markers, other_markers = other_markers)
cells_to_reassign[[cluster_here]] <- intersect(which(pd_norm$Cluster == cluster_here), which(pd_norm$cell_types_markers %in% "epithelial"))[which(apply(mean_reassign_exprs[[cluster_here]], 1, function(x){if (x[2] >= max(x)) return(1) else {return(0)}}) == 1)]
pd_norm$cell_types_cl_all[cells_to_reassign[[cluster_here]]] <- "macrophage"
#因此，clusters2由18个巨噬细胞组成。
table(pd_norm$cell_types_cl_all[which(pd_norm$Cluster == cluster_here)])

## cluster3：46 epithelial, 86 stroma, 14 endothelial, 5 undecided and 5 unknown cells.

cluster_here <- 3
table(pd_norm$cell_types_markers[which(pd_norm$Cluster == cluster_here)])
> table(pd_norm$cell_types_markers[which(pd_norm$Cluster == cluster_here)])

endothelial  epithelial      stroma   undecided     unknown 
         14          46          86           5           5 
         
mean_exprs[[cluster_here]] <- NULL
mean_reassign_exprs[[cluster_here]] <- NULL
cells_to_assign[[cluster_here]] <- NULL
cells_to_reassign[[cluster_here]] <- NULL


# cluster 4: 将unknown and undecided cells进一步细分 
cluster_here <- 4
table(pd_norm$cell_types_markers[which(pd_norm$Cluster == cluster_here)])
to_assign_cluster <- c("undecided", "unknown")
mean_exprs[[cluster_here]] <- compute_mean_expr_types(types = c("epithelial", "immune", "other"), mat_expr = mat_norm, 
                                                      cells_pos = intersect(which(pd_norm$Cluster == cluster_here), which(pd_norm$cell_types_markers %in% to_assign_cluster)),
                                                      epithelial_markers = epithelial_markers, immune_markers = immune_markers, other_markers = other_markers)
# clusters3中unknown and undecided cells共有5个细胞的epithelial markers的平均表达最强，被指定为epithelial cells

cells_to_assign[[cluster_here]] <- intersect(which(pd_norm$Cluster == cluster_here), which(pd_norm$cell_types_markers %in% to_assign_cluster))[which(apply(mean_exprs[[cluster_here]], 1, function(x){if (x[1] >= max(x)) return(1) else {return(0)}}) == 1)]
pd_norm$cell_types_cl_all[cells_to_assign[[cluster_here]]] <- "epithelial"
table(pd_norm$cell_types_cl_all[which(pd_norm$Cluster == cluster_here)])

# 确定个stroma cell是否高表达epithelial marker，若epithelial marker表达较高，则被指定为 epithelial cells
mean_reassign_exprs[[cluster_here]] <- compute_mean_expr_types(types = c("epithelial", "immune", "other"), mat_expr = mat_norm, 
                                                               cells_pos = intersect(which(pd_norm$Cluster == cluster_here), which(pd_norm$cell_types_markers %in% "stroma")),
                                                               epithelial_markers = epithelial_markers, immune_markers = immune_markers, other_markers = other_markers)
# only re-assign the stroma cells if their epithelial expression is higher
cells_to_reassign[[cluster_here]] <- intersect(which(pd_norm$Cluster == cluster_here), which(pd_norm$cell_types_markers %in% "stroma"))[which(apply(mean_reassign_exprs[[cluster_here]], 1, function(x){if (x[1] >= max(x)) return(1) else {return(0)}}) == 1)]
pd_norm$cell_types_cl_all[cells_to_reassign[[cluster_here]]] <- "epithelial"
table(pd_norm$cell_types_cl_all[which(pd_norm$Cluster == cluster_here)])

cluster 5:46 epithelial, 4 stroma, 52 T, 3 B, 11 undecided and 28 unknown cells
cluster_here <- 5
table(pd_norm$cell_types_markers[which(pd_norm$Cluster == cluster_here)])
mean_exprs[[cluster_here]] <- NULL
mean_reassign_exprs[[cluster_here]] <- NULL
cells_to_assign[[cluster_here]] <- NULL
cells_to_reassign[[cluster_here]] <- NULL
........

接下来的第6、7、8和9簇细胞，按照同样的办法细分，代码太长了，我就不贴上来了。在这里先小结一下，我们开始看到的是1189个细胞，然后将这些细胞用无聚类监督分析先分为9个簇，因为相似的细胞更容易成为一个簇，但是簇里面的细胞类型我们已经基于文献的markers分类好了，若簇里面有多种细胞类型的细胞存在，那么我们需要进一步鉴定他们跟簇里面那群最大比例的细胞是不是就是同一类细胞（不然为何他们会聚集到了同一簇），还是说有自己的独特身份？但是我们怎么定义那个最大比例的细胞，作者将这个比例定为80%，举个例子，在clusters4总共有197个细胞（185 epithelial，3 stroma, 4 undecided and 5 unknown cells），其中 epithelial细胞占了94%，超过了80%的比例，那么其他12个细胞有可能“伪装”了自己，其真实身份有可能是epithelial细胞，不然怎么解释他们跟epithelial细胞聚集到了一起呢？这时候就需要进一步确定其他12个细胞epithelial markers 平均表达量是否高于其他的markers,是的话他们的“庐山真面目”就是epithelial细胞，若没有，那么就是冤枉到它了，让它保持原来的身份。而对于clusters3（ 46 epithelial, 86 stroma, 14 endothelial, 5 undecided and 5 unknown cells）来说，没有1种细胞的比例超过80%，因此他们都可以保持自己的身份。我感觉这个过程还是蛮有趣的，在对细胞核实身份的同时，其实就是不断接近客观世界的过程。今天这些内容已经蛮多了，下一节我们再继续学习。

往期回顾

OSCA单细胞数据分析笔记9—Clustering

scPhere——用地球仪来展示降维结果

2021第一期生信入门微信群答疑精选200题

这50个ggplot2现成图表你居然没有从头到尾自己画一遍

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

  • 生信爆款入门-2021第4期
  • 数据挖掘（GEO,TCGA,单细胞）2021第4期
  • 明码标价之共享96线程384G内存服务器
geo marker na

0 人点赞