在：癌基因一定在肿瘤部位高表达吗 我们探索发现并不是使用的癌基因都在肿瘤部位高表达，也不是所有的高表达基因都是癌基因，对抑癌基因也是如此。这个就很有意思的，因为癌基因的定义就是那些在肿瘤里面过度激活的基因，而抑癌基因就是在肿瘤里面失活的基因，不过过度激活不一定要转录本大量增加，可能是其它生物学机理，比如蛋白质产物大量增加，又或者说蛋白质产物效果增强。
同理，我们会问另外一个问题，就是癌基因都是肿瘤的风险因子吗，它高表达会导致癌症比如死的越来越快吗？反之，抑癌基因一定是肿瘤的保护因子吗，它表达量越高癌症病人越受到保护吗，因为想当然的我们会认为抑癌基因能抑制癌症嘛，所以它表达量越高越好。
同样的，我们可以使用TCGA数据库的多种癌症来举例说明这两个问题：
整理表达量矩阵和生存分析
首先，我们选择同样的 TCGA-CDR-SupplementalTableS1.xlsx 文件里面的生存信息，在前面的教程给出来了下载地址：
代码语言：javascript复制library(readxl)
phe=read_excel('../MC3/TCGA-CDR-SupplementalTableS1.xlsx',sheet = 1)
survdata=as.data.frame(phe[,c(2,3,26:27)] )
survdata$OS.time = survdata$OS.time/365
boxplot(survdata$OS.time)
survdata=survdata[survdata$OS.time>0,]
table(survdata$type)
survdata=na.omit(survdata)
table(survdata$type)
可以看到，绝大部分癌症类型里面的有生存信息的肿瘤病人数量还可以：
代码语言：javascript复制> table(survdata$type)
ACC BLCA BRCA CESC CHOL COAD
91 410 1083 294 45 442
DLBC ESCA GBM HNSC KICH KIRC
47 185 595 527 112 535
KIRP LAML LGG LIHC LUAD LUSC
288 173 511 371 509 497
MESO OV PAAD PCPG PRAD READ
85 582 184 179 500 164
SARC SKCM STAD TGCT THCA THYM
261 454 417 134 506 123
UCEC UCS UVM
546 56 80
接下来就是把表达量矩阵信息跟这个临床信息相对应，这里我们仍然是 使用代码批量读取前面的教程：居然有如此多种癌症（是时候开启pan-cancer数据挖掘模式），的rdata文件里面的表达量矩阵，然后根据TCGA数据库的病人的ID来对表达量矩阵信息跟这个临床信息相对应。
代码语言：javascript复制
fs=list.files('Rdata/',
pattern = 'htseq_counts')
fs
# 首先加载每个癌症的表达量矩阵
# 在每个癌症里面与临床信息相对应。
pre_sur = lapply(fs, function(x){
# x=fs[9]
pro= gsub('TCGA-','',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)
exprSet <- pd_mat[,Group=='tumor']
dim(exprSet)
exprSet[1:4,1:4]
pid=substring(colnames(exprSet),1,12)
exprSet=exprSet[,!duplicated(colnames(exprSet))]
colnames(exprSet) = substring(colnames(exprSet),1,12)
# 1. prepare data for coxph----
## 批量生存分析使用coxph回归方法
## http://www.sthda.com/english/wiki/cox-proportional-hazards-model
phe <- survdata[survdata$type==pro,c(1,4,3)]
colnames(phe)=c('p','time', 'event')
phe=phe[!duplicated(phe$p),]
rownames(phe)=phe$p
ids= intersect( rownames(phe), colnames(exprSet) )
phe=phe[ids,]
exprSet=exprSet[,ids]
return( list(
phe=phe,
exprSet=exprSet
))
})
names(pre_sur) = gsub('TCGA-','',gsub('.htseq_counts.Rdata','',fs))
do.call(rbind,
lapply(pre_sur, function(x){table(x$phe$event)}))
可以看到每个癌症里面的病人的结局时间分布并不固定：
代码语言：javascript复制> do.call(rbind,
lapply(pre_sur, function(x){table(x$phe$event)}))
0 1
ACC 51 28
BLCA 228 178
BRCA 925 151
CESC 220 71
CHOL 18 18
COAD 340 98
DLBC 38 9
ESCA 97 64
接下来进行批量单基因cox分析
必须要保证前面的表达量矩阵信息跟这个临床信息相对应哦，这样我的代码才有可能是复制粘贴就能运行：
代码语言：javascript复制library(survminer)
library(survival)
cox_list = lapply( pre_sur , function(x){
phe=x$phe
exprSet=x$exprSet
mySurv <- with(phe, Surv(time, event))
cox_results <-apply(exprSet , 1 , function(gene){
# gene= as.numeric(exprSet[1,])
group=ifelse(gene>median(gene),'high','low')
if( length(table(group))<2)
return(NULL)
survival_dat <- data.frame(group=group,# stage=phe$stage,
stringsAsFactors = F)
m=coxph(mySurv ~ group,
# mySurv ~ stage group, # 如果有交叉变量
data = survival_dat)
beta <- coef(m)
se <- sqrt(diag(vcov(m)))
HR <- exp(beta)
HRse <- HR * se
#summary(m)
tmp <- round(cbind(coef = beta, se = se, z = beta/se, p = 1 - pchisq((beta/se)^2, 1),
HR = HR, HRse = HRse,
HRz = (HR - 1) / HRse, HRp = 1 - pchisq(((HR - 1)/HRse)^2, 1),
HRCILL = exp(beta - qnorm(.975, 0, 1) * se),
HRCIUL = exp(beta qnorm(.975, 0, 1) * se)), 3)
return(tmp['grouplow',])
})
return(t(cox_results))
})
names(cox_list)
names(cox_list)=gsub('TCGA-','',gsub('.htseq_counts.Rdata','',fs))
save(cox_list,file = 'batch_cox_list.Rdata')
我们简单的看看其中一个癌症的结果矩阵：
代码语言：javascript复制> load(file = 'batch_cox_list.Rdata')
> head(cox_list[[1]])
coef se z p HR HRse HRz HRp HRCILL HRCIUL
ZZZ3 -0.470 0.38 -1.23 0.220 0.62 0.24 -1.57 0.12 0.29 1.32
ZZEF1 -0.421 0.38 -1.10 0.273 0.66 0.25 -1.36 0.17 0.31 1.39
ZYX -0.224 0.38 -0.58 0.559 0.80 0.31 -0.66 0.51 0.38 1.70
ZYG11B -0.051 0.38 -0.13 0.894 0.95 0.36 -0.14 0.89 0.45 2.01
ZYG11A -1.235 0.40 -3.12 0.002 0.29 0.12 -6.17 0.00 0.13 0.63
ZXDC -0.212 0.38 -0.56 0.577 0.81 0.31 -0.62 0.54 0.38 1.70
要根据这个批量单基因cox分析结果来挑选保护因子和风险因子，跟前面的使用了固定阈值， logFC_cutoff = 1 以及 adj.P.Val<0.05，确定每个癌症的统计学显著的表达量上下调基因是一回事。
我们仍然是需要首先筛选p值确定统计学显著，然后根据HR值来判断保护因子和风险因子：
关于HR值：
- HR = 1: No effect
- HR < 1: Reduction in the hazard
- HR > 1: Increase in Hazard
Note that in cancer studies:
- A covariate with hazard ratio > 1 (i.e.: b > 0) is called bad prognostic factor
- A covariate with hazard ratio < 1 (i.e.: b < 0) is called good prognostic factor
值得注意的是理论上HR > 1是增加了肿瘤的风险，意思是该基因表达量升高会增加肿瘤风险，但是我上面的代码里面选取的是 return(tmp['grouplow',]) ，所以我们的结果HR > 1的含义是该基因表达量降低会增加肿瘤风险。
解释起来会有一点点绕，不过这样的结果很容易通过一个分组KM曲线去肉眼检查一下基因具体到底是保护因子还是风险因子，就跟我们肉眼检查表达量上下调基因会使用箱线图一样。
跟表达量上下调基因
跟前面的笔记：癌基因一定在肿瘤部位高表达吗 的代码大同小异：
代码语言：javascript复制load(file = 'batch_cox_list.Rdata')
head(cox_list[[1]])
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))
sum1 = lapply(names(cox_list), function(x){
# x= names(cox_list)[[11]]
print(x)
cox_results = as.data.frame(cox_list[[x]])
if(nrow(cox_results)<10){
return(NULL)
}else{
head(cox_results)
cox_results=cox_results[cox_results$p < 0.05 ,]
gene_up= rownames( cox_results[ cox_results$HR > 1 ,])
gene_down= rownames( cox_results[ cox_results$HR < 1 ,])
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)
options(digits = 2)
sum1
结果还是蛮让我意外的，如下所示：
代码语言：javascript复制> sum1
cancer group deg inter1 inter2 ratio1 ratio2
1 ACC up 606 14 37 0.0237 0.03040
2 ACC down 4338 150 219 0.2542 0.17995
3 BLCA up 2464 67 103 0.1136 0.08463
4 BLCA down 696 14 60 0.0237 0.04930
5 BRCA up 775 33 57 0.0559 0.04684
6 BRCA down 700 16 36 0.0271 0.02958
7 CESC up 2141 71 101 0.1203 0.08299
8 CESC down 245 5 20 0.0085 0.01643
9 CHOL up 427 19 21 0.0322 0.01726
10 CHOL down 103 3 7 0.0051 0.00575
11 COAD up 95 3 5 0.0051 0.00411
12 COAD down 586 14 21 0.0237 0.01726
13 DLBC up 79 1 5 0.0017 0.00411
14 DLBC down 73 1 1 0.0017 0.00082
因为每个癌症都是使用了同样的阈值，所以各自的保护因子和风险因子数量不一样。但是可以看到，跟前面的笔记：癌基因一定在肿瘤部位高表达吗 的结论类似，并没有明显的倾向性。
肯定是有一些癌基因，它在肿瘤里面高表达却不影响生存，但是并不是全部的癌基因哦，因为导致肿瘤发生发展只需要少量关键基因即可。
同理肯定是有一些抑癌基因在肿瘤里面表达量降低，但也没有影响生存。
其实生存分析受到了的干扰因素非常多，一个目标基因可能是非常有临床意义所以它有统计学显著的生存意义，但是两万多个基因总是有那么一些基因跟目标基因表达量相关性非常高，所以也有 统计学显著的生存意义 。但是我们确实不能因为某个基因有统计学显著的生存意义就判定他一定是很重要。
比如，我们随机选取几个基因看看：
代码语言：javascript复制> head(cox_list[[1]][cox_list[[1]][,4]<0.05,])
coef se z p HR HRse HRz HRp HRCILL HRCIUL
ZYG11A -1.24 0.40 -3.1 0.002 0.291 0.115 -6.2 0.000 0.134 0.63
ZWINT -2.38 0.54 -4.4 0.000 0.092 0.050 -18.1 0.000 0.032 0.27
ZWILCH -1.75 0.47 -3.8 0.000 0.174 0.081 -10.2 0.000 0.070 0.43
ZW10 -1.11 0.42 -2.6 0.008 0.329 0.138 -4.9 0.000 0.145 0.75
ZSWIM4 -0.83 0.40 -2.1 0.037 0.436 0.174 -3.2 0.001 0.200 0.95
ZSWIM3 0.92 0.41 2.3 0.024 2.505 1.016 1.5 0.139 1.131 5.55
批量绘制多个基因的KM曲线代码是：
代码语言：javascript复制
## 针对第一个癌症，具体检查几个基因
head(cox_list[[1]][cox_list[[1]][,4]<0.05,])
cg = head(rownames(cox_list[[1]][cox_list[[1]][,4]<0.05,]))
cg
x=pre_sur[[1]]
phe=x$phe
exprSet=x$exprSet
identical(colnames(exprSet), rownames(phe))
mySurv <- with(phe, Surv(time, event))
library(ggstatsplot)
overlap_list <- lapply(cg, function(i){
# i = cg[1]
survival_dat = phe
gene = as.numeric(exprSet[i,])
survival_dat$gene = ifelse(gene > median(gene),'high','low')
table(survival_dat$gene)
library(survival)
fit <- survfit(Surv(time, event) ~ gene,
data = survival_dat)
survp = ggsurvplot(fit,data = survival_dat, #这里很关键，不然会报错
legend.title = i, #定义图例的名称
# legend = "top",#图例位置
# legend.labs = c('High', 'Low'),
pval = T, #在图上添加log rank检验的p值
# pval.method = TRUE,#添加p值的检验方法
risk.table = TRUE,
risk.table.y.text = F,
xlab = "Time in years", #x轴标题
# xlim = c(0, 10), #展示x轴的范围
# break.time.by = 1, #x轴间隔
size = 1.5, #线条大小
ggtheme = theme_ggstatsplot(),
palette="nejm", #配色
)
return(survp)
})
n = length(overlap_list)
n
x=floor(n/2);y=2
all_plot <- arrange_ggsurvplots(overlap_list,print = F,ncol =x, nrow = y,
risk.table.height = 0.3,
surv.plot.height = 0.7)
all_plot
出图如下：
批量绘制多个基因的KM曲线代码
我在生信技能树多次分享过生存分析的细节；
- 人人都可以学会生存分析（学徒数据挖掘）
- 学徒数据挖掘之谁说生存分析一定要按照表达量中位值或者平均值分组呢？
- 基因表达量高低分组的cox和连续变量cox回归计算的HR值差异太大?
- 学徒作业-两个基因突变联合看生存效应
- TCGA数据库里面你的基因生存分析不显著那就TMA吧
- 对“不同数据来源的生存分析比较”的补充说明
- 批量cox生存分析结果也可以火山图可视化
- 既然可以看感兴趣基因的生存情况，当然就可以批量做完全部基因的生存分析
- 多测试几个数据集生存效应应该是可以找到统计学显著的！
- 我不相信kmplot这个网页工具的结果(生存分析免费做)
- 为什么不用TCGA数据库来看感兴趣基因的生存情况
- 200块的代码我的学徒免费送给你，GSVA和生存分析
- 集思广益-生存分析可以随心所欲根据表达量分组吗
- 生存分析时间点问题
- 寻找生存分析的最佳基因表达分组阈值
- apply家族函数和for循环还是有区别的（批量生存分析出图bug）
- TCGA数据库生存分析的网页工具哪家强
- KM生存曲线经logRNA检验后也可以计算HR值
生存分析是目前肿瘤等疾病研究领域的点睛之笔！
