在上一節(jié)，由于大部分細(xì)胞（868個(gè)）都被歸為上皮細(xì)胞群中（Fig2 c），這868個(gè)細(xì)胞可被分成5個(gè)cluster，接著對(duì)這5個(gè)cluster細(xì)胞進(jìn)行探索。我們使用一組來自對(duì)乳腺腫塊的非監(jiān)督分析的基因表達(dá)特征對(duì)5個(gè)cluster進(jìn)行了研究。這些基因表達(dá)特征通過比較三陰性乳腺癌（TNBC）的四個(gè)亞型(ERBB2 amplicon，Luminal Subtype 、Basal epithelial-cell enriched 和Luminal epithelial gene cluster containing ER)而建立。先看看這5個(gè)clusters的basal細(xì)胞來源的細(xì)胞群有多少。大多數(shù)TNBC是基底樣腫瘤，它們與多種TNBC型亞型重疊，與非固有基底TNBCs相比，與克隆異質(zhì)性增加有關(guān)。（備注：這篇文獻(xiàn)用到了很多apply循環(huán)，大家仔細(xì)琢磨，大概意思能看懂就行，然后可以把它應(yīng)用到自己的數(shù)據(jù)中）
## 讀取數(shù)據(jù) "data/genes_for_basal_vs_non_basal_tnbc_PNAS.txt" , header = TRUE, sep = "\t" )#提取Basal.epithelial.cell.enriched.cluster的基因 $Basal .epithelial.cell.enriched.cluster#合并剩下17個(gè)基因 "SOX9" "GALNT3" "CDH3" "LAMC2" "CX3CL1" "TRIM29" "KRT17" "KRT5" "CHI3L2" "SLPI" "NFIB" "MRAS" "TGFB2" "CAPN6" "DMD" "FABP7" "CXCL1" #算出17個(gè)basal_PNAS基因在1112個(gè)細(xì)胞的表達(dá)平均值 #檢查一下數(shù)據(jù) #提取868個(gè)上皮細(xì)胞群體的17個(gè)basal_PNAS基因表達(dá)平均值 which (pd_ct$cell_types_cl_all == "epithelial" )]#檢查一下數(shù)據(jù) #把17個(gè)basal_PNAS基因表達(dá)平均值賦給HSMM_allepith_clustering，以便于后續(xù)分析 $basal_PNAS_avg_exprs <- basal_PNAS_avg_exprs#畫figS9b "basal_PNAS_avg_exprs" , cell_size = 2) + facet_wrap(~patient) +"yellow" , high = "blue" )#畫figS9a "basal_PNAS_avg_exprs" , cell_size = 2) + facet_wrap(~Cluster) +"yellow" , high = "blue" )figS9a:大多數(shù)TNBC樣本都有basal gene signature的表達(dá)。figS9b:在868個(gè)上皮細(xì)胞群中，cluster2的basal gene signature表達(dá)量最豐富。
接著，使用另外一個(gè)基因表達(dá)特征數(shù)據(jù)集TNBCtype4 signatures（Lehman_signature），這個(gè)signatures根據(jù)基因表達(dá)變化將TNBC細(xì)胞分為6個(gè)類：basal_like_1、basal_like_2、immunomodulatory、mesenchymal、mesenchymal_stem_like和luminal_ar。作者將基因表達(dá)特征中上調(diào)基因的平均表達(dá)值減去下調(diào)基因的平均表達(dá)值，將差值作為每個(gè)細(xì)胞在TNBCtype4 signatures （basal_like_1、basal_like_2、mesenchymal和luminal_ar）中的每個(gè)基因表達(dá)值，挑選最高基因表達(dá)值對(duì)應(yīng)的signature，將其分配給對(duì)應(yīng)的細(xì)胞。
#讀取數(shù)據(jù) "data/Lehman_signature.txt" , sep = "\t" , header = TRUE, stringsAsFactors = FALSE)#這個(gè)for循環(huán)提取了lehman_long里面的四列g(shù)ene、regulation、no_samples和signature，建立一個(gè)data.rrame for (i in 0:5) {"gene" "regulation" "no_samples" "signature" if (i == 0) {which (lehman$signature == "" ),]if (i > 0) {"gene" , i, sep = "." )"regulation" , i, sep = "." )"no_samples" , i, sep = "." )"signature" , i, sep = "." )"gene" , "regulation" , "no_samples" , "signature" )if (length(which (is.na(mat_to_bind$no_samples ))) > 0 )which (mat_to_bind$signature == "" ),]#刪掉一些mat_ct沒有檢測(cè)到的基因 which (!is.na(match(lehman$gene , rownames(mat_ct)))),]$signature )function (x){for (s in 1:length(lehman_signatures)) {# current signature "UP" )"DOWN" )$gene , rownames(mat_ct)) $gene , rownames(mat_ct))#算上調(diào)、下調(diào)的基因在樣本中的平均表達(dá)值 return (mns)#檢查數(shù)據(jù) #只看868個(gè)上皮細(xì)胞的情況 which (pd_ct$cell_types_cl_all == "epithelial" )]#提取lehman_avg_exps前面6行，對(duì)應(yīng)的是up which (pd_ct$cell_types_cl_all == "epithelial" )]#提取lehman_avg_exps后面6行，對(duì)應(yīng)的是down which (pd_ct$cell_types_cl_all == "epithelial" )]#上調(diào)基因的平均表達(dá)值減去下調(diào)基因的平均表達(dá)值 #挑選最高基因表達(dá)值對(duì)應(yīng)的signature，將其分配給對(duì)應(yīng)的細(xì)胞。 function (x){rownames(lehman_avg_both)[which.max(x)]})which (pd_ct$cell_types_cl_all == "epithelial" )]#刪除immunomodulatory和mesenchymal_stem_like signature which (rownames(lehman_avg_both_epithelial) %in % c("immunomodulatory" , "mesenchymal_stem_like" )),]function (x){rownames(lehman_avg_both_epithelial_new)[which.max(x)]})接下來畫圖，同樣地，需要對(duì)heatmap函數(shù)代碼進(jìn)行修改。
ha_lehman_epith_pat <- list()for (i in 1:length(patients_now)) {if (i == 1)"1" = "#ee204d" , "2" = "#17806d" , "3" = "#b2ec5d" , "4" = "#cda4de" , "5" = "#1974d2" )),"left" , annotation_name_gp = gpar(fontsize = 12),"topcenter" , title = "cluster" )),"mm" ),if (i > 1 && i != 5 )"1" = "#ee204d" , "2" = "#17806d" , "3" = "#b2ec5d" , "4" = "#cda4de" , "5" = "#1974d2" )),"left" , annotation_name_gp = gpar(fontsize = 12),"topcenter" , title = "cluster" )),"mm" ),if (i == 5)"1" = "#ee204d" , "2" = "#17806d" , "3" = "#b2ec5d" , "4" = "#cda4de" , "5" = "#1974d2" )),"right" , annotation_name_gp = gpar(fontsize = 12),"topcenter" ,title = "cluster" )),"mm" ),#檢查數(shù)據(jù) #將basal signature添加進(jìn)去，以便后續(xù)作圖 for (i in 1:length(patients_now)) {$basal_PNAS_avg_exprs [which (HSMM_allepith_clustering$patient == patients_now[i])])"intrinsic_basal" # 畫圖 "blue" ,"white" , "red" )),"blue" ,"white" , "red" )),"blue" ,"white" , "red" )),"blue" ,"white" , "red" )),"blue" ,"white" , "red" )),"blue" ,"white" , "red" )),"right" ,#畫fig3d print (draw(ht_sep_lehmans_both_wbasal_new, annotation_legend_side = "right" ))我們只需要把右邊注釋條PS一下，就可以達(dá)到跟文獻(xiàn)的圖片一模一樣了。
#檢查數(shù)據(jù) $assignments_lehman_both_new <- assignments_lehman_both_epithelial_new"assignments_lehman_both_new" , cell_size = 2) + facet_wrap(~patient)#畫figS8 "assignments_lehman_both_new" , cell_size = 2) + facet_wrap(~Cluster)cluster4也富集了 Basal-Like 1 signature，而cluster3高度富集了TNBCtype 中的“Luminal Androgen Receptor” signature。為了更清楚的看到上皮細(xì)胞群的5個(gè)cluster對(duì)應(yīng)的TNBCtype signatures的平均表達(dá)量，接著繼續(xù)探索下去...
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 "clust" , c(1:length(unique(HSMM_allepith_clustering$Cluster ))), sep = "" )#行名：basal_like_1、basal_like_2、mesenchymal、luminal_ar" #算出每個(gè)cluster的signatures 平均值 for (c in 1:length(unique(HSMM_allepith_clustering$Cluster ))) {which (HSMM_allepith_clustering$Cluster == c)], 1, mean)#增加一列cluster $Cluster <- rownames(clust_avg_lehman_both_new)#拆分?jǐn)?shù)據(jù) "Cluster" )#畫fig3e #guides(colour = guide_legend(override.aes = list(size=3))) + "average expression of signature in cluster" ) +"cluster" ) +可以看到：cluster2和4強(qiáng)烈地表達(dá)Basal-Like 1 signature，而cluster3顯著表達(dá)Basal-Like 2 signature和 luminal_AR signature
接著，讀取另外一個(gè)ML_signature(mature luminal signature)，將上調(diào)基因的平均表達(dá)量中減去下調(diào)基因的平均表達(dá)量，計(jì)算出three normal breast signatures下每個(gè)細(xì)胞的表達(dá)量(Lim et al., 2009a)，然后將每個(gè)細(xì)胞分配給其具有最高表達(dá)量的signatures。這三個(gè)normal breast signatures 是：mature luminal (ML)，basal和luminal progenitor (LP)，在每個(gè)signatures中，都有對(duì)應(yīng)的上調(diào)基因和下調(diào)基因。
"data/ML_signature.txt" , sep = "\t" , header = TRUE)if (length(which (ml_signature_long$Symbol == "" )) > 0)#將空白行去掉 which (ml_signature_long$Symbol == "" ),]#按照基因字母進(jìn)行排序，如果基因字母有一樣的，則按照Average.log.fold.change絕對(duì)值的負(fù)數(shù)進(jìn)行從小到大排序 $Symbol , -abs(ml_signature_long$Average .log.fold.change) ), ]#對(duì)基因取唯一值，去重復(fù) $Symbol ), ]#總共有818個(gè)基因 which (!is.na(match(ml_signature_long$Symbol , rownames(mat_ct)))), ]#上調(diào)基因有384個(gè) which (ml_signature$Average .log.fold.change > 0), ]#下調(diào)基因有197個(gè) which (ml_signature$Average .log.fold.change < 0), ]#匹配一下 $Symbol , rownames(mat_ct))$Symbol , rownames(mat_ct))#讀取basal signature，處理過程跟上面的一樣的。 "data/basal_signature.txt" , sep = "\t" , header = TRUE)if (length(which (basal_signature_long$Symbol == "" )) > 0)which (basal_signature_long$Symbol == "" ),]$Symbol , -abs(basal_signature_long$Average .log.fold.change) ), ]$Symbol ), ]#總共有1335個(gè)基因 which (!is.na(match(basal_signature_long$Symbol , rownames(mat_ct)))), ]#上調(diào)基因有588個(gè) which (basal_signature$Average .log.fold.change > 0), ]#下調(diào)基因有757個(gè) which (basal_signature$Average .log.fold.change < 0), ]$Symbol , rownames(mat_ct))$Symbol , rownames(mat_ct))#讀取LP signature,還是同樣的操作 "data/Lp_signature.txt" , sep = "\t" , header = TRUE)if (length(which (lp_signature_long$Symbol == "" )) > 0)which (lp_signature_long$Symbol == "" ),]$Symbol , -abs(lp_signature_long$Average .log.fold.change) ), ]$Symbol ), ]which (!is.na(match(lp_signature_long$Symbol , rownames(mat_ct)))), ]which (lp_signature$Average .log.fold.change > 0), ]which (lp_signature$Average .log.fold.change < 0), ]$Symbol , rownames(mat_ct))$Symbol , rownames(mat_ct))#對(duì)ML、basal和LP 3個(gè)signatures基因，將上調(diào)基因的表達(dá)值減去下調(diào)基因表達(dá)值，并分別返回結(jié)果。 function (x){return (c(avg_ml_up, avg_basal_up, avg_lp_up, avg_ml_both, avg_basal_both, avg_lp_both))"avg_ml_up" , "avg_basal_up" , "avg_lp_up" , "avg_ml_both" , "avg_basal_both" , "avg_lp_both" )#檢查數(shù)據(jù) #只看上皮細(xì)胞群 which (pd_ct$cell_types_cl_all == "epithelial" )]which (pd_ct$cell_types_cl_all == "epithelial" )]which (pd_ct$cell_types_cl_all == "epithelial" )]#挑選最大值=上調(diào)基因的平均表達(dá)值最大數(shù)值分配給對(duì)應(yīng)的細(xì)胞類型 function (x){rownames(normsig_avg_ups)[which.max(x)]})which (pd_ct$cell_types_cl_all == "epithelial" )]#上調(diào)基因的平均表達(dá)值-下調(diào)基因的平均表達(dá)值的最大數(shù)值分配給對(duì)應(yīng)的細(xì)胞類型 function (x){rownames(normsig_avg_both)[which.max(x)]})which (pd_ct$cell_types_cl_all == "epithelial" )]# heatmaps on normal signatures per patient which (pd_ct$cell_types_cl_all == "epithelial" ),]for (i in 1:length(patients_now)) {which (pd_ct_epith$patient == patients_now[i])]which (pd_ct_epith$patient == patients_now[i])]which (pds_ct[[i]]$cell_types_cl_all == "epithelial" ),]"blue" ,"white" , "red" )),"blue" ,"white" , "red" )),"blue" ,"white" , "red" )),"cm" ),"blue" ,"white" , "red" )),"cm" ),"blue" ,"white" , "red" )),"cm" ),"blue" ,"white" , "red" )),"right" ,"cm" ),#畫fig3b.pdf print (draw(ht_sep_normsig_both, annotation_legend_side = "bottom" ))# 每個(gè)樣本的normal signatures 數(shù)目 $assignments_normsig_both <- assignments_normsig_both_epithelial$assignments_normsig_ups <- assignments_normsig_ups_epithelial#畫fig3f "assignments_normsig_both" , cell_size = 2) + facet_wrap(~patient)#每個(gè)clusters 的normal signatures 數(shù)目 "assignments_normsig_both" , cell_size = 2) + facet_wrap(~Cluster)#檢查數(shù)據(jù) $Cluster , clustering_allepith)$Cluster )), ncol = nrow(normsig_avg_both_epithelial))"clust" , c(1:length(unique(HSMM_allepith_clustering$Cluster ))), sep = "" )#算上皮細(xì)胞群的avg_both 平均表達(dá)值，接下來還是同樣的操作 for (c in 1:length(unique(HSMM_allepith_clustering$Cluster ))) {which (HSMM_allepith_clustering$Cluster == c)], 1, mean)$Cluster <- rownames(clust_avg_normsig_both)"Cluster" )#畫fig3e "average expression of signature in cluster" ) +"cluster" ) +fig3e：Clusters 2和Clusters4 強(qiáng)烈表達(dá) LP signature, 而cluster 3則高表達(dá) ML signature.
接著為了進(jìn)一步探究臨床相關(guān)性，作者使用了三個(gè)臨床相關(guān)的gene signatures，進(jìn)一步探究這868個(gè)上皮細(xì)胞的特征，這868個(gè)上皮細(xì)胞真的被研究到很徹底，真的佩服，這工作量好大....
第一個(gè)gene signatures:70-gene prognostic signature ，該signatures最初是從對(duì)有無轉(zhuǎn)移復(fù)發(fā)患者的原發(fā)腫瘤之間差異表達(dá)基因的分析中得出的，總共70個(gè)基因。
"data/mammaprint_sig_new.txt" , header = TRUE, sep = "\t" )function (x){return (intersect(x, rownames(mat_ct)))})[,1]which (pd_ct$cell_types_cl_all == "epithelial" )]$mammaprint_avg_exprs <- mammaprint_avg_exprs# 畫figS13b "mammaprint_avg_exprs" , cell_size = 2) + facet_wrap(~patient) +"yellow" , high = "blue" )# 畫figS13a "mammaprint_avg_exprs" , cell_size = 2) + facet_wrap(~Cluster) +"yellow" , high = "blue" )第二個(gè)gene signatures:49-gene metastatic burden signature.該signatures可以區(qū)分了患者來源的小鼠TNBC異種移植模型中單個(gè)循環(huán)轉(zhuǎn)移細(xì)胞所產(chǎn)生的高轉(zhuǎn)移負(fù)荷和低轉(zhuǎn)移負(fù)荷，總共包括49個(gè)基因。
zenawerb_long <- read.table("data/werb_49_metastasis_sig.txt" , header = TRUE, sep = "\t" )function (x){return (intersect(x, rownames(mat_ct)))})[,1]which (pd_ct$cell_types_cl_all == "epithelial" )]$zenawerb_avg_exprs <- zenawerb_avg_exprs#畫figS14b "zenawerb_avg_exprs" , cell_size = 2) + facet_wrap(~patient) +"yellow" , high = "blue" )#畫figS14a "zenawerb_avg_exprs" , cell_size = 2) + facet_wrap(~Cluster) +"yellow" , high = "blue" )第三個(gè)gene siganatures:從接受手術(shù)前化療治療的原發(fā)性乳腺癌患者的殘存腫瘤群體中富集的基因中獲得的，包括354個(gè)基因。
artega_long <- read.table("data/artega_sig.txt" , header = TRUE, sep = "\t" )function (x){return (intersect(x, rownames(mat_ct)))})[,1]which (pd_ct$cell_types_cl_all == "epithelial" )]$artega_avg_exprs <- artega_avg_exprs"artega_avg_exprs" , cell_size = 2) + facet_wrap(~patient) +"yellow" , high = "blue" )#畫figS15b "artega_avg_exprs" , cell_size = 2) + facet_wrap(~Cluster) +"yellow" , high = "blue" )#將3個(gè)gene signatures的表達(dá)值合成一個(gè)數(shù)據(jù)框 #取行名 "mammaprint" , "zenawerb" , "artega" )for (i in 1:length(patients_now)) {which (pd_ct_epith$patient == patients_now[i])]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]])))"blue" ,"white" , "red" )),"blue" ,"white" , "red" )),"blue" ,"white" , "red" )),"blue" ,"white" , "red" )),"blue" ,"white" , "red" )),"blue" ,"white" , "red" )),"right" ,#畫fig4a print (draw(ht_sep_prognosis, annotation_legend_side = "right" ))
