所有的肿瘤都有恶性增殖的特性吗

生信技能树

发布于 2022-07-26 10:21:53

2660

发布于 2022-07-26 10:21:53

文章被收录于专栏：生信技能树

前面我们在：癌基因一定在肿瘤部位高表达吗 我们针对每个癌症都在各种内部做了肿瘤组织和正常对照的差异表达量分析，然后在癌基因都是肿瘤的风险因子吗 我们针对每个癌症的全部基因批量了做了单基因的cox分析。

然后发现很多癌症都有MKI67和TOP2A这样的基因的高表达，而且它的高表达是坏的预后。我们就有一个自然而然的假设，就是：是否所有的肿瘤都有恶性增殖的特性呢？

批量针对差异分析后的基因排序进行GSEA分析

在前面的教程：居然有如此多种癌症（是时候开启pan-cancer数据挖掘模式），我们把全部的TCGA的33种癌症的表达量矩阵区拆分成为蛋白编码基因和非编码基因这两个不同的表达量矩阵，并且保存成为了rdata文件。然后在：癌基因一定在肿瘤部位高表达吗 我们仅仅是针对normal样品数量大于30的癌症做了差异分析，12个癌症的样品数量满足要求，每个差异分析都是得到如下所示的矩阵：

> 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这一列进行排序，现在我们使用最常用的KEGG数据库的约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函数里面的 pvalueCutoff = 0.9, 所以并不是每次差异分析的gsea分析都会返回全部的KEGG数据库的约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

需要把结果对齐后，才能成为一个行是通路，列是癌症的，数值是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

如下所示：

行是通路，列是癌症的，数值是gsea的es打分的矩阵

虽然这个时候的矩阵gsea_df里面确实是有全部的KEGG数据库的 333 条通路，但是可以看到因为我的gseKEGG函数里面的 pvalueCutoff = 0.9, 所以其实去除NA后，就100多通路是所有的癌症返回了的。

不过不重要，我们这个时候就是看看细胞增殖相关的通路。在 https://www.genome.jp/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里面是完全反过来的规律，蛮奇怪的。

代码也很简单的：

pheatmap::pheatmap(gsea_df[cg,],display_numbers = T,number_format =  "%.2f")

也仅仅是，基本上可以确定所有的肿瘤都有恶性增殖的特性，那么同理，我们可以查看肿瘤是否都有其它特性：

2000年，麻省理工学院的Robert A. Weinberg教授和Douglas Hanahan教授共同在Cell 杂志上发表了第一版“The Hallmarks of Cancer”。
- 1、自给自足的生长信号（self-sufficiency in growth signals）：如H-RAS致癌基因的激活。
- 2、对生长抑制信号不敏感（insensitivity to anti-growth signals）：如视网膜母细胞瘤（Rb）基因（肿瘤抑制基因）的缺失。
- 3、逃避凋亡（evading apoptosis）：如生成IGF因子。
- 4、无限复制的潜力（limitless replicative potential）：如激活端粒酶。
- 5、持续的血管生成（sustained angiogenesis）：如生成VEGF诱导物。
- 6、组织侵袭和转移（tissue invasion & metastasis）：如E-钙黏蛋白的失活。
2011年，两位教授根据十年间对肿瘤研究的最新进展和研究范式的更迭，与时俱进地对原来的6大癌症标志性特征进行了更新和拓展，撰写的“Hallmarks of cancer: the next generation”仍旧发表在Cell 杂志上，增加了4个新的标志性特征，分别是：
- 细胞能量代谢的失控（deregulating cellular energetics）
- 逃避免疫清除（avoiding immune destruction）
- 肿瘤促炎症作用（tumor-promoting inflammation）
- 基因组的不稳定性和突变（genome instability and mutation）
2022年的第三版“Hallmarks of Cancer”，发表在国际顶级刊物Cancer Discovery，总结的14大癌症标志性特征，就是说增加了4个新的标志性特征，分别是：
- 解锁表型可塑性（Unlocking phenotypic plasticity）
- 非突变表观遗传重编程（Non-mutational epigenetic reprogramming）
- 多态性的微生物组（Polymorphic microbiomes）
- 衰老细胞（Senescent cells）。