NC单细胞文章复现(六)：Gene expression signatures(1) / 开普饭

在上一节，由于大部分细胞（868个）都被归为上皮细胞群中（Fig2 c），这868个细胞可被分成5个cluster，接着对这5个cluster细胞进行探索。我们使用一组来自对乳腺肿块的非监督分析的基因表达特征对5个cluster进行了研究。这些基因表达特征通过比较三阴性乳腺癌（TNBC）的四个亚型(ERBB2 amplicon，Luminal Subtype 、Basal epithelial-cell enriched 和Luminal epithelial gene cluster containing ER)而建立。先看看这5个clusters的basal细胞来源的细胞群有多少。大多数TNBC是基底样肿瘤，它们与多种TNBC型亚型重叠，与非固有基底TNBCs相比，与克隆异质性增加有关。（备注：这篇文献用到了很多apply循环，大家仔细琢磨，大概意思能看懂就行，然后可以把它应用到自己的数据中）

## 读取数据 basal_PNAS_all <- read.table("data/genes_for_basal_vs_non_basal_tnbc_PNAS.txt", header = TRUE, sep = "\t") #提取Basal.epithelial.cell.enriched.cluster的基因 basal_PNAS_long <- basal_PNAS_all$Basal.epithelial.cell.enriched.cluster #合并剩下17个基因 basal_PNAS <- intersect(basal_PNAS_long, rownames(mat_ct)) > basal_PNAS [1] "SOX9" "GALNT3" "CDH3" "LAMC2" "CX3CL1" "TRIM29" "KRT17" "KRT5" "CHI3L2" [10] "SLPI" "NFIB" "MRAS" "TGFB2" "CAPN6" "DMD" "FABP7" "CXCL1" #算出17个basal_PNAS基因在1112个细胞的表达平均值 basal_PNAS_avg_exprs <- apply(mat_ct[match(basal_PNAS, rownames(mat_ct)),], 2, mean) #检查一下数据 all.equal(names(basal_PNAS_avg_exprs), colnames(mat_ct)) #提取868个上皮细胞群体的17个basal_PNAS基因表达平均值 basal_PNAS_avg_exprs <- basal_PNAS_avg_exprs[which(pd_ct$cell_types_cl_all == "epithelial")]


#检查一下数据

all.equal(colnames(HSMM_allepith_clustering), names(basal_PNAS_avg_exprs))

#把17个basal_PNAS基因表达平均值赋给HSMM_allepith_clustering，以便于后续分析

pData(HSMM_allepith_clustering)$basal_PNAS_avg_exprs <- basal_PNAS_avg_exprs

#画figS9b

plot_cell_clusters(HSMM_allepith_clustering, 1, 2, color = "basal_PNAS_avg_exprs", cell_size = 2) + facet_wrap(~patient) +

  scale_color_continuous(low = "yellow", high = "blue")

#画figS9a plot_cell_clusters(HSMM_allepith_clustering, 1, 2, color = "basal_PNAS_avg_exprs", cell_size = 2) + facet_wrap(~Cluster) + scale_color_continuous(low = "yellow", high = "blue")

figS9a:大多数TNBC样本都有basal gene signature的表达。figS9b:在868个上皮细胞群中，cluster2的basal gene signature表达量最丰富。

接着，使用另外一个基因表达特征数据集TNBCtype4 signatures（Lehman_signature），这个signatures根据基因表达变化将TNBC细胞分为6个类：basal_like_1、basal_like_2、immunomodulatory、mesenchymal、mesenchymal_stem_like和luminal_ar。作者将基因表达特征中上调基因的平均表达值减去下调基因的平均表达值，将差值作为每个细胞在TNBCtype4 signatures （basal_like_1、basal_like_2、mesenchymal和luminal_ar）中的每个基因表达值，挑选最高基因表达值对应的signature，将其分配给对应的细胞。

