单细胞RNA-seq揭示TNBC的异质性（图表复现04)

 后起之秀奔涌而至，欢迎大家在《生信技能树》的舞台分享自己的心得体会！（文末有惊喜）

下面是《单细胞天地》小编日行一膳的投稿

前面的单细胞RNA-seq揭示TNBC的异质性（图表复现03）教程里面我们一起复现了文章“ Unravelling subclonal heterogeneity and aggressive disease states in TNBC through single-cell RNA-seq  ” 的Figure 2 ，这周继续来复现一下Figure 3的相关内容

复现图表

Fig4a

代码重复

1.操作准备流程

## 70个基因预后特征(PS)， 49个基因转移负担特征(MBS)， 354个基因残留肿瘤特征(RTS)
mammaprint_long <- read.table("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

zenawerb_long <- read.table("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

artega_long <- read.table("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

2.可视化输入文件准备

#预后特征热图
prognosis_sig <- cbind(mammaprint_avg_exprs, zenawerb_avg_exprs, artega_avg_exprs)
colnames(prognosis_sig) <- c("PS", "MBS", "RTS")

#热图输入文件准备
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]])))
}

3.可视化

#热图可视化
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],
          show_heatmap_legend = FALSE,
          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,
          show_heatmap_legend = 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,
          show_heatmap_legend = 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,
          show_heatmap_legend = 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,
          show_heatmap_legend = 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 = "Expression",title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9)))
#pdf("fig4a.pdf", onefile = FALSE, width = 20)
print(draw(ht_sep_prognosis, annotation_legend_side = "right"))
#dev.off()

图片展示

Fig4b

代码重复

1.操作准备流程

# 预后特征的小提琴排序图
clust_avg_prognosis <- matrix(NA, nrow = length(unique(HSMM_allepith_clustering$Cluster)), ncol = ncol(prognosis_sig))
rownames(clust_avg_prognosis) <- paste("clust", c(1:length(unique(HSMM_allepith_clustering$Cluster))), sep = "")
colnames(clust_avg_prognosis) <- colnames(prognosis_sig)
for (c in 1:length(unique(HSMM_allepith_clustering$Cluster))) {
  clust_avg_prognosis[c,] <- apply(prognosis_sig[which(HSMM_allepith_clustering$Cluster == c),], 2, mean)}

prognosis_sig <- as.data.frame(prognosis_sig)
all.equal(rownames(prognosis_sig), colnames(HSMM_allepith_clustering))
prognosis_sig$Cluster <- as.numeric(HSMM_allepith_clustering$Cluster)

prognosis_melt <- melt(prognosis_sig, id.vars = c("Cluster"))
prognosis_melt$value <- as.numeric(prognosis_melt$value)
prognosis_melt <- prognosis_melt %>%
  filter(Cluster %in% c(2,3,4))
prognosis_melt$Cluster<-as.factor(prognosis_melt$Cluster)

2.可视化操作

#pdf("fig4b.pdf", width = 9)
p <- ggplot(prognosis_melt, aes(Cluster, value, fill = Cluster)) +
  scale_fill_manual(values = c("1" = "#ee204d", "2" = "#17806d", "3" = "#b2ec5d", "4" = "#cda4de", "5" = "#1974d2"))
p + geom_violin(adjust = .5) + facet_wrap(~variable) + stat_summary(fun.y = mean, geom = "point", shape = 22, size = 3)
#dev.off()

图片展示

Fig4d

代码重复

1.操作准备流程

# 通路热图
path1 <- c("ST3GAL4", "ST3GAL6", "ST8SIA1", "FUT1", "FUT2", "FUT3", "FUT4", "FUT6", "FUT9", "GLTP", "SPTLC1", "KDSR", "SPTLC2", "CERK", "ARSA")
idx_path1 <- match(path1, rownames(HSMM_allepith_clustering))

path2 <- c("CCL20", "CCL22", "CCL4", "CCR6", "IL11", "IL12RB1", "IL6R", "CSF2RA", "BMP7", "GLMN", "GPI", "INHBA", "TNF", 
           "TNFSF15", "GHR", "LEPR", "TLR1", "TLR2", "TLR5", "TLR7", "TLR10", "F11R")
idx_path2 <- match(path2, rownames(HSMM_allepith_clustering))

path3 <- c("ERBB2", "AKT1", "AKT3", "PIK3R3", "PIK3R4", "RPS6KB2", "TRIB3", "BTK", "GRB10", "GRB2", "ILK", "PAK1", "PRKCZ", "CSNK2A1",
           "IRS1", "IRAK1", "MYD88", "MAP2K1", "MAPK8", "MAPK1", "PTPN11", "EIF4E", "EIF4EBP1", "EIF4G1", "FKBP1A", "PDK1", "RHEB", "RPS6KB1")
idx_path3 <- match(path3, rownames(HSMM_allepith_clustering))

all_idx_paths <- c(idx_path1, idx_path2)

