NC单细胞文章复现(三)：复杂热图

生信技能树jimmy

发布于 2021-07-02 18:20:19

2.8K2

发布于 2021-07-02 18:20:19

分享是一种态度

我们这次直接拿GSE118390上已经normalized 的数据进行下游分析。首先我们先看看文献的这张复杂热图，哈哈，这张热图画得真是好看。左边是不同的markers基因对应的细胞类型，上边是6个TNBC病人总共1189个细胞，cycling指的是处于细胞分裂期的样本，下边还有1189个细胞对应的是CD45阳性还是CD阴性的。代码作者都已经有提供，但是如果按照作者的代码跑，仍然会出现error，我把代码一些参数修改了一下并进行翻译，再结合文献的描述，跟大家一起学习。我们看看怎么画呗。

一.复杂热图

作者首先将已建立的用于定义免疫细胞、内皮细胞和基质细胞以及关键乳腺上皮亚群（basal cells, luminal progenitors, andmature luminal cells）的特定基因集去注释1189个细胞的类型，然后使用聚类将1189个细胞中的1112个分为非上皮（non-epithelial：n=244）和上皮（epithelial：n=868）两种类型。已知TNBC的一个亚类会表现出大量的免疫细胞浸润，CD45阴性细胞（CD45-unselected）中含有大量的免疫细胞亚群，其中大部分是T淋巴细胞，其次是巨噬细胞。间质细胞(stromal cells)成分在TNBC病例之间也存在差异，这种细胞在每个肿瘤中占总细胞的比例<15%。内皮细胞（endothelial cells）是少数群体，最多占肿瘤细胞的4%

#加载包
library(here)
source(here::here("code", "1.libraries.R"))
#funcs.R和funcs_markers是已经包装好的函数，我们接下来会有稍微讲解一下。
source(here::here("code", "funcs.R"))
source(here::here("code", "funcs_markers.R"))
#读取数据和表型信息
mat_norm <- read.table("./data/norm_data.txt", sep = "\t", header = TRUE)#1189个细胞，13280个基因
pd_norm <- read.table("./data/pd_norm.txt", sep = "\t", header = TRUE, stringsAsFactors = FALSE)
fd_norm <- read.table("./data/fd_norm.txt", sep = "\t", header = TRUE, stringsAsFactors = FALSE)
#创建SingleCellExperiment对象
sceset_final <- SingleCellExperiment(assays = list(exprs = as.matrix(mat_norm)),
                                     colData = pd_norm, 
                                     rowData = fd_norm)

#对各个指标创建创建图像矢量
epithelial_col <- brocolors("crayons")["Maroon"]
basal_epithelial_col <- brocolors("crayons")["Red"]
luminal_epithelial_col <- brocolors("crayons")["Sunset Orange"]
luminal_progenitor_col <- brocolors("crayons")["Salmon"]

stroma_col <- brocolors("crayons")["Aquamarine"]
endothelial_col <- brocolors("crayons")["Wisteria"]

PTPRC_col <- brocolors("crayons")["Inchworm"]
t_cell_col <- brocolors("crayons")["Screamin' Green"]
b_cell_col <- brocolors("crayons")["Fern"]
macrophage_col <- brocolors("crayons")["Tropical Rain Forest"]

marker_cols <- c("epithelial" = unname(epithelial_col), "basal epithelial" = unname(basal_epithelial_col), 
                 "luminal epithelial" = unname(luminal_epithelial_col), "luminal progenitor" = unname(luminal_progenitor_col),
                 "stroma" = unname(stroma_col), "endothelial" = unname(endothelial_col), 
                 "immune" = unname(PTPRC_col), "T cell" = unname(t_cell_col), "B cell" = unname(b_cell_col), 
                 "macrophage" = unname(macrophage_col))
cycling_mel_cols <- c("non-cycling" = "gainsboro",
                      "cycling" = unname(brocolors("crayons")["Mulberry"]))
depletion_cols <- c("depleted" = unname(brocolors("crayons")["White"]), "not depleted" = unname(brocolors("crayons")["Red"]))
pats_cols <- c("PT039" = unname(brocolors("crayons")["Orange Red"]), "PT058" = unname(brocolors("crayons")["Orange"]), 
               "PT081" = unname(brocolors("crayons")["Pink Flamingo"]), "PT084" = unname(brocolors("crayons")["Fern"]), 
               "PT089" = unname(brocolors("crayons")["Blue Violet"]), "PT126" = unname(brocolors("crayons")["Sky Blue"]))
