一个基因引发的血案
大家好，我是老米，学习生信一个月，这是我的第二篇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的样本。
过程不详述，大概如下：
从USCS下载了TCGA的前列腺癌基因表达数据（550个样本）和临床病人数据。
批量处理了Normal组和Tumor组表达量，计算他们的差异，并取t.test P < 0.05的基因（得到1万3千多个基因）。
再和临床数据merge，批量分析统计OS P < 0.05的基因，最终得到586个表达有显著差异且OS显著的基因(Gepia 才十几个[捂脸])。
用refGenome 处理GTF数据对基因类别进行了注释。大致看了下，然后挑了个lncRNA等待实验验证。哈哈，二区文章，我来啦。。。
总结：数据库有风险，挖掘需谨慎！新年快乐！
最后按照惯例，我们应该是有生信技能树的友情推广，但是我们长沙站招生已经满了，而且未来的两个月并没有外出巡讲计划，所以就不宣传了，大家可以继续添加小助手，说出自己的城市需求，我们会优先安排明年巡讲场次：广州专场（全年无休）GEO数据挖掘课，带你飞（1.11-1.12）和 生信入门课全国巡讲2019收官--长沙站