NC单细胞文章复现(三)：复杂热图 / 开普饭

我们这次直接拿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个条件之一：

至少2个上皮性markers,
若该细胞只有一个乳腺上皮标志物：EpCAM、KRT8、KRT18、KRT19，对于该标志物，该标志物在该患者中须有超过50%的细胞有表达。

免疫成分细胞:如果一个细胞表现出特殊的免疫类型（T细胞，B细胞，巨噬细胞）：

只有一种类型的特异性免疫标记物（至少2个）；
有PTPRC markers和该类型的唯一特异性免疫标记物（至少1个）；
至少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降维聚类等等。哈哈哈，学无止境，我们继续结合作者的代码和文献内容，慢慢学吧.....

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

一.复杂热图

二.identifying cell types

相关推荐