tsne_cols <- c("epithelial" = unname(basal_epithelial_col), "stroma" = unname(stroma_col), "endothelial" = unname(endothelial_col),
               "Tcell" = unname(t_cell_col), "Bcell" = unname(b_cell_col), "macrophage" = unname(macrophage_col))
#tsne_cols_unknown <- c(tsne_cols, "unknown" = "black", "undecided" = "black")
anno_colors <- list("marker" = marker_cols, "cycling" = cycling_mel_cols, "immune depletion" = depletion_cols, 
                    "patient" = pats_cols, "tsne" = tsne_cols)

接下来注释细胞周期基因，采用了黑色素瘤的细胞周期markers基因（Tirosch et al 2016）

#读取黑色素瘤细胞周期基因列表
melanoma_cellcycle <- read.table("./data/melanoma_cellcycle.txt", sep = "\t", header = TRUE, stringsAsFactors = FALSE)
melanoma_g1s <- melanoma_cellcycle$G1S #G1S期markers有43个

#接下来这个match_clean_vector_genes函数是包装在funs.R中，
#具体来说是先删除#melanoma_g1s的重复基因，删除空白字符，转换成data.frame，匹配#mat_norm和 melanoma_g1s的marker基因，再判断一下是否有缺失值。
melanoma_g1s <- match_clean_vector_genes(mat_norm, melanoma_g1s)
#avg_expr_genes是用apply包装的一个函数，即算出mat_norm中1189个细胞的G1S期评分（即G1S期marker表达量占全部细胞的比例）
scores_g1s <- avg_expr_genes(mat_norm, melanoma_g1s$index)
#同样地，算出G2M期的评分
melanoma_g2m <- melanoma_cellcycle$G2M
melanoma_g2m <- match_clean_vector_genes(mat_norm, melanoma_g2m)
scores_g2m <- avg_expr_genes(mat_norm, melanoma_g2m$index)

我们看看文献是怎么定义cycling和non-cycling

cycling cells： 1、细胞的G1S期评分≥其中位数+2MAD和细胞的G2M期评分≤其中位数+2 MAD。 2、细胞的G1S期评分≤其中位数2 MADs 和细胞的G2M期评分≥2 MADs。 3、细胞的G1S期评分≥其中位数+2倍的绝对中位差和细胞的G2M期评分≥其中位数+2倍的绝对中位差。
non-cycling cells：细胞的G1S期评分<其中位数+2MAD和细胞的G2M期评分<其中位数+2 MAD

#这个for循环值得好好学习，以后要是我们要是有自己的数据拿来注释细胞周期，也许可以用到。
cycling_mel <- rep(NA, length(scores_g1s))
for (i in 1:length(cycling_mel)) {
  if (scores_g1s[i] >= (median(scores_g1s) + 2 * mad(scores_g1s)) && scores_g2m[i] < (median(scores_g2m) + 2 * mad(scores_g2m)))
    cycling_mel[i] <- "cycling"
  if (scores_g1s[i] < (median(scores_g1s) + 2 * mad(scores_g1s)) && scores_g2m[i] >= (median(scores_g2m) + 2 * mad(scores_g2m)))
    cycling_mel[i] <- "cycling"
  if (scores_g1s[i] < (median(scores_g1s) + 2 * mad(scores_g1s)) && scores_g2m[i] < (median(scores_g2m) + 2 * mad(scores_g2m)))
    cycling_mel[i] <- "non-cycling"
  if (scores_g1s[i] >= (median(scores_g1s) + 2 * mad(scores_g1s)) && scores_g2m[i] >= (median(scores_g2m) + 2 * mad(scores_g2m)))
    cycling_mel[i] <- "cycling"
}
table(cycling_mel)
> table(cycling_mel)
cycling_mel
    cycling non-cycling  
        216         973

此时我们看到，cycling cells有216个，non-clcling cells有973个，咦？怎么文献中说cycling cells有214个，non-clcling cells有898个，难道中间哪一步有瑕疵？别急，我们慢慢往下看。