#读取数据 lehman_long <- read.table("data/Lehman_signature.txt", sep = "\t", header = TRUE, stringsAsFactors = FALSE) #这个for循环提取了lehman_long里面的四列gene、regulation、no_samples和signature，建立一个data.rrame for (i in 0:5) { gene <- "gene" regulation <- "regulation" no_samples <- "no_samples" signature <- "signature" if (i == 0) { lehman <- lehman_long[, 1:4] lehman <- lehman[-which(lehman$signature == ""),] } if (i > 0) { gene <- paste("gene", i, sep = ".") regulation <- paste("regulation", i, sep = ".") no_samples <- paste("no_samples", i, sep = ".") signature <- paste("signature", i, sep = ".") mat_to_bind <- lehman_long[, c(gene, regulation, no_samples, signature)] colnames(mat_to_bind) <- c("gene", "regulation", "no_samples", "signature") if (length(which(is.na(mat_to_bind$no_samples))) > 0 ) mat_to_bind <- mat_to_bind[-which(mat_to_bind$signature == ""),] lehman <- rbind(lehman, mat_to_bind) } } #删掉一些mat_ct没有检测到的基因 lehman <- lehman[which(!is.na(match(lehman$gene, rownames(mat_ct)))),] lehman_signatures <- unique(lehman$signature) lehman_avg_exps <- apply(mat_ct, 2, function(x){ mns <- matrix(NA, nrow = length(lehman_signatures), ncol = 2) rownames(mns) <- lehman_signatures for (s in 1:length(lehman_signatures)) { sign <- lehman_signatures[s] # current signature lehman_here <- lehman %>% dplyr::filter(signature == sign) lehman_here_up <- lehman_here %>% dplyr::filter(regulation == "UP") lehman_here_down <- lehman_here %>% dplyr::filter(regulation == "DOWN") idx_genes_up <- match(lehman_here_up$gene, rownames(mat_ct)) idx_genes_down <- match(lehman_here_down$gene, rownames(mat_ct)) mns[s,] <- c(mean(x[idx_genes_up]), mean(x[idx_genes_down])) #算上调、下调的基因在样本中的平均表达值 } return(mns) }) #检查数据 all.equal(colnames(lehman_avg_exps), rownames(pd_ct)) #只看868个上皮细胞的情况 lehman_avg_exprs_epithelial <- lehman_avg_exps[,which(pd_ct$cell_types_cl_all == "epithelial")] #提取lehman_avg_exps前面6行，对应的是up lehman_avg_ups <- lehman_avg_exps[c(1:6), ] rownames(lehman_avg_ups) <- lehman_signatures all.equal(colnames(lehman_avg_ups), rownames(pd_ct)) lehman_avg_ups_epithelial <- lehman_avg_ups[,which(pd_ct$cell_types_cl_all == "epithelial")] #提取lehman_avg_exps后面6行，对应的是down lehman_avg_downs <- lehman_avg_exps[c(7:12),] rownames(lehman_avg_downs) <- lehman_signatures all.equal(colnames(lehman_avg_downs), rownames(pd_ct)) lehman_avg_downs_epithelial <- lehman_avg_downs[,which(pd_ct$cell_types_cl_all == "epithelial")] #上调基因的平均表达值减去下调基因的平均表达值 lehman_avg_both <- lehman_avg_ups - lehman_avg_downs all.equal(colnames(lehman_avg_both), rownames(pd_ct)) #挑选最高基因表达值对应的signature，将其分配给对应的细胞。 assignments_lehman_both <- apply(lehman_avg_both, 2, function(x){rownames(lehman_avg_both)[which.max(x)]}) assignments_lehman_both_epithelial <- assignments_lehman_both[which(pd_ct$cell_types_cl_all == "epithelial")] #删除immunomodulatory和mesenchymal_stem_like signature lehman_avg_both_epithelial_new <- lehman_avg_both_epithelial[-which(rownames(lehman_avg_both_epithelial) %in% c("immunomodulatory", "mesenchymal_stem_like")),] assignments_lehman_both_epithelial_new <- apply(lehman_avg_both_epithelial_new, 2, function(x){rownames(lehman_avg_both_epithelial_new)[which.max(x)]})

接下来画图，同样地，需要对heatmap函数代码进行修改。

