癌基因是基因的一类，激活后可促使正常细胞癌变、侵袭及转移。通常我们以为的激活是它表达量上升，尤其是肿瘤部位和正常组织对比的时候，而癌基因激活的方式包括点突变、基因扩增、染色体重排、病毒感染等。
而抑癌基因（tumor suppressor genes），也称肿瘤抑制基因，或俗称抗癌基因，是一类存在于正常细胞内可抑制细胞生长并具有潜在抑癌作用的基因。抑癌基因在控制细胞生长、增殖及分化过程中起着十分重要的负调节作用，它与原癌基因相互制约，维持正负调节信号的相对稳定。当这类基因在发生突变、缺失或失活时可引起细胞恶性转化而导致肿瘤的发生。所以通常我们以为他的失活是在肿瘤部位表达量下降。
那么，很自然的一个需要探索的是癌基因一定在肿瘤部位高表达吗？以及抑癌基因一定在肿瘤部位低表达吗？
下面，我们以TCGA数据库的多种癌症来举例说明这两个问题：
批量差异分析
在前面的教程：居然有如此多种癌症（是时候开启pan-cancer数据挖掘模式），我们把全部的TCGA的33种癌症的表达量矩阵区拆分成为蛋白编码基因和非编码基因这两个不同的表达量矩阵，并且保存成为了rdata文件。
在不同癌症里面，蛋白质编码相关基因数量一直在一万八附近，而非编码基因数量跨度比较大，从一万二到两万七不等。因为非编码基因的组织病人特异性都很强，但是绝大部分的癌基因和抑癌基因都是蛋白质编码基因，所以我们仅仅是针对于每个癌症的不到两万个的蛋白质编码基因在肿瘤部位和正常组织进行转录组测序表达量矩阵的差异分析。
代码语言：javascript复制 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算法差异分析。
代码语言：javascript复制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个癌症的样品数量满足要求，每个差异分析都是得到如下所示的矩阵：
代码语言：javascript复制> 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.kcl.ac.uk/index.php
共收录了 591 2756 个驱动基因，分成以下两个部分
- Canonical cancer drivers (591)
- Candidate cancer drivers (2756)
抑癌基因（tumor suppressor genes），也称肿瘤抑制基因，这里推荐TSGene数据库，它的网址如下：
https://bioinfo.uth.edu/TSGene/
最新版本为v2.0版本， 1217 human TSGs (1018 protein-coding and 199 non-coding genes).
比较出名的癌基因和抑癌基因大家需要有一个认知，比如我们读取下载的两个文件, Canonical cancer drivers (591) 和 1217 human TSGs ；
代码语言：javascript复制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))
可以看到，它们首先就有交集了：
代码语言：javascript复制> 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，所以得到基因数量会比较多：
代码语言：javascript复制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个癌症的统计，如下所示：
代码语言：javascript复制> 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
可以看到，肉眼看上去，很难说清楚每个癌症的统计学显著的表达量上下调基因，跟癌基因以及抑癌基因的交集有什么偏好性。当然了， 如果需要用 数值说话而不是肉眼看，这里可以使用卡方检验。
我试了试切换阈值去修改每个癌症的统计学显著的表达量上下调基因数量，也试过换奇怪数据库拿到癌基因以及抑癌基因列表，但是结论都是一样的。没有很明显的倾向性。
结论显而易见
肯定是有一些癌基因在肿瘤里面高表达，但是并不是全部的癌基因哦，因为导致肿瘤发生发展只需要少量关键基因即可。
同理肯定是有一些抑癌基因在肿瘤里面表达量降低，但不会是全部。
所以，并不是说在肿瘤里面表达量升高，它就有可能是癌基因，因为也有很多抑癌基因也是跟癌症的统计学显著的表达量上调基因有交集。
生命科学就是这样的复杂。
