|
作者：酸菜 转载请注明：解螺旋·临床医生科研成长平台 话说天下大势，粉久必冷。为了温暖你们的小心脏，今日祭出“不做实验却想发SCI”的一篇生信+统计的文章教学。是不是看到题目就兴奋的不要不要的？ 不做实验的理由千千万，核心一条就是没钱。医生做科研，做实验，学统计，敲代码，你总得会一样。统计学深了那是数学问题，酸菜自从高考数学差一分及格之后就受到了永久性的心理创伤，用SPSS做做单因素和多因素分析已到极限，剩下的就只有生信这一条路来尝试无本生意了。 今天的教学案例2017年新发表于Oncotarget （IF = 5.168），一本让人爱恨交织的杂志，一年收稿三千多篇，没发过的人都说它很水，内心却是弱水三千，吾亦想取一瓢的挣扎纠结。 能发一篇Oncotarget不错了，要啥自行车。 这篇文章的创新在于揭示了EB病毒基因与人类胃癌转录组的交互作用。EB病毒同胃癌的关系文献早有报道，研究分子机制是常规思路，而用生物信息学、统计学方法从TCGA （The Cancer Genome Atlas） 数据库挖掘到特征性分子标记进行EB病毒与人基因的交互网络分析，则是新颖的切入点，这个模式借鉴一下，微生物+癌症基因组合，潜在的资源颇丰。 文章总体思路是先找到EBV(+)和EBV(-)的胃癌中有差异表达的基因和通路，随后分析它们和EBV基因的相关性，构建交互作用网络。技术路线分三大步骤，作者用了一张图来概括： Figure4 数据分析流程 第一步是数据特征提取（A），第二步是单因素相关分析（B），第三步是多因素相关分析（C）。 单因素和多因素分析是临床研究的入门统计技术了，关键性难度在于测序数据分析方面，基因组数据本来就复杂，这里又需要分析人和病毒两个转录组，其中的多重检验会导致很多的假阳性率，所以需要降维和校正。 第一步：提取特征分子标志 从TCGA数据库下载EBV(+)和EBV(-)的原始测序fastq数据，有25个EBV(+)和260个EBV(-)样本。然后构建一个临床特征表格，包括性别、年龄、肿瘤部位、组织学病理学诊断等等，接着用R语言的dist()函数计算两组样本间临床特征的欧式距离，找出距离最近的样本进行匹配。 因为有几个EBV(-)样本跟EBV(+)样本的距离是一样的，所以最终选到了20个EBV(-)样本，来跟那25个EBV(+)样本配对。这样每个配对之间的距离，就比25个EBV(+)对260个EBV(-)之间的平均距离小多了，控制了混杂因素，使组间具有可比性。 该文章附件中的表格展示了配对样本的一部分，左边4列分别为EBV(+)样本编号、组织学诊断、EBV(-)样本编号、组织学诊断。右边2列数字就是欧式距离，左边的是EBV(+)样本与它所匹配的EBV(-)样本的距离，右边的是该EBV(+)样本与原260个EBV(-)样本的平均距离，可见前者是比后者要小多，消除了样本间差异的偏倚因素。 接下来对EBV基因进行筛选，从原始的88个基因中，将在超过5个样本中原始计数（raw count）为0的基因去除，剩下19个基因，这就获得了EBV相关的特征分子标志。 部分EBV基因在EBV(+)样本中的表达 然后提取人类基因组的特征分子标志。用Bioconductor的DEseq2包找出两组样本间差异表达的基因（DEx gene），这里就要用到Bonferroni法对P值进行校正，原始P值达到1E-6，即校正后的P < 0.05，才算有统计学意义。共找到939个DEx基因，其中189个上调，750个下调。 再用Bioconductor中的GAGE (Generally Applicable Gene-set Enrichment for Pathway Analysis)分析法，找到两组样本中差异表达的通路，共计29个，在EBV(+)样本中全是上调的。 同时用R语言的MEGENA包构建那些DEx基因的共表达网络，找到它们的聚类基因模块（Gene modules）和其中的关键节点基因（hub genes）。前者modules是指一组具有较强相关性的基因，而hub则是指modules中具有全局性影响的枢纽基因。这一步是利用人类基因组中本来就有的交互网络数据，来减少数据维度，聚焦到关键信息。这一步找到了91个节点基因和27个基因模块 (一个节点基因可出现在不同的模块中)。 然而，这些数据维度还是高，怎么进一步抓取核心要素呢？作者又用了主成分分析法（PCA）提取其中起主要作用的信息，他们只提取了第一个主成分（PC1）来代表该模块或通路。这样，网络构建和相关性分析的素材提取完毕，最终获得19个EBV基因、91个人类胃癌相关的hub基因、27个基因模块和29个通路（的PC1）。 对于生信有基础的同学，可按照上述方法利用工具自行探索，零基础的同学我们还有后续教学。 第二步，单因素相关分析 单因素相关分析就是要从19个EBV基因、91个人类胃癌相关的hub基因、27个基因模块和29个通路中再进一步获取相关性最高的组合，采用配对的Pearson相关性检验对三组数据分别进行分析：1）基因模块的PC1和EBV基因表达数据；2）差异通路的PC1和EBV基因表达数据；3）hub基因和EBV基因表达数据。用置换检验法获取相关系数（用C++）。 至于相关系数的显著性，则没有采用通用的p<0.05，因为它假设两个变量都是正态分布，这对RNA-seq数据来说通常不现实。所以研究者采用p < 0.1 作为统计学显著性的门槛。 这一步得到了7个EBV基因和12个模块的PC1呈显著相关，其中相关性最强的是BLAF1基因和12号模块，相关系数达到0.662，如下表。 3个EBV基因和4个人类胃癌的hub基因相关，其中最强的是LMP-1和Clorf115，相关系数0.754。 又有2个EBV基因和4个通路的PC1相关，其中BLAF4和4个通路都相关，最强的是BLAF4和磷脂酰肌醇信号通路，相关系数达到0.709。 好了，单因素分析这步到这里做完。 第三步，多因素关联分析 这步要用到稀疏典型相关分析法（sparse Canonical Correlation Analysis, sCCA）。典型相关分析法（CCA）是分析两个数据矩阵的经典方法，不需要降维提取主成分。但普通CCA只适用于行比列多的矩阵，即样本比基因多，这情况显然不对，是样本少基因多，所以要用sCCA。这步用R语言的PMA包来做。 其中用CCA()函数可得到每个元素（比如每个EBV基因）的典型系数（canonical coefficient），它表示该元素对整个矩阵典型相关系数的贡献，如果某个元素的典型系数非零，则表示它对全局相关性有很大贡献，被选为核心（essential）基因。 这步得到22个基因模块与19个EBV基因计数矩阵相关（P< 0.1，理由同上），如下表。其中核心EBV基因为LMP-1，BALF1，BALF2，BARF1，BNRF1，LF1和BZLF1。 此外，人类hub基因计数矩阵和EBV基因计数矩阵的典型相关系数是0.806。对这个相关性贡献最大的核心EBV基因和核心人类hub基因就不一一列出了。所有的差异表达的通路都和EBV基因计数矩阵达到了显著相关的典型系数。最强几个的如下表所示。 至此，数据分析的主要步骤就做完了。但文章没完呀，还要做个总结，比较一下Pearson和sCCA两种相关性分析的结果。咱们接着看。 第四步，数据汇总 先总结一下EBV基因的情况。如下表，各个重要EBV基因在Pearson和sCCA分析中，分别与人类胃癌基因模块、hub基因、通路相关的情况。
人类重要基因模块在两种相关性分析中的情况则如下表：
比较有意思的是模块5和模块12。模块5在DAVID数据库的注释中，有膜、跨膜、离子通道等关键词；而模块12则和饥饿激素通路相关，在前人的研究中，它可能会促进胃肠道和胰腺疾病的恶化。由此，文章数据分析的统计学结果终于与生物学意义产生了交汇！ 下图展示了模块5和EBV基因相互作用的关系图，作为上述相关性的一个示例：
再来看看重要的人类hub基因，它们是C1orf115，CNTD2和VANGL2。它们在Pearson分析中，和EBV的BALF1和LMP-1相关，在sCCA分析中也被选为核心基因。下图展示了人类hub基因和EBV基因的互作网络。
人类重要的通路则是凋亡、溶酶体、Jak-STAT信号通路和磷脂酰肌醇信号通路，它们在Pearson和sCCA分析中都与EBV的BALF4和/或BALF5基因相关。值得注意的是，BALF4基因与Jak-STAT信号通路、磷脂酰肌醇信号通路中的大多数基因都呈现正相关。 总体而言，在第一步差异表达分析中所找到的DEx基因、模块和通路，与筛选出来的EBV基因，大多都有相关性，这也验证了前人的观察和实验的结果，即EBV和胃癌确实存在基因层面的相互作用。 生信和统计分析完成后，得到了一堆具有统计学相关性的基因、通路，这时候需要找出它们的生物学相关性，才是拔高文章水平的点睛之笔。作者分析了他们找到的每个具有统计学相关性的元素在既往文献报道中的相关性，文献检索的工作相当细致。例如LMP-1编码的潜伏膜蛋白1可导致细胞生长失调，是一个致癌蛋白；BALF1编码抗凋亡蛋白Bcl2同源体，它的表达也和几种癌症的发生有相关性，等等等等。 当然除了验证前人的研究，BALF2可谓一个新发现，之前没有报道过它和癌症有什么关系，本文是第一次提出，如果能够经过实验验证这文章绝非五分这般平庸啊。 TCGA的挖掘潜力日益受到科研界的重视，搞肿瘤的这一优势得天独厚，应善加利用。生信的技能学习是一个由浅入深的渐进过程，解螺旋将推出系列课程，今天放出第一波，由我们的合作伙伴GCBI制作的TCGA数据库基础应用教学，课件免费下载，在服务号完成分享任务后获取16mins视频。 懒癌晚期患者，请加入解螺旋钻石会员，不但能够学习临床研究，基金申请和SCI写作的深度课程，当下还可享受“TCGA基因万能筛”活动，1000元，给你用生信方式筛选一个靠谱的研究靶标基因，实验验证不成功免费重筛，把专业的事交给专业的人办。 下面是“TCGA基因万能筛”活动获得的靶标基因示例： 服务特色 针对研究癌症却没有方向的会员，整理TCGA数据库中特定癌症的癌/癌旁数据，包括miRNA或mRNA数据，筛选癌/癌旁差异基因，并结合多种生物信息学分析内容（包括差异基因功能、信号通路、生存以及基因网络分析等），为会员研究特定癌症的发生发展机制提供有价值的靶基因及相应的理论基础。 分析结果展示 差异分析结果
差异基因功能分析
差异基因信号通路分析 基于KEGG数据库，对差异基因利用Fisher精确检验和卡方检验，把目标基因参与的Pathway进行显著性分析，按照P value<0.05进行筛选，得到显著性的Pathway。
生存分析 利用单基因在所有样本中的表达情况，按照中位数进行分组，利用时序检验(log-rank)检验进行比较两组的总生存期的差异性，筛选条件为：p<0.05，研究基因表达对预后的影响。
调控关系 每个基因与基因之间都有非常多的调控关系，我们从基因相关的转录因子、miRNA、lncRNA和上下游相关基因的四个角度来展示基因与基因间的关系。
基因CFH中，已经验证的上游靶向结合的miRNA为has-miR-146a-5p。在基因网络中，CFH能够抑制基因CFB和C3的表达。
文献报道中和CFH相关的lncRNA的71个。
转录因子预测 根据位置结合关系，预测结合CFH启动子区域（转录起始位点上游2000bp，下游500bp）的转录因子。并依据Transfac数据库的预测分值和COSMIC数据库以及dbSNP数据库是否存在对应注释信息整合出推荐度；以推荐度来表示预测出来的转录因子的科研价值，推荐度越高越好。
|
|
|