ha_lehman_epith_pat <- list() for (i in 1:length(patients_now)) { if (i == 1) ha_lehman_epith_pat[[i]] <- HeatmapAnnotation(df=data.frame(cluster_all = clusterings_sep_allepith[[i]]), col = list(cluster_all = c("1" = "#ee204d", "2" = "#17806d", "3" = "#b2ec5d", "4" = "#cda4de", "5" = "#1974d2")), annotation_name_side = "left", annotation_name_gp = gpar(fontsize = 12), annotation_legend_param = list(list(title_position = "topcenter", title = "cluster")), show_annotation_name = FALSE, gap = unit(c(2), "mm"), show_legend = FALSE) if (i > 1 && i != 5 ) ha_lehman_epith_pat[[i]] <- HeatmapAnnotation(df=data.frame(cluster_all = clusterings_sep_allepith[[i]]), col = list(cluster_all = c("1" = "#ee204d", "2" = "#17806d", "3" = "#b2ec5d", "4" = "#cda4de", "5" = "#1974d2")), annotation_name_side = "left", annotation_name_gp = gpar(fontsize = 12), annotation_legend_param = list(list(title_position = "topcenter", title = "cluster")), show_annotation_name = FALSE, gap = unit(c(2), "mm"), show_legend = FALSE) if (i == 5) ha_lehman_epith_pat[[i]] <- HeatmapAnnotation(df=data.frame(cluster_all = clusterings_sep_allepith[[i]]), col = list(cluster_all = c("1" = "#ee204d", "2" = "#17806d", "3" = "#b2ec5d", "4" = "#cda4de", "5" = "#1974d2")), annotation_name_side = "right", annotation_name_gp = gpar(fontsize = 12), annotation_legend_param = list(list(title_position = "topcenter",title = "cluster")), show_annotation_name = FALSE, gap = unit(c(2), "mm"), show_legend = TRUE) } #检查数据 all.equal(names(lehmans_epith_pat_both), patients_now) #将basal signature添加进去，以便后续作图 lehmans_epith_pat_both_wbasal_new <- lehmans_epith_pat_both_new for (i in 1:length(patients_now)) { lehmans_epith_pat_both_wbasal_new[[i]] <- rbind(lehmans_epith_pat_both_new[[i]], pData(HSMM_allepith_clustering)$basal_PNAS_avg_exprs[which(HSMM_allepith_clustering$patient == patients_now[i])]) rownames(lehmans_epith_pat_both_wbasal_new[[i]])[5] <- "intrinsic_basal" }

# 画图 ht_sep_lehmans_both_wbasal_new <- Heatmap(lehmans_epith_pat_both_wbasal_new[[1]], col = colorRamp2(c(-0.7, 0, 1), c("blue","white", "red")), cluster_rows = FALSE, show_column_names = FALSE, column_title = patients_now[1], column_title_gp = gpar(fontsize = 12), top_annotation = ha_lehman_epith_pat[[1]], name = patients_now[1], show_row_names = FALSE, heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) + Heatmap(lehmans_epith_pat_both_wbasal_new[[2]], col = colorRamp2(c(-0.7, 0, 1), c("blue","white", "red")), cluster_rows = FALSE, show_column_names = FALSE, column_title = patients_now[2], column_title_gp = gpar(fontsize = 12), top_annotation = ha_lehman_epith_pat[[2]], name = patients_now[2], show_row_names = FALSE, heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) + Heatmap(lehmans_epith_pat_both_wbasal_new[[3]], col = colorRamp2(c(-0.7, 0, 1), c("blue","white", "red")), cluster_rows = FALSE, show_column_names = FALSE, column_title = patients_now[3], column_title_gp = gpar(fontsize = 12), top_annotation = ha_lehman_epith_pat[[3]], name = patients_now[3], show_row_names = FALSE, heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) + Heatmap(lehmans_epith_pat_both_wbasal_new[[4]], col = colorRamp2(c(-0.7, 0, 1), c("blue","white", "red")), cluster_rows = FALSE, show_column_names = FALSE, column_title = patients_now[4], column_title_gp = gpar(fontsize = 12), top_annotation = ha_lehman_epith_pat[[4]], name = patients_now[4], show_row_names = FALSE, heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) + Heatmap(lehmans_epith_pat_both_wbasal_new[[5]], col = colorRamp2(c(-0.7, 0, 1), c("blue","white", "red")), cluster_rows = FALSE, show_column_names = FALSE, column_title = patients_now[5], column_title_gp = gpar(fontsize = 12), top_annotation = ha_lehman_epith_pat[[5]], name = patients_now[5], show_row_names = FALSE, heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) + Heatmap(lehmans_epith_pat_both_wbasal_new[[6]], col = colorRamp2(c(-0.7, 0, 1), c("blue","white", "red")), cluster_rows = FALSE, row_names_side = "right", column_title = patients_now[6], column_title_gp = gpar(fontsize = 12), top_annotation = ha_lehman_epith_pat[[6]], name = patients_now[6], show_column_names = FALSE, heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) #画fig3d print(draw(ht_sep_lehmans_both_wbasal_new, annotation_legend_side = "right"))

我们只需要把右边注释条PS一下，就可以达到跟文献的图片一模一样了。

#检查数据 all.equal(colnames(HSMM_allepith_clustering), names(assignments_lehman_both_epithelial_new)) pData(HSMM_allepith_clustering)$assignments_lehman_both_new <- assignments_lehman_both_epithelial_new 画fig3g plot_cell_clusters(HSMM_allepith_clustering, 1, 2, color = "assignments_lehman_both_new", cell_size = 2) + facet_wrap(~patient)