#把cycling和non-cycling的信息添加到sceset_final中,方便后续处理
colData(sceset_final)$mel_scores_g1s <- scores_g1s
colData(sceset_final)$mel_scores_g2m <- scores_g2m
colData(sceset_final)$cycling_mel <- cycling_mel
pd_norm <- colData(sceset_final)
pd_norm <- as.data.frame(pd_norm)
#我们可以顺便看看1189个细胞的cycling和non-cycling的信息
b=colData(sceset_final)
cell_cycling_non_cycling=cbind(rownames(b),b$cycling_mel)
cell_cycling_non_cycling=as.data.frame(cell_cycling_non_cycling)
cell_cycling_non_cycling
#                  V1          V2
#1        PT089_P1_A01 non-cycling
#2        PT089_P1_A02 non-cycling
#3        PT089_P1_A03 non-cycling
#4        PT089_P1_A04 non-cycling
#         ........................
#初始化部分矩阵
patients_now <- c()
mats_now <- list()
pds_now <- list()

for (i in 1:length(unique(pd_norm$patient))) {  #
  patients_now[i] <- sort(unique(pd_norm$patient))[i] #病人ID按照字母和数值从小到大排序
  mats_now[[i]] <- mat_norm[, pd_norm$patient == patients_now[i]]
  pds_now[[i]] <- pd_norm[pd_norm$patient == patients_now[i],]
}
names(mats_now) <- patients_now #包含了6个样本各自的13280个基因表达量
names(pds_now) <- patients_now #包含了6个样本各自的表型信息

二.identifying cell types

接着这篇文章作者开始进行细胞注释，分为两个步骤。

使用于文献的特定表达标记的基因列表来定义细胞类型。
clustering

