【原】癌基因一定在肿瘤部位高表达吗

健明 2022-06-19 发布于广东

展开全文

癌基因是基因的一类，激活后可促使正常细胞癌变、侵袭及转移。通常我们以为的激活是它表达量上升，尤其是肿瘤部位和正常组织对比的时候，而癌基因激活的方式包括点突变、基因扩增、染色体重排、病毒感染等。

而抑癌基因（tumor suppressor genes），也称肿瘤抑制基因，或俗称抗癌基因，是一类存在于正常细胞内可抑制细胞生长并具有潜在抑癌作用的基因。抑癌基因在控制细胞生长、增殖及分化过程中起着十分重要的负调节作用，它与原癌基因相互制约，维持正负调节信号的相对稳定。当这类基因在发生突变、缺失或失活时可引起细胞恶性转化而导致肿瘤的发生。所以通常我们以为他的失活是在肿瘤部位表达量下降。

那么，很自然的一个需要探索的是癌基因一定在肿瘤部位高表达吗？以及抑癌基因一定在肿瘤部位低表达吗？

下面，我们以TCGA数据库的多种癌症来举例说明这两个问题：

批量差异分析

在前面的教程：居然有如此多种癌症（是时候开启pan-cancer数据挖掘模式），我们把全部的TCGA的33种癌症的表达量矩阵区拆分成为蛋白编码基因和非编码基因这两个不同的表达量矩阵，并且保存成为了rdata文件。

在不同癌症里面，蛋白质编码相关基因数量一直在一万八附近，而非编码基因数量跨度比较大，从一万二到两万七不等。因为非编码基因的组织病人特异性都很强，但是绝大部分的癌基因和抑癌基因都是蛋白质编码基因，所以我们仅仅是针对于每个癌症的不到两万个的蛋白质编码基因在肿瘤部位和正常组织进行转录组测序表达量矩阵的差异分析。

 ls -lh |cut -d" " -f 7-|head 

  4.6M 10  1  2021 TCGA-ACC.htseq_counts.Rdata
   26M 10  1  2021 TCGA-BLCA.htseq_counts.Rdata
   78M 10  1  2021 TCGA-BRCA.htseq_counts.Rdata
   19M 10  1  2021 TCGA-CESC.htseq_counts.Rdata
  2.7M 10  1  2021 TCGA-CHOL.htseq_counts.Rdata
   29M 10  1  2021 TCGA-COAD.htseq_counts.Rdata
  2.9M 10  1  2021 TCGA-DLBC.htseq_counts.Rdata
   13M 10  1  2021 TCGA-ESCA.htseq_counts.Rdata
   11M 10  1  2021 TCGA-GBM.htseq_counts.Rdata

我们需要使用代码批量读取前面的教程：居然有如此多种癌症（是时候开启pan-cancer数据挖掘模式），的rdata文件里面的表达量矩阵，并且最简单的limma的voom算法差异分析。

rm(list=ls())
options(stringsAsFactors = F)
library(stringr)  
library(data.table)  

###### step1: 批量读取前面的Rdata  ###### 
fs=list.files('Rdata/',
              pattern = 'htseq_counts')
fs

# 首先加载每个癌症的表达量矩阵（提取非蛋白编码基因部分）
# 成为一个 list 
deg_list = lapply(fs, function(x){
  # x=fs[2]
  pro=gsub('.htseq_counts.Rdata','',x)
  print(pro)
  
  load(file =  file.path('Rdata/',x))  
  
  Group = ifelse(
    as.numeric(substring(colnames(pd_mat),14,15)) < 10,
    'tumor','normal'
  )
  table(Group)
  
  if(sum(Group=='normal') < 30){
    return(NULL)
  } else{
    library(limma)
    Group = factor(Group,levels = c('normal','tumor'))
    design = model.matrix(~ Group)
    head(design)
    v = voom(pd_mat,design,normalize="quantile")
    fit = lmFit(v,design)
    fit = eBayes(fit)
    colnames(coef(fit))
    #[1] "(Intercept)" "GroupTP"
    tT = topTable(fit, adjust="fdr", sort.by="p", number=Inf)
    head(tT)
    return(tT)
  } 
  
})
names(deg_list)
names(deg_list)=gsub('TCGA-','',gsub('.htseq_counts.Rdata','',fs))
deg_list=deg_list[!unlist(lapply(deg_list, is.null))]