#画figS8 plot_cell_clusters(HSMM_allepith_clustering, 1, 2, color = "assignments_lehman_both_new", cell_size = 2) + facet_wrap(~Cluster)

cluster4也富集了 Basal-Like 1 signature，而cluster3高度富集了TNBCtype 中的“Luminal Androgen Receptor” signature。为了更清楚的看到上皮细胞群的5个cluster对应的TNBCtype signatures的平均表达量，接着继续探索下去...

clust_avg_lehman_both_new <- matrix(NA, nrow = length(unique(HSMM_allepith_clustering$Cluster)), ncol = nrow(lehman_avg_both_epithelial_new)) #列名：cluster1,cluster2,cluster3,cluster4,cluster5 rownames(clust_avg_lehman_both_new) <- paste("clust", c(1:length(unique(HSMM_allepith_clustering$Cluster))), sep = "") #行名：basal_like_1、basal_like_2、mesenchymal、luminal_ar" colnames(clust_avg_lehman_both_new) <- rownames(lehman_avg_both_epithelial_new) #算出每个cluster的signatures 平均值 for (c in 1:length(unique(HSMM_allepith_clustering$Cluster))) { clust_avg_lehman_both_new[c,] <- apply(lehman_avg_both_epithelial_new[,which(HSMM_allepith_clustering$Cluster == c)], 1, mean) }


clust_avg_lehman_both_new <- as.data.frame(clust_avg_lehman_both_new)

#增加一列cluster

clust_avg_lehman_both_new$Cluster <- rownames(clust_avg_lehman_both_new)

#拆分数据

clust_avg_lehman_melt_new <- melt(clust_avg_lehman_both_new, "Cluster")

#画fig3e ggplot(clust_avg_lehman_melt_new, aes(Cluster, value, fill = factor(variable), color = factor(variable), shape = factor(variable))) + geom_point(size = 3, stroke = 1) + scale_shape_discrete(solid = T) + #guides(colour = guide_legend(override.aes = list(size=3))) + ylab("average expression of signature in cluster") + xlab("cluster") + ylim(c(-0.35, 0.5))

可以看到：cluster2和4强烈地表达Basal-Like 1 signature，而cluster3显著表达Basal-Like 2 signature和 luminal_AR signature

接着，读取另外一个ML_signature(mature luminal signature)，将上调基因的平均表达量中减去下调基因的平均表达量，计算出three normal breast signatures下每个细胞的表达量(Lim et al., 2009a)，然后将每个细胞分配给其具有最高表达量的signatures。这三个normal breast signatures 是：mature luminal (ML)，basal和luminal progenitor (LP)，在每个signatures中，都有对应的上调基因和下调基因。

ml_signature_long <- read.table("data/ML_signature.txt", sep = "\t", header = TRUE) if (length(which(ml_signature_long$Symbol == "")) > 0) #将空白行去掉 ml_signature_long <- ml_signature_long[-which(ml_signature_long$Symbol == ""),] #按照基因字母进行排序，如果基因字母有一样的，则按照Average.log.fold.change绝对值的负数进行从小到大排序 ml_signature_long <- ml_signature_long[order(ml_signature_long$Symbol, -abs(ml_signature_long$Average.log.fold.change) ), ] #对基因取唯一值，去重复 ml_signature_long <- ml_signature_long[ !duplicated(ml_signature_long$Symbol), ] #总共有818个基因 ml_signature <- ml_signature_long[which(!is.na(match(ml_signature_long$Symbol, rownames(mat_ct)))), ] #上调基因有384个 ml_up <- ml_signature[which(ml_signature$Average.log.fold.change > 0), ] #下调基因有197个 ml_down <- ml_signature[which(ml_signature$Average.log.fold.change < 0), ] #匹配一下 idx_ml_up <- match(ml_up$Symbol, rownames(mat_ct)) idx_ml_down <- match(ml_down$Symbol, rownames(mat_ct)) #读取basal signature，处理过程跟上面的一样的。 basal_signature_long <- read.table("data/basal_signature.txt", sep = "\t", header = TRUE) if (length(which(basal_signature_long$Symbol == "")) > 0) basal_signature_long <- basal_signature_long[-which(basal_signature_long$Symbol == ""),] basal_signature_long <- basal_signature_long[order(basal_signature_long$Symbol, -abs(basal_signature_long$Average.log.fold.change) ), ] basal_signature_long <- basal_signature_long[ !duplicated(basal_signature_long$Symbol), ] #总共有1335个基因 basal_signature <- basal_signature_long[which(!is.na(match(basal_signature_long$Symbol, rownames(mat_ct)))), ] #上调基因有588个 basal_up <- basal_signature[which(basal_signature$Average.log.fold.change > 0), ] #下调基因有757个 basal_down <- basal_signature[which(basal_signature$Average.log.fold.change < 0), ] idx_basal_up <- match(basal_up$Symbol, rownames(mat_ct)) idx_basal_down <- match(basal_down$Symbol, rownames(mat_ct))


