前面我們在：癌基因一定在腫瘤部位高表達嗎 我們針對每個癌癥都在各種內(nèi)部做了腫瘤組織和正常對照的差異表達量分析，然后在癌基因都是腫瘤的風險因子嗎 我們針對每個癌癥的全部基因批量了做了單基因的cox分析。
然后發(fā)現(xiàn)很多癌癥都有MKI67和TOP2A這樣的基因的高表達，而且它的高表達是壞的預后。我們就有一個自然而然的假設，就是：是否所有的腫瘤都有惡性增殖的特性呢？
批量針對差異分析后的基因排序進行GSEA分析
在前面的教程：居然有如此多種癌癥（是時候開啟pan-cancer數(shù)據(jù)挖掘模式），我們把全部的TCGA的33種癌癥的表達量矩陣區(qū)拆分成為蛋白編碼基因和非編碼基因這兩個不同的表達量矩陣，并且保存成為了rdata文件。然后在：癌基因一定在腫瘤部位高表達嗎 我們僅僅是針對normal樣品數(shù)量大于30的癌癥做了差異分析，12個癌癥的樣品數(shù)量滿足要求，每個差異分析都是得到如下所示的矩陣：
> head(deg_list[[1]])
logFC AveExpr t P.Value adj.P.Val B
FIGF -5.987848 -0.6548567 -58.33261 0.000000e+00 0.000000e+00 803.4724
PAMR1 -3.961586 2.4372420 -49.05257 7.740891e-291 7.070143e-287 655.9085
CD300LG -6.577422 -0.0413564 -48.39557 4.067331e-286 2.476598e-282 644.9867
LYVE1 -4.777174 1.4285749 -47.82048 5.728237e-282 2.615943e-278 635.4879
CA4 -6.961439 -2.5818286 -46.47851 3.145555e-272 1.149197e-268 612.8327
SCARA5 -6.458426 0.3467161 -45.88121 7.217392e-268 2.197335e-264 603.0621
這12個差異分析都是針對原始的counts矩陣做了最簡單的limma的voom算法，得到的兩萬個基因都是可以通過logFC這一列進行排序，現(xiàn)在我們 使用最常用的KEGG數(shù)據(jù)庫的約300條通路對12個癌癥的差異分析結果進行批量GSEA分析。
load(file = 'batch_deg_list.Rdata')
gsea_list = lapply(deg_list, function(DEG_limma_voom){
nrDEG=DEG_limma_voom[,c(1,4)]
colnames(nrDEG)=c('log2FoldChange','pvalue')
head(nrDEG)
library(org.Hs.eg.db)
library(clusterProfiler)
gene <- bitr(rownames(nrDEG), fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)
gene$logfc <- nrDEG$log2FoldChange[match(gene$SYMBOL,rownames(nrDEG))]
head(gene)
geneList=gene$logfc
names(geneList)=gene$ENTREZID
geneList=sort(geneList,decreasing = T)
head(geneList)
library(clusterProfiler)
kk_gse <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 1000,
minGSSize = 10,
pvalueCutoff = 0.9,
verbose = FALSE)
tmp=kk_gse@result
kk=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
return(kk@result)
})
因為要使用clusterProfiler包，所以不得不進行繁瑣的基因ID轉換，然后因為我的gseKEGG函數(shù)里面的 pvalueCutoff = 0.9, 所以并不是每次差異分析的gsea分析都會返回全部的KEGG數(shù)據(jù)庫的約300條通路。如下所示：
> head(gsea_list[[1]][,1:6])
ID Description setSize enrichmentScore NES pvalue
hsa04151 hsa04151 PI3K-Akt signaling pathway 332 -0.3782177 -1.540106 0.001242236
hsa04080 hsa04080 Neuroactive ligand-receptor interaction 326 -0.3951545 -1.605126 0.001245330
hsa04020 hsa04020 Calcium signaling pathway 233 -0.4464945 -1.752872 0.001310616
hsa04014 hsa04014 Ras signaling pathway 224 -0.4180371 -1.633703 0.001326260
hsa04024 hsa04024 cAMP signaling pathway 212 -0.4554115 -1.777863 0.001328021
hsa04015 hsa04015 Rap1 signaling pathway 202 -0.4137399 -1.603427 0.001340483
> do.call(rbind,lapply(gsea_list, dim))
[,1] [,2]
BRCA 291 11
COAD 282 11
HNSC 280 11
KIRC 262 11
KIRP 306 11
LIHC 271 11
LUAD 296 11
LUSC 268 11
PRAD 270 11
STAD 272 11
THCA 247 11
UCEC 265 11
需要把結果對齊后，才能成為一個行是通路，列是癌癥的，數(shù)值是gsea的es打分的矩陣。
allkid= unique(unlist(lapply(gsea_list,'[',2)))
gsea_df = do.call(cbind,lapply(gsea_list,function(x){
x[match(allkid,x[,2]),4]
}))
rownames(gsea_df)=allkid
如下所示：
行是通路，列是癌癥的，數(shù)值是gsea的es打分的矩陣雖然這個時候的矩陣gsea_df里面確實是有全部的KEGG數(shù)據(jù)庫的 333 條通路，但是可以看到因為我的gseKEGG函數(shù)里面的 pvalueCutoff = 0.9, 所以其實去除NA后，就100多通路是所有的癌癥返回了的。
不過不重要，我們這個時候就是看看細胞增殖相關的通路。在 https://www./kegg/pathway.html 可以找到一下，包括：
2.4 Replication and repair
03030 DNA replication
03410 Base excision repair
03420 Nucleotide excision repair
03430 Mismatch repair
03440 Homologous recombination
03450 Non-homologous end-joining
03460 Fanconi anemia pathway
4.2 Cell growth and death
04110 Cell cycle
那，我們就單獨提取它們并且可視化：
cg = c(
'DNA replication','Base excision repair',
'Nucleotide excision repair','Mismatch repair',
'Homologous recombination','Non-homologous end-joining',
'Fanconi anemia pathway','Cell cycle'
)
gsea_df[cg,1:4]
可以看到：
> gsea_df[cg,1:4]
BRCA COAD HNSC KIRC
DNA replication 0.6542816 0.6613777 0.6909022 0.5131581
Base excision repair 0.5654879 0.5235520 0.4837935 0.4240750
Nucleotide excision repair 0.4617563 0.5185588 0.4338332 0.3560308
Mismatch repair 0.5843370 0.6097009 0.5569984 0.4988459
Homologous recombination 0.5899100 0.6590182 0.5875136 0.5413358
Non-homologous end-joining 0.4208269 0.4679819 0.4294577 NA
Fanconi anemia pathway 0.5920658 0.6461166 0.5570080 0.5225555
Cell cycle 0.5683101 0.5176032 0.5948943 0.4815887
都是顯而易見的正向GSEA富集，因為它們的ES值都大于 0.5 , 當然了，可視化可能是效果更好一點：
可視化可能是效果更好一點不知道為什么Non-homologous end-joining 這個 通路在THCA里面是完全反過來的規(guī)律，蠻奇怪的。
代碼也很簡單的：
pheatmap::pheatmap(gsea_df[cg,],display_numbers = T,number_format = "%.2f")
也僅僅是，基本上可以確定所有的腫瘤都有惡性增殖的特性，那么同理，我們可以查看腫瘤是否都有其它特性：
- 2000年，麻省理工學院的Robert A. Weinberg教授和Douglas Hanahan教授共同在Cell 雜志上發(fā)表了第一版“The Hallmarks of Cancer”。
- 1、自給自足的生長信號（self-sufficiency in growth signals）：如H-RAS致癌基因的激活。
- 2、對生長抑制信號不敏感（insensitivity to anti-growth signals）：如視網(wǎng)膜母細胞瘤（Rb）基因（腫瘤抑制基因）的缺失。
- 3、逃避凋亡（evading apoptosis）：如生成IGF因子。
- 4、無限復制的潛力（limitless replicative potential）：如激活端粒酶。
- 5、持續(xù)的血管生成（sustained angiogenesis）：如生成VEGF誘導物。
- 6、組織侵襲和轉移（tissue invasion & metastasis）：如E-鈣黏蛋白的失活。
- 2011年，兩位教授根據(jù)十年間對腫瘤研究的最新進展和研究范式的更迭，與時俱進地對原來的6大癌癥標志性特征進行了更新和拓展，撰寫的“Hallmarks of cancer: the next generation”仍舊發(fā)表在Cell 雜志上，增加了4個新的標志性特征，分別是：
- 細胞能量代謝的失控（deregulating cellular energetics）
- 逃避免疫清除（avoiding immune destruction）
- 腫瘤促炎癥作用（tumor-promoting inflammation）
- 基因組的不穩(wěn)定性和突變（genome instability and mutation）
- 2022年的第三版“Hallmarks of Cancer”，發(fā)表在國際頂級刊物Cancer Discovery，總結的14大癌癥標志性特征，就是說增加了4個新的標志性特征，分別是：
- 解鎖表型可塑性（Unlocking phenotypic plasticity）
- 非突變表觀遺傳重編程（Non-mutational epigenetic reprogramming）
- 多態(tài)性的微生物組（Polymorphic microbiomes）
學徒作業(yè)
查一下14大癌癥標志性特征對應的基因列表，然后在泛癌里面看看是不是它們確實在腫瘤組織里面的相對于正常來說是激活的，也可以做一下分類模型，看看哪個癌癥標志性特征預測癌癥和正常是效果最好的。
