一个基因引发的血案

大家好,我是老米,学习生信一个月,这是我的第二篇Markdown。不知道多少人还记得我的第一个作品:原来一个星期真的可以零基础入门TCGA数据挖掘,甚至markdown写作公众号投稿 (感兴趣的自己点击查看哦)

今天给大家讲个鬼故事。呃,这大过年的。。。最近翻到技能树的一个帖子:中国制造:碉堡的TCGA可视化网站GEPIA。该文章才发表两年,基于Gepia数据挖掘发表了很多文章:

分数都不低,二区文章一大堆。

这么强大的工具,忍不住在其网站逛了好久。。。各各基因之间跳来跳去。。。琢磨着看能不能找到几个有表达及生存都有差异的基因?。。。然后突然想到一个伟大的想法,能不能把所有这些基因都选出来?该死的,网站不提供查询结果的下载,怎么办? 一个个网页复制粘贴?NONNO,网络上有一堆强大的爬虫工具,比如八爪鱼,后羿…安装之后倒腾了好久,终于设定好爬虫流程,妈呀,采集效率太慢了,初步估算(后来证实)也要两天时间才能爬完!刚好又是星期五,静静的开着电脑让他爬吧。可是晚上回到家忍不住,实在是想快点拿到数据,突然想到之前偷学了点python爬虫,能不能改改拿来用用?说干就干,花了一个小时调整好了python代码。16万条数据,十来分钟就下好了。python爬虫效率真高!我要学python!

喂,快醒醒,跑题了!
说回正题,拿到了自己想要的数据,马上就进行了匹配处理,代码如下:

rm(list = ls())
survSigDat <- read.csv("../GEPIA-survival.csv",sep = ",",header = T)
geneExpr <- read.table("../All.txt",sep = "\t",header = F)

colnames(geneExpr) <-c("CancerType","GeneName","GeneID","Median_Tumor","Median_Normal","Log2FC","adjp")
geneExpr$GeneName <- gsub(" ","",geneExpr$GeneName)
geneExpr$GeneID <- gsub(" ","",geneExpr$GeneID)

library(dplyr)
geneAndSurvSig <- left_join(geneExpr,survSigDat,by= c("CancerType","GeneID","GeneName"))
geneAndSurvSig <- geneAndSurvSig[(which(geneAndSurvSig$pValueOS != "")),]

geneAndSurvSig_prad <- geneAndSurvSig[(which(geneAndSurvSig$CancerType == "PRAD")),]
write.table(geneAndSurvSig_prad,file = "geneAndSurvSig_prad.txt")

大概查看了下前列腺癌中,表达有差异且癌症病人中位值生存分析有差异的基因。结果如下:

看起来不错嘛,不过好像都是在PRAD里面下调的。pubmed上查一下,剩下没报道的应该就是新课题了吧!哈哈哈。。。

不查不得了,一查吓一跳!SNHG12 这个基因刚刚被发表了:

进文章里面一看,姨?好像不对啊!作者说这个基因在癌症里是上调的,和Gepia出来的结果不相同啊。
SNHG12文章图:

Gepia上面是说下调的:

一个说上调,一个说下调。肿么会这样?
细看,SNHG12作者的Normal组为N=52,Gepia的Normal组为151(是不是这100个正常人的前列腺有毒)。Tumor组差别不大,SNHG12为497,Gepia为492。而Gepia生存分析用到的也是490个病人。

看来可能是因为Normal组引起的差异。查看了Gepia数据库的来源:

除了TCGA,Gepia还采用了GTEx的样本,GTEx是个啥?各种器官组织捐献者的的表达数据库。而SNHG12文章只用了的TCGA样本。

那前列腺癌里,到底这个基因上调还是下调?  再看一下oncomine (Oncomine 整合了GEO、TCGA和已发表的文献来源的RNA和DNA-seq数据),不同文章对这个基因的表达情况?结果如下:

大牛们的数据基本支持上调的趋势!

独秀同学继续提问,为啥加入GTEx正常标本之后,表达差异不一样了?

为此,独秀同学google了一下,发现了这个文献:Unifying cancer and normal RNA sequencing data from different sources。没细读,文章大意是:当不同数据源进行差异比较的时候(如对GTEx和TCGA数据),不能直接拿来比,需要用某种方法进行uniform进行标准化,balabalabala…

文章作者最后说我们成功的整合了GTEx和TCGA的数据,现在可以做比较了!

这篇文章发表在2018年,后于Gepia的开发。我猜测目前的Gepia版本没做normalization,所以直接加入GTEx数据进入Normal组还真的是有毒!或许,Gepia可以出version 2.0了?

如果真是这样,那些基于Gepia数据挖掘发了文章的童鞋们会不会有些心慌慌,哈哈哈。。。

那么问题来了, 今后该选哪个数据库进行数据挖掘? 可靠吗?会不会经常见鬼?

鬼故事讲完了,回到最开始那个伟大的想法。用Gepia数据可能有出入【捂脸】,只能自己动手写代码实现了,统一用TCGA的样本。

过程不详述,大概如下:

  1. 从USCS下载了TCGA的前列腺癌基因表达数据(550个样本)和临床病人数据。

  2. 批量处理了Normal组和Tumor组表达量,计算他们的差异,并取t.test P < 0.05的基因(得到1万3千多个基因)。

  3. 再和临床数据merge,批量分析统计OS P < 0.05的基因,最终得到586个表达有显著差异且OS显著的基因(Gepia 才十几个[捂脸])。

  4. 用refGenome 处理GTF数据对基因类别进行了注释。大致看了下,然后挑了个lncRNA等待实验验证。哈哈,二区文章,我来啦。。。

总结:数据库有风险,挖掘需谨慎!新年快乐!

最后按照惯例,我们应该是有生信技能树的友情推广,但是我们长沙站招生已经满了,而且未来的两个月并没有外出巡讲计划,所以就不宣传了,大家可以继续添加小助手,说出自己的城市需求,我们会优先安排明年巡讲场次:广州专场(全年无休)GEO数据挖掘课,带你飞(1.11-1.12)和 生信入门课全国巡讲2019收官--长沙站

(0)

相关推荐