#读取LP signature,还是同样的操作

lp_signature_long <- read.table("data/Lp_signature.txt", sep = "\t", header = TRUE)

if (length(which(lp_signature_long$Symbol == "")) > 0)

  lp_signature_long <- lp_signature_long[-which(lp_signature_long$Symbol == ""),]

lp_signature_long <- lp_signature_long[order(lp_signature_long$Symbol, -abs(lp_signature_long$Average.log.fold.change) ), ]

lp_signature_long <- lp_signature_long[ !duplicated(lp_signature_long$Symbol), ]

lp_signature <- lp_signature_long[which(!is.na(match(lp_signature_long$Symbol, rownames(mat_ct)))), ]

lp_up <- lp_signature[which(lp_signature$Average.log.fold.change > 0), ]

lp_down <- lp_signature[which(lp_signature$Average.log.fold.change < 0), ]

idx_lp_up <- match(lp_up$Symbol, rownames(mat_ct))

idx_lp_down <- match(lp_down$Symbol, rownames(mat_ct))

#对ML、basal和LP 3个signatures基因，将上调基因的表达值减去下调基因表达值，并分别返回结果。

normsig_avg_exprs <- apply(mat_ct, 2, function(x){

  

  avg_ml_up <- mean(x[idx_ml_up])

  avg_ml_down <- mean(x[idx_ml_down])

  avg_ml_both <- avg_ml_up - avg_ml_down

  

  avg_basal_up <- mean(x[idx_basal_up])

  avg_basal_down <- mean(x[idx_basal_down])

  avg_basal_both <- avg_basal_up - avg_basal_down

  

  avg_lp_up <- mean(x[idx_lp_up])

  avg_lp_down <- mean(x[idx_lp_down])

  avg_lp_both <- avg_lp_up - avg_lp_down

  

  return(c(avg_ml_up, avg_basal_up, avg_lp_up, avg_ml_both, avg_basal_both, avg_lp_both))

})

rownames(normsig_avg_exprs) <- c("avg_ml_up", "avg_basal_up", "avg_lp_up", "avg_ml_both", "avg_basal_both", "avg_lp_both")

#检查数据

all.equal(colnames(normsig_avg_exprs), rownames(pd_ct))

#只看上皮细胞群

normsig_avg_exprs_epithelial <- normsig_avg_exprs[,which(pd_ct$cell_types_cl_all == "epithelial")]
normsig_avg_ups <- normsig_avg_exprs[c(1:3), ]

all.equal(colnames(normsig_avg_ups), rownames(pd_ct))

normsig_avg_ups_epithelial <- normsig_avg_ups[,which(pd_ct$cell_types_cl_all == "epithelial")]
normsig_avg_both <- normsig_avg_exprs[c(4:6),]

all.equal(colnames(normsig_avg_both), rownames(pd_ct))

normsig_avg_both_epithelial <- normsig_avg_both[,which(pd_ct$cell_types_cl_all == "epithelial")]

#挑选最大值=上调基因的平均表达值最大数值分配给对应的细胞类型

assignments_normsig_ups <- apply(normsig_avg_ups, 2, function(x){rownames(normsig_avg_ups)[which.max(x)]})

assignments_normsig_ups_epithelial <- assignments_normsig_ups[which(pd_ct$cell_types_cl_all == "epithelial")]

#上调基因的平均表达值-下调基因的平均表达值的最大数值分配给对应的细胞类型

assignments_normsig_both <- apply(normsig_avg_both, 2, function(x){rownames(normsig_avg_both)[which.max(x)]})

assignments_normsig_both_epithelial <- assignments_normsig_both[which(pd_ct$cell_types_cl_all == "epithelial")]
# heatmaps on normal signatures per patient

pd_ct_epith <- pd_ct[which(pd_ct$cell_types_cl_all == "epithelial"),]

normsig_epith_pat_both <- list()

normsig_epith_pat_ups <- list()

pds_epith_ct <- list()