names_path <- c(rep("glycosphigolipid metabolism", length(idx_path1)), rep("innate immunity", length(idx_path2)))
anno_cols_path <- c("glycosphigolipid metabolism" = "#ff9baa", "innate immunity" = "#17806d")
ha_path_rows <- HeatmapAnnotation(df = data.frame(pathway = names_path),
                                  annotation_legend_param = list(pathway = list(ncol = 1, title = "pathway", title_position = "topcenter")),
                                  which = "row", col = list("pathway" = anno_cols_path), annotation_width = unit(3, "mm"))
# 每个簇分开矩阵，只提取相关基因
mat_epith_allpats <- exprs(HSMM_allepith_clustering)
mats_epith_patient <- list()
pds_epith_patient <- list()
mats_epith_patient_clusters <- list()
for (i in 1:length(patients_now)) {
  pds_epith_patient[[i]] <- pData(HSMM_allepith_clustering)[which(pData(HSMM_allepith_clustering)$patient == patients_now[i]), ]
  mats_epith_patient[[i]] <- mat_epith_allpats[all_idx_paths, which(pData(HSMM_allepith_clustering)$patient == patients_now[i])]
  mats_epith_patient_clusters[[i]] <- list()
  for (c in 1:length(unique(HSMM_allepith_clustering$Cluster))) {
    mats_epith_patient_clusters[[i]][[c]] <- mats_epith_patient[[i]][, which(pds_epith_patient[[i]]$Cluster == c)]
  }
  names(mats_epith_patient_clusters[[i]]) <- paste0("clust", c(1:5))
}
names(mats_epith_patient) <- patients_now
names(pds_epith_patient) <- patients_now

2.可视化

# 热图绘制
ht_paths <-
  ha_path_rows + 
  Heatmap(mats_epith_patient_clusters[[1]][[2]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE,
          name = "cluster2", column_title = "2",
          show_row_names = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE,
          row_names_gp = gpar(fontsize = 9, col = anno_cols_path),
          split = names_path,
          heatmap_legend_param = list(title = "Expression",title_gp = gpar(fontsize = 9), labels_gp = gpar(fontsize = 9)))+
  Heatmap(mats_epith_patient_clusters[[1]][[3]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE, 
          name = "cluster3", column_title = "3",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
  Heatmap(mats_epith_patient_clusters[[1]][[4]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE, 
          name = "cluster4", column_title = "4",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
  Heatmap(mats_epith_patient_clusters[[2]][[2]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE,
          name = "cluster2_2", column_title = "2",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
  Heatmap(mats_epith_patient_clusters[[2]][[3]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE, 
          name = "cluster3_2", column_title = "3",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
  Heatmap(mats_epith_patient_clusters[[2]][[4]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE, 
          name = "cluster4_2", column_title = "4",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
 Heatmap(mats_epith_patient_clusters[[3]][[2]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE,
          name = "cluster2_3", column_title = "2",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
 Heatmap(mats_epith_patient_clusters[[3]][[4]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE, 
          name = "cluster4_3", column_title = "4",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
Heatmap(mats_epith_patient_clusters[[4]][[2]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE,
          name = "cluster2_4", column_title = "2",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
  Heatmap(mats_epith_patient_clusters[[4]][[3]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE, 
          name = "cluster3_4", column_title = "3",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
  Heatmap(mats_epith_patient_clusters[[4]][[4]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE, 
          name = "cluster4_4", column_title = "4",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
  Heatmap(mats_epith_patient_clusters[[4]][[5]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE, 
          name = "cluster5_4", column_title = "5",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
 Heatmap(mats_epith_patient_clusters[[5]][[1]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE,
          name = "cluster1_5", column_title = "1",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
  Heatmap(mats_epith_patient_clusters[[5]][[2]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE,
          name = "cluster2_5", column_title = "2",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
  Heatmap(mats_epith_patient_clusters[[5]][[3]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE, 
          name = "cluster3_5", column_title = "3",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
  Heatmap(mats_epith_patient_clusters[[5]][[5]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE, 
          name = "cluster5_5", column_title = "5",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
  Heatmap(mats_epith_patient_clusters[[6]][[2]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE,
          name = "cluster2_6", column_title = "2",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
  Heatmap(mats_epith_patient_clusters[[6]][[3]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE, 
          name = "cluster3_6", column_title = "3",
          show_row_names = FALSE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE) +
  Heatmap(mats_epith_patient_clusters[[6]][[4]],
          col = colorRamp2(c(-0.5, 1, 3), c("blue", "white", "red")),
          cluster_rows = TRUE, show_row_dend = FALSE, 
          name = "cluster4_6", column_title = "4",
          show_row_names = TRUE,show_heatmap_legend = FALSE,
          show_column_dend = TRUE, show_column_names = FALSE)
#pdf("fig4d.pdf", onefile = FALSE, width = 25)
draw(ht_paths)
#dev.off()

图片展示

(完)

全部的代码和数据，都可以在我们《生信技能树》公众号后台回复“tnbc”获取，还是很有学习的价值。其实我更希望看到学徒跟着我的《CNS图表复现专栏》来整理这些代码；