因为我仅仅是挑选了33个癌症里面，那些正常对照组织样品数量超过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

定位到癌基因列表和抑癌基因列表

NCG是一个肿瘤驱动基因的数据库，目前最新版本为v7.0, 网址如下

http://ncg./index.php

共收录了 591+2756 个驱动基因，分成以下两个部分

Canonical cancer drivers (591)
Candidate cancer drivers (2756)

抑癌基因（tumor suppressor genes），也称肿瘤抑制基因，这里推荐TSGene数据库，它的网址如下：

https://bioinfo./TSGene/

最新版本为v2.0版本， 1217 human TSGs (1018 protein-coding and 199 non-coding genes).

比较出名的癌基因和抑癌基因大家需要有一个认知，比如我们读取下载的两个文件, Canonical cancer drivers (591) 和 1217 human TSGs ；

library(data.table)
drivers = fread('canonical_drivers.txt',data.table = F)[,1]
tsg = fread('Human_TSGs.txt',data.table = F)[,2]
length(intersect(drivers,tsg))
head(intersect(drivers,tsg))

可以看到，它们首先就有交集了：

> length(intersect(drivers,tsg))
[1] 174
> head(intersect(drivers,tsg))
[1] "AMER1"    "APC"      "ARHGEF12"
[4] "ARID1A"   "ARID2"    "ASXL1"   
>

批量把表达量上下调基因和癌基因以及抑癌基因做交集

这里定义每个癌症的统计学显著的表达量上下调基因，我使用了固定阈值， logFC_cutoff = 1 以及 adj.P.Val<0.05，所以得到基因数量会比较多：

sum1 = lapply(names(deg_list), function(x){
  deg_limma_voom = deg_list[[x]]
  head(deg_limma_voom)
  logFC_cutoff <- with(deg_limma_voom,mean(abs( logFC)) + 2*sd(abs( logFC )) );logFC_cutoff
  logFC_cutoff = 1
  gene_up= rownames( deg_limma_voom[with(deg_limma_voom,logFC > logFC_cutoff & adj.P.Val<0.05),])
  gene_down=rownames( deg_limma_voom[with(deg_limma_voom,logFC < -logFC_cutoff & adj.P.Val<0.05),]) 
  
  df=data.frame(
    cancer=x,
    group=c('up','down'),
    deg=c(length(gene_up),length(gene_down)),
    inter1 = c(length(intersect(gene_up,drivers)),length(intersect(gene_down,drivers))),
    inter2 = c(length(intersect(gene_up,tsg)),length(intersect(gene_down,tsg))) 
  )
  return(df)
})
sum1 = do.call(rbind,sum1 )
sum1$ratio1 = sum1$inter1/length(drivers)
sum1$ratio2 = sum1$inter2/length(tsg)

这12个癌症的统计，如下所示：

> options(digits = 2)
> sum1
   cancer group  deg inter1 inter2 ratio1 ratio2
1    BRCA    up 2171     64    101  0.108  0.083
2    BRCA  down 2431     79    220  0.134  0.181
3    COAD    up 2467     62    115  0.105  0.094
4    COAD  down 2614     87    185  0.147  0.152
5    HNSC    up 1944     61    104  0.103  0.085
6    HNSC  down 2010     47    148  0.080  0.122
7    KIRC    up 2870     98    150  0.166  0.123
8    KIRC  down 2492     66    170  0.112  0.140
9    KIRP    up 2771     72    125  0.122  0.103
10   KIRP  down 2541     80    226  0.136  0.186
11   LIHC    up 1782     58     81  0.098  0.067
12   LIHC  down 2088     63    181  0.107  0.149
13   LUAD    up 2142     60    101  0.102  0.083
14   LUAD  down 2563     88    235  0.149  0.193
15   LUSC    up 3161    100    154  0.169  0.127
16   LUSC  down 3510    115    269  0.195  0.221
17   PRAD    up  987     21     40  0.036  0.033
18   PRAD  down 1527     43    131  0.073  0.108
19   STAD    up 1756     69    100  0.117  0.082
20   STAD  down 1932     54    161  0.092  0.132
21   THCA    up 1528     43    107  0.073  0.088
22   THCA  down 1432     49    112  0.083  0.092
23   UCEC    up 2684     74    138  0.125  0.113
24   UCEC  down 2803    105    249  0.178  0.205