for (i in 1:length(patients_now)) {

  normsig_epith_pat_both[[i]] <- normsig_avg_both_epithelial[,which(pd_ct_epith$patient == patients_now[i])]

  normsig_epith_pat_ups[[i]] <- normsig_avg_ups_epithelial[,which(pd_ct_epith$patient == patients_now[i])]

  pds_epith_ct[[i]] <- pds_ct[[i]][which(pds_ct[[i]]$cell_types_cl_all == "epithelial"),]

}

names(normsig_epith_pat_both) <- patients_now

names(normsig_epith_pat_ups) <- patients_now

names(pds_epith_ct) <- patients_now
ht_sep_normsig_both <-

  Heatmap(normsig_epith_pat_both[[1]],

          col = colorRamp2(c(-0.7, -0.2, 0.7), c("blue","white", "red")),

          cluster_rows = FALSE,

          show_column_names = FALSE,

          column_title = patients_now[1],

          top_annotation = ha_lehman_epith_pat[[1]],

          column_title_gp = gpar(fontsize = 12),

          show_row_names = FALSE,

          name = patients_now[1], 

        

          heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) +

  Heatmap(normsig_epith_pat_both[[2]],

          col = colorRamp2(c(-0.7, -0.2, 0.7), c("blue","white", "red")),

          cluster_rows = FALSE,

          show_column_names = FALSE,

          column_title = patients_now[2],

          column_title_gp = gpar(fontsize = 12),

          top_annotation = ha_lehman_epith_pat[[2]],

          name = patients_now[2], 

          show_row_names = FALSE,

       

          heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) +

  Heatmap(normsig_epith_pat_both[[3]],

          col = colorRamp2(c(-0.7, -0.2, 0.7), c("blue","white", "red")),

          cluster_rows = FALSE,

          show_column_names = FALSE,

          column_title = patients_now[3],

          column_title_gp = gpar(fontsize = 12),

          top_annotation = ha_lehman_epith_pat[[3]],

          name = patients_now[3], 

          show_row_names = FALSE,

          top_annotation_height = unit(c(2), "cm"),

          heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) +

  Heatmap(normsig_epith_pat_both[[4]],

          col = colorRamp2(c(-0.7, -0.2, 0.7), c("blue","white", "red")),

          cluster_rows = FALSE,

          show_column_names = FALSE,

          column_title = patients_now[4],

          column_title_gp = gpar(fontsize = 12),

          top_annotation = ha_lehman_epith_pat[[4]],

          name = patients_now[4], 

          show_row_names = FALSE,

          top_annotation_height = unit(c(2), "cm"),

          heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) +

  Heatmap(normsig_epith_pat_both[[5]],

          col = colorRamp2(c(-0.7, -0.2, 0.7), c("blue","white", "red")),

          cluster_rows = FALSE,

          show_column_names = FALSE,

          column_title = patients_now[5],

          column_title_gp = gpar(fontsize = 12),

          top_annotation = ha_lehman_epith_pat[[5]],

          name = patients_now[5], 

          show_row_names = FALSE,

          top_annotation_height = unit(c(2), "cm"),

          heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) +

  Heatmap(normsig_epith_pat_both[[6]],

          col = colorRamp2(c(-0.7, -0.2, 0.7), c("blue","white", "red")),

          cluster_rows = FALSE,

          row_names_side = "right",

          column_title = patients_now[6],

          column_title_gp = gpar(fontsize = 12),

          top_annotation = ha_lehman_epith_pat[[6]],

          name = patients_now[6], 

          show_column_names = FALSE,

        

          top_annotation_height = unit(c(2), "cm"),

          heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9)))

#画fig3b.pdf

print(draw(ht_sep_normsig_both, annotation_legend_side = "bottom"))

# 每个样本的normal signatures 数目 all.equal(colnames(HSMM_allepith_clustering), names(assignments_normsig_both_epithelial)) pData(HSMM_allepith_clustering)$assignments_normsig_both <- assignments_normsig_both_epithelial pData(HSMM_allepith_clustering)$assignments_normsig_ups <- assignments_normsig_ups_epithelial

#画fig3f plot_cell_clusters(HSMM_allepith_clustering, 1, 2, color = "assignments_normsig_both", cell_size = 2) + facet_wrap(~patient)

#每个clusters 的normal signatures 数目 plot_cell_clusters(HSMM_allepith_clustering, 1, 2, color = "assignments_normsig_both", cell_size = 2) + facet_wrap(~Cluster)

