我們這次直接拿GSE118390上已經(jīng)normalized 的數(shù)據(jù)進行下游分析。首先我們先看看文獻的這張復(fù)雜熱圖褒颈，哈哈寨典，這張熱圖畫得真是好看。左邊是不同的markers基因?qū)?yīng)的細胞類型认然，上邊是6個TNBC病人總共1189個細胞补憾，cycling指的是處于細胞分裂期的樣本，下邊還有1189個細胞對應(yīng)的是CD45陽性還是CD陰性的卷员。代碼作者都已經(jīng)有提供盈匾，但是如果按照作者的代碼跑，仍然會出現(xiàn)error毕骡，我把代碼一些參數(shù)修改了一下并進行翻譯削饵，再結(jié)合文獻的描述，跟大家一起學(xué)習(xí)未巫。我們看看怎么畫唄窿撬。
Fig1.b
作者首先將已建立的用于定義免疫細胞、內(nèi)皮細胞和基質(zhì)細胞以及關(guān)鍵乳腺上皮亞群（basal cells, luminal progenitors, andmature luminal cells）的特定基因集去注釋1189個細胞的類型叙凡，然后使用聚類將1189個細胞中的1112個分為非上皮（non-epithelial：n=244）和上皮（epithelial：n=868）兩種類型尤仍。已知TNBC的一個亞類會表現(xiàn)出大量的免疫細胞浸潤，CD45陰性細胞（CD45-unselected）中含有大量的免疫細胞亞群狭姨，其中大部分是T淋巴細胞，其次是巨噬細胞苏遥。間質(zhì)細胞(stromal cells)成分在TNBC病例之間也存在差異饼拍，這種細胞在每個腫瘤中占總細胞的比例<15%。內(nèi)皮細胞（endothelial cells）是少數(shù)群體田炭，最多占腫瘤細胞的4%师抄。
#加載包
library(here)
source(here::here("code", "1.libraries.R"))
#funcs.R和funcs_markers是已經(jīng)包裝好的函數(shù)，我們接下來會有稍微講解一下教硫。
source(here::here("code", "funcs.R"))
source(here::here("code", "funcs_markers.R"))
#讀取數(shù)據(jù)和表型信息
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)
#創(chuàng)建SingleCellExperiment對象
sceset_final <- SingleCellExperiment(assays = list(exprs = as.matrix(mat_norm)),
colData = pd_norm,
rowData = fd_norm)
#對各個指標(biāo)創(chuàng)建創(chuàng)建圖像矢量
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函數(shù)是包裝在funs.R中瞬矩，
#具體來說是先刪除#melanoma_g1s的重復(fù)基因茶鉴，刪除空白字符，轉(zhuǎn)換成data.frame景用，匹配#mat_norm和 melanoma_g1s #marker基因涵叮，再判斷一下是否有缺失值。
melanoma_g1s <- match_clean_vector_genes(mat_norm, melanoma_g1s)
#avg_expr_genes是用apply包裝的一個函數(shù)伞插，即算出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期評分≥其中位數(shù)+2MAD和細胞的G2M期評分≤其中位數(shù)+2 MAD媚污。2舀瓢、細胞的G1S期評分≤其中位數(shù)2 MADs 和細胞的G2M期評分≥2 MADs。3耗美、細胞的G1S期評分≥其中位數(shù)+2倍的絕對中位差和細胞的G2M期評分≥其中位數(shù)+2倍的絕對中位差京髓。
non-cycling cells：細胞的G1S期評分<其中位數(shù)+2MAD和細胞的G2M期評分<其中位數(shù)+2 MAD
#先創(chuàng)建一個1189列的數(shù)據(jù)
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個诚些，難道中間哪一步有瑕疵飞傀？別急，我們慢慢往下看诬烹。
Identifying cycling cells
#把cycling和non-cycling的信息添加到sceset_final中,方便后續(xù)處理
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按照字母和數(shù)值從小到大排序
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個樣本各自的表型信息
接著這篇文章作者開始進行細胞注釋砸烦，分為兩個步驟。
1.使用于文獻的特定表達標(biāo)記的基因列表來定義細胞類型绞吁。
2.clustering
Identifying cell types
對于第1步驟幢痘，使用的基因列表是來自(Tirosh et al., 2016）整理好的用來注釋四種細胞類型（epithelial，stroma家破，endothelial颜说，immune）的53個marker基因。
#我們讀取已經(jīng)整理好的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個基因里面
#若沒有的話门粪，將其刪除。經(jīng)過這一步過濾烹困，將原來53個marker玄妈，過濾成49個marker。
markers <- unique(all_markers[all_markers$gene %in% rowData(sceset_final)$hgnc_symbol, ])
length(markers$gene)
> length(markers$gene) #剩下49個markers
[1] 49
接下是畫熱圖髓梅，這張圖的信息量很大拟蜻，我們慢慢來品。
#首先是先將熱圖的注釋信息弄好枯饿，這一步驟是將49個marker基因在熱圖中所對應(yīng)的細胞類型表弄好酝锅。
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循環(huán)是用一個replace函數(shù)奢方，將49個markers的細胞類型附上對應(yīng)的顏
#顏色的注釋信息在我們之前的創(chuàng)建圖像矢量的時候就已經(jīng)建立好
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])
}
##調(diào)整將熱圖左邊的細胞類型注釋條順序屈张，從上到下依次為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)])
##再調(diào)整將熱圖下邊的細胞類型注釋條順序
#順序依次為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)])
##注釋條的顏色璧尸，大小，標(biāo)題等熬拒。
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"))
接下來這一步是將復(fù)雜熱圖各自注釋信息例如聚類方向爷光、字體位置、題目澎粟、邊框等信息整理好蛀序，但是如果按照作者的代碼跑，會出現(xiàn)莫名其妙的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.
對于剛學(xué)單細胞測序的我來說活烙，這個問題的解決的過程蠻“崎嶇"的徐裸，得耗費時間不斷去嘗試各種方法。好在最后順利解決啸盏。
1.首先我們安裝最新版本的ComplexHeatmap package重贺，如果你只是用bioconductor的官方代碼去安裝，沒法安裝到最新版本的ComplexHeatmap package回懦，解決方法可以是利用devtools package去github上安裝气笙。
2.HeatmapAnnotation函數(shù)有些參數(shù)需要修改，我把修改好的貼上來了怯晕。至于具體的畫圖參數(shù)潜圃，我就不細講了，有興趣的小伙伴可以去看一看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"))), #熱圖下邊關(guān)于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"))
}
接著可以開始正式畫圖了
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"))
#pdf(("Fig1.b.pdf"), onefile = FALSE, width = 10, height = 15)
draw(ht_list, gap = unit(0.1, "cm"), heatmap_legend_side = "bottom", annotation_legend_side = "bottom")
#dev.off()
Fig1.b
這樣子一張復(fù)雜熱圖總算完成了蛉谜，基本跟文獻是一樣的稚晚。但是我們的學(xué)習(xí)還沒有結(jié)束。
接著需要定義markers的表達量大于多少才算是有表達型诚，設(shè)置的閾值是1.
Expression markers
接著需要確定這些細胞是否是：
1.上皮性來源細胞客燕，需滿足一下2個條件之一：1.至少2個上皮性markers,2.若該細胞只有一個乳腺上皮標(biāo)志物：EpCAM、KRT8狰贯、KRT18也搓、KRT19，對于該標(biāo)志物涵紊，該標(biāo)志物在該患者中須有超過50%的細胞有表達傍妒。
2.免疫成分細胞:如果一個細胞表現(xiàn)出特殊的免疫類型（T細胞，B細胞摸柄，巨噬細胞）：1.只有一種類型的特異性免疫標(biāo)記物（至少2個）颤练；2.有PTPRC markers和該類型的唯一特異性免疫標(biāo)記物（至少1個）；3.至少3個這種類型的免疫標(biāo)記驱负，最多1個其他類型的免疫標(biāo)記嗦玖。
比較特殊的是若該細胞有PTRC患雇，則再看有沒有別的marker，若還有1個以上的免疫細胞宇挫，
那么需要繼續(xù)細分免疫細胞亞群（macrophage苛吱，T cell，B cell）器瘪，若沒有其他marker翠储，則籠統(tǒng)地定義為免疫細胞。
3.基質(zhì)成分細胞：該細胞僅有基質(zhì)標(biāo)志物娱局；該細胞有最多3個基質(zhì)標(biāo)志物彰亥，且至少有1個內(nèi)皮標(biāo)志物。
4.內(nèi)皮成分細胞：該細胞僅有內(nèi)皮標(biāo)志物衰齐；或者該細胞有至少3個內(nèi)皮標(biāo)志物和最多1個基質(zhì)標(biāo)志物任斋。
具體的話可以參考
image.png
接下來的lists_markers、decide_is_epithelial 和decide_is_immune這3個函數(shù)是包裝在了funcs_markers.R里面耻涛，都是用for循環(huán)實現(xiàn)的伟端。簡單來說：lists_markers函數(shù)先把細胞分為epithelial_cells、immune_cells和other_cells终佛，然后再進行細分趟径。decide_is_epithelial函數(shù)進一步確定epithelial_cells是否有2個上皮性markers，如果有卓研，將該細胞標(biāo)記為1趴俘，如果沒有，將該細胞標(biāo)記為0奏赘。decide_is_immune函數(shù)比較復(fù)雜寥闪，如果該細胞沒有immune markers或者僅僅只有一個markers，標(biāo)記為0磨淌；如果該細胞至少有2個immune markers疲憋，并且沒有PTPRC markers，那么返回該markers對應(yīng)的細胞類型梁只；如果該細胞至少有1種類型的2個immune markers,并且有PTPRC markers缚柳，那么返回該markers對應(yīng)的細胞類型
如果該細胞有1種細胞類型的3個immune markers，而且還含有另外1種類型的1個markers搪锣，則返回有3個immune markers的類型秋忙。如果大家看不懂，可結(jié)合上面那張圖和代碼去看构舟。
thresh <- 1 #表達閾值設(shè)置為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)
#作者在這一步里繼續(xù)細化注釋上皮細胞的流程翰绊，如果這個細胞只表達1個epithelial marker，
#先查看這個markers在樣本中的表達分布,如果分布在超過50%的細胞，則標(biāo)記為1监嗜，若沒有標(biāo)記為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標(biāo)記為0
is_immune_simple <- is_immune
is_immune_simple[which(is_immune == "immune_mix")] <- 0
#將我們之前分好的有多種細胞類型的other_cell標(biāo)記為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)
#若出現(xiàn)0-0-0谐檀，則該細胞是unknown，若出現(xiàn)epithelial-stroma-0,該細胞定義為epithelial cells裁奇，理由是：
#同時具備epithelial和stroma markers的細胞一般跟EMT有關(guān)桐猬，可歸入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"
這里面for循環(huán)很多溃肪，我也只是看懂個大概意思，可能寫得也不怎么清楚音五，具體里面的內(nèi)在邏輯還需要結(jié)合文章的方法去不斷理解消化惫撰。不過我們也是學(xué)習(xí)到了一些細胞注釋的知識，例如在注釋上皮細胞類型的時候躺涝，需要注意EPCAM, KRT8, KRT18,和KRT19這些基因厨钻，這四個基因在上皮細胞高度表達，可單獨拿出來作為上皮細胞標(biāo)志坚嗜。PTPRC也稱CD45夯膀，最初被稱為白細胞共同抗原，可作為免疫細胞亞群的第一次分群的標(biāo)志苍蔬。若細胞同時具備epithelial和stroma兩種markers诱建，則應(yīng)該被歸為上皮細胞亞群。
我們下一篇會開始學(xué)習(xí)作者怎么又把1189個細胞減到1112個有“名字”的細胞碟绑，然后開始t-SNE降維聚類等等俺猿，哈哈哈，對于我來說格仲，這知識量真的是爆炸押袍，學(xué)無止境，結(jié)合作者的代碼和文獻內(nèi)容抓狭，慢慢學(xué)吧.....
