在:癌基因一定在肿瘤部位高表达吗 我们针对每个癌症都在各种内部做了肿瘤组织和正常对照的差异表达量分析,然后在癌基因都是肿瘤的风险因子吗 我们针对每个癌症的全部基因批量了做了单基因的cox分析。
但是前面我们一直在使用癌基因列表和抑癌基因列表,去跟每个癌症的统计学显著的表达量上下调基因,以及每个癌症的统计学显著的保护因子风险因子取交集,发现它们并没有特别的偏好性。我们一直没有直对比上下调基因和保护因子风险因子,而这个往往是大家初次接触肿瘤学的时候比较困惑的。
我就不下十次的接到了这样的求助,为什么它的目标基因在肿瘤里面高表达,但是它没有统计学显著的生存意义呢,或者为什么它的一个目标基因明明是有生存意义但是却并不会差异表达呢?甚至有一些目标基因明明是在肿瘤里面高表达,但是它居然生存分析的时候是保护因子,也就是说它表达量越高的那些病人反而死的越慢。
任意选一个癌症看表达量变化和生存情况的相关性
因为在:癌基因一定在肿瘤部位高表达吗 我们仅仅是针对normal样品数量大于30的癌症做了差异分析,所以首先需要把表达量变化和生存情况的癌症数据分析结果对齐。
然后我们选取两者共有的第一个癌症,绘制散点图:
代码语言:javascript复制cg=intersect(names(deg_list),names(cox_list))
cg
deg_list=deg_list[cg]
cox_list=cox_list[cg]
deg1=deg_list[[1]]
cox1=as.data.frame(cox_list[[1]])
ids=intersect(rownames(deg1),rownames(cox1))
df=data.frame(
deg=deg1[ids,'logFC'],
cox=cox1[ids,'HR']
)
head(df)
library(ggplot2)
library(ggpubr)
library(ggstatsplot)
p <- ggplot(df, aes(x=df$deg,
df$cox))
geom_point()
labs(x = "deg",
y = "cox")
scale_color_manual(values=c("blue", "grey","red"))
geom_vline(xintercept = c( -1.5,0,1.5),lty=4,col="#ec183b",lwd=0.8)
geom_hline(yintercept = c(-0,1,2 ),lty=4,col="#18a5ec",lwd=0.8)
#xlim(-5,5)
#ylim(-3,3)
theme_ggstatsplot()
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank())
p
可以看到:
两者基本上没有关系,在肿瘤组织里面相当于正常组织的基因表达量改变,基本上并不会影响它是否在肿瘤病人队列里面的的预后作用。
我上面展现的仅仅是一个癌症看表达量变化和生存情况的相关性,如果大家认真学习了我们的笔记,很容易对全部的癌症做一个批量处理,而且上面的相关性散点图也很容易加上线性回归曲线及方程式。那么作为一个练习题吧,大家把全部的癌症做完线性回归方程式。
九宫格看表达量变化和生存情况的交集
看表达量变化和生存情况的相关性我在横纵坐标加上了两条虚线,作为判断统计学显著的表达量上下调基因的阈值,我使用了固定阈值, logFC_cutoff = 1.5 以及 adj.P.Val<0.05,所以得到表达量上下调基因数量会比较多。但是里面没有引入P值的筛选,需要修正。
接下来我们就加入p值筛选,然后判断一下 表达量变化和生存情况的交集的九宫格里面的每个基因数量。代码如下所示:
代码语言:javascript复制head(deg1)
deg1$change=ifelse(deg1$adj.P.Val>0.05,'stable',
ifelse( deg1$logFC > 1.5 ,'up',
ifelse( deg1$logFC < -1.5,'down','stable') )
)
table(deg1$change)
cox1$sur=ifelse(cox1$p>0.05,'stable',
ifelse( cox1$HR > 1 ,'protective', 'risk')
)
table(cox1$sur)
ids=intersect(rownames(deg1),rownames(cox1))
df=data.frame(
deg=deg1[ids,'change'],
cox=cox1[ids,'sur']
)
table(df)
可以看到:
代码语言:javascript复制 cox
deg protective risk stable
down 174 9 1293
stable 580 646 14474
up 38 58 995
在肿瘤组织里面表达量下调的基因,我们通常是会以为是抑癌基因,它确实倾向于是保护因子。
在肿瘤组织里面表达量上调的基因,我们通常是认为是癌基因,它也确实是倾向于是风险因子。
参考:
- https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6915928/
批量检查全部的癌症的表达量变化和生存情况的交集
也是很简单的循环即可:
代码语言:javascript复制sum = do.call(rbind,
lapply(names(deg_list), function(x){
deg1=deg_list[[x]]
cox1=as.data.frame(cox_list[[x]])
head(deg1)
deg1$change=ifelse(deg1$adj.P.Val>0.05,'stable',
ifelse( deg1$logFC > 1.5 ,'up',
ifelse( deg1$logFC < -1.5,'down','stable') )
)
table(deg1$change)
cox1$sur=ifelse(cox1$p>0.05,'stable',
ifelse( cox1$HR > 1 ,'protective', 'risk')
)
table(cox1$sur)
ids=intersect(rownames(deg1),rownames(cox1))
df=data.frame(
deg=deg1[ids,'change'],
cox=cox1[ids,'sur']
)
tmp = as.data.frame(table(df))
tmp$cancer = x
tmp
}))
head(sum)
library(reshape2)
dcast(sum, deg cox ~ cancer,value.var ="Freq" )
如下所示的结果:
代码语言:javascript复制> dcast(sum, deg cox ~ cancer,value.var ="Freq" )
deg cox BRCA COAD HNSC KIRC KIRP LIHC LUAD LUSC PRAD STAD THCA UCEC
1 down protective 174 43 71 356 64 116 464 0 3 1 1 21
2 down risk 9 45 75 351 482 40 3 609 3 194 87 484
3 down stable 1293 1536 1099 927 1198 1060 1083 1857 600 902 700 1295
4 stable protective 580 42 498 4575 947 160 1374 56 66 230 185 368
5 stable risk 646 441 2613 2438 1151 5228 401 1371 143 835 445 3395
6 stable stable 14474 14512 12897 7943 12747 10292 13762 12658 16831 15968 15640 11214
7 up protective 38 10 28 174 109 3 42 39 2 33 48 82
8 up risk 58 113 141 645 232 438 258 46 9 53 4 389
9 up stable 995 1097 831 890 1148 440 870 1775 393 779 778 1047
可以看到在各个癌症里面了,up基因都倾向于是risk基因,而不是protective基因,而down基因都倾向于是protective基因,而不是risk基因。
有意思的是, LUAD 和 LUSC 呈现出来了完全相反的结果,我都怀疑是不是自己代码有问题,因为是批量处理,仓促写笔记教程。(也有可能是固定阈值惹的祸,其实有很多很好玩的现象,但是我们的时间精力都很有限,而且我们还得竭尽全力的去追逐sci文章,渐渐地放弃自己的兴趣爱好。)
很简单的结论
大家的感觉是正确的,那些在肿瘤里面异常高表达的基因,确实容易是风险基因,它表达量越高病人死的越快。
但是大数据告诉我们,也有很多基因虽然在肿瘤里面异常高表达,但是没有生存意义,或者说是反而是保护因子,这个该如何解释呢?
其实从逻辑推理来说,表达量的高低判定是肿瘤组织和正常组织的差异分析,但是生存分析仅仅是肿瘤内部的,两者不应该是有数学层面的联系。但是因为肿瘤组织高表达的通常是我们常规意义的癌基因,而癌基因大多就是风险基因,所以它们就这样在生物学层面被关联起来了。