#检查数据 all.equal(HSMM_allepith_clustering$Cluster, clustering_allepith) all.equal(colnames(normsig_avg_both_epithelial), colnames(HSMM_allepith_clustering)) clust_avg_normsig_both <- matrix(NA, nrow = length(unique(HSMM_allepith_clustering$Cluster)), ncol = nrow(normsig_avg_both_epithelial)) rownames(clust_avg_normsig_both) <- paste("clust", c(1:length(unique(HSMM_allepith_clustering$Cluster))), sep = "") colnames(clust_avg_normsig_both) <- rownames(normsig_avg_both_epithelial) #算上皮细胞群的avg_both 平均表达值，接下来还是同样的操作 for (c in 1:length(unique(HSMM_allepith_clustering$Cluster))) { clust_avg_normsig_both[c,] <- apply(normsig_avg_both_epithelial[,which(HSMM_allepith_clustering$Cluster == c)], 1, mean) }

clust_avg_normsig_both <- as.data.frame(clust_avg_normsig_both) clust_avg_normsig_both$Cluster <- rownames(clust_avg_normsig_both) clust_avg_normsig_melt <- melt(clust_avg_normsig_both, "Cluster") #画fig3e ggplot(clust_avg_normsig_melt, aes(Cluster, value, fill = factor(variable), color = factor(variable), shape = factor(variable))) + geom_point(size = 3, stroke = 1) + scale_shape_discrete(solid = T) + ylab("average expression of signature in cluster") + xlab("cluster") + ylim(c(-0.35, 0.5))

fig3e：Clusters 2和Clusters4 强烈表达 LP signature, 而cluster 3则高表达 ML signature.

接着为了进一步探究临床相关性，作者使用了三个临床相关的gene signatures，进一步探究这868个上皮细胞的特征，这868个上皮细胞真的被研究到很彻底，真的佩服，这工作量好大....

第一个gene signatures:70-gene prognostic signature ，该signatures最初是从对有无转移复发患者的原发肿瘤之间差异表达基因的分析中得出的，总共70个基因。

mammaprint_long <- read.table("data/mammaprint_sig_new.txt", header = TRUE, sep = "\t") mammaprint <- apply(mammaprint_long, 2, function(x){return(intersect(x, rownames(mat_ct)))})[,1] mammaprint_avg_exprs <- apply(mat_ct[match(mammaprint, rownames(mat_ct)),], 2, mean) all.equal(names(mammaprint_avg_exprs), colnames(mat_ct)) mammaprint_avg_exprs <- mammaprint_avg_exprs[which(pd_ct$cell_types_cl_all == "epithelial")]


all.equal(colnames(HSMM_allepith_clustering), names(mammaprint_avg_exprs))

pData(HSMM_allepith_clustering)$mammaprint_avg_exprs <- mammaprint_avg_exprs

# 画figS13b plot_cell_clusters(HSMM_allepith_clustering, 1, 2, color = "mammaprint_avg_exprs", cell_size = 2) + facet_wrap(~patient) + scale_color_continuous(low = "yellow", high = "blue")

# 画figS13a plot_cell_clusters(HSMM_allepith_clustering, 1, 2, color = "mammaprint_avg_exprs", cell_size = 2) + facet_wrap(~Cluster) + scale_color_continuous(low = "yellow", high = "blue")

第二个gene signatures:49-gene metastatic burden signature.该signatures可以区分了患者来源的小鼠TNBC异种移植模型中单个循环转移细胞所产生的高转移负荷和低转移负荷，总共包括49个基因。

zenawerb_long <- read.table("data/werb_49_metastasis_sig.txt", header = TRUE, sep = "\t") zenawerb <- apply(zenawerb_long, 2, function(x){return(intersect(x, rownames(mat_ct)))})[,1] zenawerb_avg_exprs <- apply(mat_ct[match(zenawerb, rownames(mat_ct)),], 2, mean) all.equal(names(zenawerb_avg_exprs), colnames(mat_ct)) zenawerb_avg_exprs <- zenawerb_avg_exprs[which(pd_ct$cell_types_cl_all == "epithelial")]


all.equal(colnames(HSMM_allepith_clustering), names(zenawerb_avg_exprs))

pData(HSMM_allepith_clustering)$zenawerb_avg_exprs <- zenawerb_avg_exprs

#画figS14b plot_cell_clusters(HSMM_allepith_clustering, 1, 2, color = "zenawerb_avg_exprs", cell_size = 2) + facet_wrap(~patient) + scale_color_continuous(low = "yellow", high = "blue")

#画figS14a plot_cell_clusters(HSMM_allepith_clustering, 1, 2, color = "zenawerb_avg_exprs", cell_size = 2) + facet_wrap(~Cluster) + scale_color_continuous(low = "yellow", high = "blue")

第三个gene siganatures:从接受手术前化疗治疗的原发性乳腺癌患者的残存肿瘤群体中富集的基因中获得的，包括354个基因。