对于第1步骤，使用的基因列表是来自(Tirosh et al., 2016）整理好的用来注释四种细胞类型（epithelial，stroma，endothelial，immune）的53个marker基因。

#我们读取已经整理好的markers，开始是有53个markers
all_markers <- read.table("data/markers_clean.txt", sep = "\t", header = TRUE, stringsAsFactors = FALSE)
length(all_markers$gene)
> length(all_markers$gene) #53个markers
[1] 53
#j接着看看53个marker有没有包含在1189个细胞中的13280个基因里面
#若没有的话，将其删除。经过这一步过滤，将原来53个marker，过滤成49个marker。
markers <- unique(all_markers[all_markers$gene %in% rowData(sceset_final)$hgnc_symbol, ])
length(markers$gene)
> length(markers$ge

接下是画热图，这张图的信息量很大，我们慢慢来品。

#首先是先将热图的注释信息弄好，这一步骤是将49个marker基因在热图中所对应的细胞类型表弄好。
markers$type_heatmap <- markers$type
markers$type_heatmap[which(markers$type_heatmap == "luminalprogenitor")] <- "luminal progenitor"
markers$type_heatmap[which(markers$type_heatmap == "luminalepithelial")] <- "luminal epithelial"
markers$type_heatmap[which(markers$type_heatmap == "basalepithelial")] <- "basal epithelial"
markers$type_heatmap[which(markers$type_heatmap %in% c("EPCAM", "EGFR", "CDH1"))] <- "epithelial"
markers$type_heatmap[which(markers$type_heatmap == "Bcell")] <- "B cell"
markers$type_heatmap[which(markers$type_heatmap == "Tcell")] <- "T cell"
markers$type_long_heatmap <- markers$type_long
markers$type_long_heatmap[which(markers$type == "stroma")] <- "stroma"
markers$type_long_heatmap[which(markers$type == "endothelial")] <- "endothelial"
#接下来画 fig 1b，第一个for循环是用一个replace函数，将49个markers的细胞类型附上对应颜色
#颜色的注释信息在我们之前的创建图像矢量的时候就已经建立好
colors_markers_ch <- markers$type_heatmap
for (i in c(1:length(names(anno_colors$marker)))) {
  colors_markers_ch <- replace(colors_markers_ch, colors_markers_ch == names(anno_colors$marker)[i], anno_colors$marker[i])
}
##调整将热图左边的细胞类型注释条顺序，从上到下依次为epithelial，stroma，endothelial，immune
splits_ch <- as.factor(markers$type_long_heatmap)
splits_ch <- factor(splits_ch, levels(splits_ch)[c(2,4,1,3)])
##再调整将热图下边的细胞类型注释条顺序
#顺序依次为Epithelial、Basal epithelial、Luminal epithelial、Luminal progenitor
#StromaEndothelial、Immune、T cell、B cell 和 Macrophage
colors_anno_markers_ch <- as.factor(markers$type_heatmap)
colors_anno_markers_ch <- factor(colors_anno_markers_ch, levels(colors_anno_markers_ch)[c(4,2,5,6,9,3,8,10,1,7)])
##注释条的颜色，大小，标题等。
b=data.frame (colors_anno_markers_ch)
ha_rows <- HeatmapAnnotation(df = data.frame(type=colors_anno_markers_ch),
                             annotation_legend_param = list(type = list(ncol = 2, title = "cell type", title_position = "topcenter")),
                             which = "row", col = list("type"=anno_colors$marker),annotation_width = unit(3, "mm"))

接下来这一步是将复杂热图各自注释信息例如聚类方向、字体位置、题目、边框等信息整理好，但是如果按照作者的代码跑，会出现莫名其妙的error.

cycling_now <- list()
depletion_now <- list()
ha_cols_up <- list()
ha_cols_bottom <- list()
for (i in 1:length(patients_now)) {
  cycling_now[[i]] <- pds_now[[i]][,"cycling_mel"]
  
  if (i == 1)
    ha_cols_up[[i]] <- HeatmapAnnotation(data.frame(cycling = cycling_now[[i]]), 
                                         col = list(cycling = anno_colors$cycling), 
                                         show_annotation_name = TRUE, 
                                         annotation_name_side = "left", annotation_name_gp = gpar(fontsize = 11),
                                         annotation_legend_param = list(list(title_position = "topcenter",
                                                                             title = c("cycling status"))),
                                         gap = unit(c(1, 1), "mm"))
  if (i > 1)
    ha_cols_up[[i]] <- HeatmapAnnotation(data.frame(cycling = cycling_now[[i]]), 
                                         col = list(cycling = anno_colors$cycling), 
                                         show_legend = FALSE,
                                         gap = unit(c(1, 1), "mm"))
}
> 错误: annotations should have names.

我不断去尝试各种方法。好在最后顺利解决。

首先我们安装最新版本的ComplexHeatmap package，如果你只是用bioconductor的官方代码去安装，没法安装到最新版本的ComplexHeatmap package，解决方法可以是利用devtools package去github上安装。
HeatmapAnnotation函数有些参数需要修改，我把修改好的贴上来了。至于具体的画图参数，我就不细讲了，有兴趣的小伙伴可以去看一看ComplexHeatmap package的使用说明书。

library(devtools)
install_github("jokergoo/ComplexHeatmap")
cycling_now <- list()
depletion_now <- list()
ha_cols_up <- list()
ha_cols_bottom <- list()
#先注释cycing和non-cycling的信息
for (i in 1:length(patients_now)) {
  cycling_now[[i]] <- pds_now[[i]][,"cycling_mel"] 
  
  if (i == 1)
    ha_cols_up[[i]] <- HeatmapAnnotation(df=data.frame(cycling = cycling_now[[i]]), 
                                         col = list(cycling = anno_colors$cycling), 
                                         which = "column", 
                                         show_annotation_name = TRUE,  #第一个病人ID的cycing字体展示出来
                                         annotation_name_side = "right", #将cycling这个单词放在最右边
                                         annotation_name_gp = gpar(fontsize = 11),#字体大小
                                         annotation_legend_param = list(list(title_position = "topcenter",
                                                                             title = c("cycling status"))), #热图下边关于cycling的注释条信息
                                         annotation_name_align = F,
                                         gap = unit(c(1, 1), "mm")) #注释之间的距离
  if (i > 1)
    ha_cols_up[[i]] <- HeatmapAnnotation(df=data.frame(cycling = cycling_now[[i]]), 
                                         col = list(cycling = anno_colors$cycling), 
                                         show_legend = F,#从第2个病人样本ID开始，不要展示cycling字体
                                         
                                         gap = unit(c(1, 1), "mm"))
  
  }

#注释CD45 status的信息
for (i in 1:length(patients_now)) {
  depletion_now[[i]] <- pds_now[[i]][,"depletion_batch"]
  
  if (i == 1)
    ha_cols_bottom[[i]] <- HeatmapAnnotation(df=data.frame('CD45' = depletion_now[[i]]), 
                                             col = list('CD45' = c("depleted_yes" = "gainsboro", "depleted_no" = "gray54")),
                                             show_annotation_name = TRUE,
                                             which = "column",
                                             annotation_name_side = "right", annotation_name_gp = gpar(fontsize = 11),
                                             annotation_legend_param = list(title = "CD45 status", 
                                                                            title_position = "topcenter", 
                                                                            at = c("depleted_yes", "depleted_no"),
                                                                            labels = c("CD45 depleted","CD45 unselected")),
                                             annotation_name_align = T,
                                             show_legend = TRUE,
                                             gap = unit(c(1), "mm"))
  
  
  if (i == 2)
    ha_cols_bottom[[i]] <- HeatmapAnnotation(df=data.frame('CD45' = depletion_now[[i]]), 
                                             col = list('CD45' = c("depleted_yes" = "gainsboro", "depleted_no" = "gray54")),
                                             show_annotation_name = FALSE, 
                                             show_legend = FALSE,
                                             which = "column",
                                             annotation_name_side = "left", annotation_name_gp = gpar(fontsize = 11),
                                             annotation_legend_param = list(title = "CD45 status", 
                                                                            title_position = "topcenter", 
                                                                            at = c("depleted_yes", "depleted_no"),
                                                                            labels = c("CD45 depleted","CD45 unselected")),
                                             gap = unit(c(1), "mm"))
  if (i > 2)
    ha_cols_bottom[[i]] <- HeatmapAnnotation(df=data.frame('CD45' = depletion_now[[i]]), 
                                             col = list('CD45' = c("depleted_yes" = "gainsboro", "depleted_no" = "gray54")),
                                             show_legend = FALSE,
                                             which = "column",
                                             gap = unit(c(1), "mm")) 
  
}
#画图，具体参数大家自己有兴趣可以去查看Heatmap函数
names(cycling_now) <- patients_now
names(depletion_now) <- patients_now
names(ha_cols_up) <- patients_now
names(ha_cols_bottom) <- patients_now



a=ha_rows + 
  Heatmap(mats_now[[1]][match(markers$gene,rownames(mats_now[[1]])),], 
          cluster_rows = FALSE, show_column_dend = TRUE, show_column_names = FALSE,
          name = patients_now[1], clustering_distance_columns = "euclidean", row_names_side = "right", 
          row_names_gp = gpar(fontsize = 10, col = colors_markers_ch),
          split = splits_ch, gap = unit(1, "mm"), column_title = patients_now[1], 
          column_title_gp = gpar(fontsize = 11),
          row_title_gp = gpar(font = 11), top_annotation = ha_cols_up[[1]],
          heatmap_legend_param = list(title = "expression", title_position = "topcenter", color_bar = "continuous"),
          bottom_annotation = ha_cols_bottom[[1]])


ht_list <- ha_rows + 
  Heatmap(mats_now[[1]][match(markers$gene,rownames(mats_now[[1]])),], 
          cluster_rows = FALSE, show_column_dend = TRUE, show_column_names = FALSE,
          name = patients_now[1], clustering_distance_columns = "euclidean", row_names_side = "right", 
          row_names_gp = gpar(fontsize = 10, col = colors_markers_ch),
          split = splits_ch, gap = unit(1, "mm"), column_title = patients_now[1], 
          column_title_gp = gpar(fontsize = 11),
          row_title_gp = gpar(font = 11), top_annotation = ha_cols_up[[1]],
          heatmap_legend_param = list(title = "expression", title_position = "topcenter", color_bar = "continuous"),
          bottom_annotation = ha_cols_bottom[[1]]) + 
  Heatmap(mats_now[[2]][match(markers$gene,rownames(mats_now[[1]])),],
          cluster_rows = FALSE, show_column_dend = TRUE, show_column_names = FALSE,
          name = patients_now[2], clustering_distance_columns = "euclidean",
          show_row_names = FALSE, column_title = patients_now[2],
          column_title_gp = gpar(fontsize = 11),
          show_heatmap_legend = FALSE,
          row_title_gp = gpar(font = 11), top_annotation = ha_cols_up[[2]], bottom_annotation = ha_cols_bottom[[2]],
          gap = unit(1, "mm")) +
  Heatmap(mats_now[[3]][match(markers$gene,rownames(mats_now[[3]])),],
          cluster_rows = FALSE, show_column_dend = TRUE, show_column_names = FALSE,
          name = patients_now[3], clustering_distance_columns = "euclidean",
          show_row_names = FALSE, column_title = patients_now[3],
          column_title_gp = gpar(fontsize = 11),
          show_heatmap_legend = FALSE,
          row_title_gp = gpar(font = 11), top_annotation = ha_cols_up[[3]], bottom_annotation = ha_cols_bottom[[3]],
          gap = unit(1, "mm")) +
  Heatmap(mats_now[[4]][match(markers$gene,rownames(mats_now[[4]])),],
          cluster_rows = FALSE, show_column_dend = TRUE, show_column_names = FALSE,
          name = patients_now[4], clustering_distance_columns = "euclidean",
          show_row_names = FALSE, column_title = patients_now[4], 
          column_title_gp = gpar(fontsize = 11),
          show_heatmap_legend = FALSE,
          row_title_gp = gpar(font = 11), top_annotation = ha_cols_up[[4]], bottom_annotation = ha_cols_bottom[[4]],
          gap = unit(1, "mm")) +
  Heatmap(mats_now[[5]][match(markers$gene,rownames(mats_now[[5]])),], 
          cluster_rows = FALSE, show_column_dend = TRUE, show_column_names = FALSE,
          name = patients_now[5], clustering_distance_columns = "euclidean",
          show_row_names = FALSE, column_title = patients_now[5], 
          column_title_gp = gpar(fontsize = 11), 
          show_heatmap_legend = FALSE,
          row_title_gp = gpar(font = 11), top_annotation = ha_cols_up[[5]], bottom_annotation = ha_cols_bottom[[5]],
          gap = unit(1, "mm")) + 
  Heatmap(mats_now[[6]][match(markers$gene,rownames(mats_now[[6]])),], 
          cluster_rows = FALSE, show_column_dend = TRUE, show_column_names = FALSE,
          name = patients_now[6], clustering_distance_columns = "euclidean",
          row_names_gp = gpar(fontsize = 10, col = colors_markers_ch), 
          split = splits_ch, column_title = patients_now[6], 
          column_title_gp = gpar(fontsize = 11),
          row_title_gp = gpar(font = 11), top_annotation = ha_cols_up[[6]], bottom_annotation = ha_cols_bottom[[6]],
          show_heatmap_legend = FALSE,
          gap = unit(1, "mm"))

draw(ht_list, gap = unit(0.1, "cm"), heatmap_legend_side = "bottom", annotation_legend_side = "bottom")

这样子一张复杂热图总算完成了，基本跟文献是一样的。但是我们的学习还没有结束。接着需要定义markers的表达量大于多少才算是有表达，设置的阈值是1.

接着需要确定这些细胞是否是：

上皮性来源细胞，需满足一下2个条件之一：
1. 至少2个上皮性markers,
2. 若该细胞只有一个乳腺上皮标志物：EpCAM、KRT8、KRT18、KRT19，对于该标志物，该标志物在该患者中须有超过50%的细胞有表达。
免疫成分细胞:如果一个细胞表现出特殊的免疫类型（T细胞，B细胞，巨噬细胞）：
1. 只有一种类型的特异性免疫标记物（至少2个）；
2. 有PTPRC markers和该类型的唯一特异性免疫标记物（至少1个）；
3. 至少3个这种类型的免疫标记，最多1个其他类型的免疫标记。比较特殊的是若该细胞有PTRC，则再看有没有别的marker，若还有1个以上的免疫细胞，那么需要继续细分免疫细胞亚群（macrophage，T cell，B cell），若没有其他marker，则笼统地定义为免疫细胞。
基质成分细胞：该细胞仅有基质标志物；该细胞有最多3个基质标志物，且至少有1个内皮标志物。
内皮成分细胞：该细胞仅有内皮标志物；或者该细胞有至少3个内皮标志物和最多1个基质标志物。

具体的话可以参考

接下来的lists_markers、decide_is_epithelial 和decide_is_immune这3个函数是包装在了funcs_markers.R里面，都是用for循环实现的。简单来说：lists_markers函数先把细胞分为epithelial_cells、immune_cells和other_cells，然后再进行细分。decide_is_epithelial函数进一步确定epithelial_cells是否有2个上皮性markers，如果有，将该细胞标记为1，如果没有，将该细胞标记为0。decide_is_immune函数比较复杂，如果该细胞没有immune markers或者仅仅只有一个markers，标记为0；如果该细胞至少有2个immune markers，并且没有PTPRC markers，那么返回该markers对应的细胞类型；如果该细胞至少有1种类型的2个immune markers,并且有PTPRC markers，那么返回该markers对应的细胞类型

如果该细胞有1种细胞类型的3个immune markers，而且还含有另外1种类型的1个markers，则返回有3个immune markers的类型。如果大家看不懂，可结合上面那张图和代码去看。

thresh <- 1 #表达阈值设置为1
cells_markers <- lists_markers(mat_norm, thresh, markers)
ncol(mat_norm)
#定义上皮细胞
epithelial_markers <- cells_markers$epithelial_cells
is_epithelial <- decide_is_epithelial(epithelial_markers)
#定义免疫细胞
immune_markers <- cells_markers$immune_cells
is_immune <- decide_is_immune(immune_markers)
 #定义其他类型的细胞类型：stroma, endothelial or adipocytes cells
other_markers <- cells_markers$other_cells
is_other <- decide_is_other(other_markers)
#作者在这一步里继续细化注释上皮细胞的流程，如果这个细胞只表达1个epithelial marker，
#先查看这个markers在样本中的表达分布,如果分布在超过50%的细胞，则标记为1，若没有标记为0
one_epithelial_marker <- expression_one_epithelial_marker(mat_norm, pd_norm, is_epithelial, epithelial_markers, "pats", 0.5)
is_epithelial[which(one_epithelial_marker$is_epithelial_extra == 1)] <- 1
one_epithelial_marker <- expression_one_epithelial_marker(mat_norm, pd_norm, is_epithelial, epithelial_markers, "pats", 0.5)
is_epithelial[which(one_epithelial_marker$is_epithelial_extra == 1)] <- 1

is_epithelial_simple <- is_epithelial
is_epithelial_simple[which(is_epithelial == 1)] <- "epithelial"
#将我们之前分好的有多种细胞类型的immune_mix标记为0
is_immune_simple <- is_immune
is_immune_simple[which(is_immune == "immune_mix")] <- 0
#将我们之前分好的有多种细胞类型的other_cell标记为0
is_other_simple <- is_other
is_other_simple[which(is_other == "other_mix")] <- 0

cells_types <- paste(is_epithelial_simple, is_immune_simple, is_other_simple, sep = "_")
names(cells_types) <- names(is_epithelial)


#若出现0-0-0，则该细胞是unknown，若出现epithelial-stroma-0,该细胞定义为epithelial cells，理由是：
#同时具备epithelial和stroma markers的细胞一般跟EMT有关，可归入epithelial一类。
cell_types <- sapply(strsplit(cells_types, "_"), function(x){
  if (sum(x == 0) == 3) return("unknown") else 
    if (sum(x == 0) == 2) return(setdiff(x, "0")) else
      if (sum(c("epithelial", "stroma", "0") %in% x) == 3) return("epithelial") else
        return(paste(setdiff(x, "0"),collapse = "_"))})
cell_types_simple <- cell_types
#若某个细胞的被注释到的细胞类型超过1，则该细胞定义为undecided.
cell_types_simple[which(sapply(strsplit(cell_types, "_"), length) > 1)] <- "undecided"
#此时我们再看，跟文献注释到的细胞数目一模一样。
> table(cell_types_simple)
cell_types_simple
      Bcell endothelial  epithelial  macrophage      stroma       Tcell   undecided 
         19          14         831          51          95          53          31 
    unknown 
         95

这里面for循环很多，我也只是看懂个大概意思，可能写得也不怎么清楚，具体里面的内在逻辑还需要结合文章的方法去不断理解消化。不过我们也是学习到了一些细胞注释的知识，例如在注释上皮细胞类型的时候，需要注意EPCAM, KRT8, KRT18,和KRT19这些基因，这四个基因在上皮细胞高度表达，可单独拿出来作为上皮细胞标志。PTPRC也称CD45，最初被称为白细胞共同抗原，可作为免疫细胞亚群的第一次分群的标志。若细胞同时具备epithelial和stroma两种markers，则应该被归为上皮细胞亚群。

我们下一篇会开始学习作者怎么又把1189个细胞减到1112个有“名字”的细胞，然后开始t-SNE降维聚类等等。哈哈哈，学无止境，我们继续结合作者的代码和文献内容，慢慢学吧.....