artega_long <- read.table("data/artega_sig.txt", header = TRUE, sep = "\t") artega <- apply(artega_long, 2, function(x){return(intersect(x, rownames(mat_ct)))})[,1] artega_avg_exprs <- apply(mat_ct[match(artega, rownames(mat_ct)),], 2, mean) all.equal(names(artega_avg_exprs), colnames(mat_ct)) artega_avg_exprs <- artega_avg_exprs[which(pd_ct$cell_types_cl_all == "epithelial")]


all.equal(colnames(HSMM_allepith_clustering), names(artega_avg_exprs))

pData(HSMM_allepith_clustering)$artega_avg_exprs <- artega_avg_exprs

画figS15a plot_cell_clusters(HSMM_allepith_clustering, 1, 2, color = "artega_avg_exprs", cell_size = 2) + facet_wrap(~patient) + scale_color_continuous(low = "yellow", high = "blue")

#画figS15b plot_cell_clusters(HSMM_allepith_clustering, 1, 2, color = "artega_avg_exprs", cell_size = 2) + facet_wrap(~Cluster) + scale_color_continuous(low = "yellow", high = "blue")

#将3个gene signatures的表达值合成一个数据框 prognosis_sig <- cbind(mammaprint_avg_exprs, zenawerb_avg_exprs, artega_avg_exprs) #取行名 colnames(prognosis_sig) <- c("mammaprint", "zenawerb", "artega")

prognosis_epith_pat <- list() for (i in 1:length(patients_now)) { prognosis_epith_pat[[i]] <- t(prognosis_sig)[,which(pd_ct_epith$patient == patients_now[i])] } names(prognosis_epith_pat) <- patients_now for (i in 1:length(patients_now)) { print(all.equal(colnames(prognosis_epith_pat[[1]]), rownames(pds_epith_ct[[1]]))) print(all.equal(names(clusterings_sep_allepith[[1]]), colnames(prognosis_epith_pat[[1]]))) } ht_sep_prognosis <- Heatmap(prognosis_epith_pat[[1]], cluster_rows = FALSE, col = colorRamp2(c(-0.2, 0.2, 1), c("blue","white", "red")), show_column_names = FALSE, column_title = patients_now[1], top_annotation = ha_lehman_epith_pat[[1]], column_title_gp = gpar(fontsize = 12), show_row_names = FALSE, name = patients_now[1], heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) + Heatmap(prognosis_epith_pat[[2]], col = colorRamp2(c(-0.2, 0.2, 1), c("blue","white", "red")), cluster_rows = FALSE, show_column_names = FALSE, column_title = patients_now[2], column_title_gp = gpar(fontsize = 12), top_annotation = ha_lehman_epith_pat[[2]], name = patients_now[2], show_row_names = FALSE, heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) + Heatmap(prognosis_epith_pat[[3]], col = colorRamp2(c(-0.2, 0.2, 1), c("blue","white", "red")), cluster_rows = FALSE, show_column_names = FALSE, column_title = patients_now[3], column_title_gp = gpar(fontsize = 12), top_annotation = ha_lehman_epith_pat[[3]], name = patients_now[3], show_row_names = FALSE, heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) + Heatmap(prognosis_epith_pat[[4]], col = colorRamp2(c(-0.2, 0.2, 1), c("blue","white", "red")), cluster_rows = FALSE, show_column_names = FALSE, column_title = patients_now[4], column_title_gp = gpar(fontsize = 12), top_annotation = ha_lehman_epith_pat[[4]], name = patients_now[4], show_row_names = FALSE, heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) + Heatmap(prognosis_epith_pat[[5]], col = colorRamp2(c(-0.2, 0.2, 1), c("blue","white", "red")), cluster_rows = FALSE, show_column_names = FALSE, column_title = patients_now[5], column_title_gp = gpar(fontsize = 12), top_annotation = ha_lehman_epith_pat[[5]], name = patients_now[5], show_row_names = FALSE, heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) + Heatmap(prognosis_epith_pat[[6]], col = colorRamp2(c(-0.2, 0.2, 1), c("blue","white", "red")), cluster_rows = FALSE, row_names_side = "right", column_title = patients_now[6], column_title_gp = gpar(fontsize = 12), top_annotation = ha_lehman_epith_pat[[6]], name = patients_now[6], show_column_names = FALSE, heatmap_legend_param = list(title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9))) #画fig4a print(draw(ht_sep_prognosis, annotation_legend_side = "right"))

NC单细胞文章复现(六)：Gene expression signatures(1)

相关推荐