美女教授带你从统计学视角看转录组分析

引言

分子生物学的中心法则自1958年由Francis Crick提出到今年正好60周年，它描述了“DNA制造RNA，RNA制造蛋白质”的遗传信息的标准流程 [1]。十年前，第二代RNA测序技术(RNA-seq)的诞生及其迅速发展使得研究者可以在对RNA序列没有任何先验信息的情况下高通量地对全转录组进行测序 [2]。现如今第二代RNA测序技术已经成为了研究基因和RNA表达最常用的手段之一，它的广泛应用极大地促进了生物和医学领域的各类研究，包括对基因表达与调控，RNA可变剪切以及蛋白质翻译等多项生物过程的了解 [3]。具体见生信老司机以中心法则为主线讲解组学技术的应用和生信分析心得。

为了回答各种生物问题，十年来不同领域的研究者已为第二代RNA-seq数据分析提出了超过2000种计算与分析方法 (39个转录组分析工具，120种组合评估(转录组分析工具哪家强-导读版))。近日，加州大学洛杉矶分校（UCLA）的李婧翌（Jingyi Jessica Li）教授和她的学生李维（Wei Vivian Li）第一次从统计建模与分析的角度对第二代RNA测序数据的计算方法进行了总结和讨论，发表在最新一期的Quantitative Biology期刊中（Modeling and analysis of RNA-seq data：a review from a statistical perspective）4。该综述文章从四个层面（样本，基因，转录本，和外显子选择）对RNA-seq数据的分析方法进行了总结，旨在归纳看似不同的方法背后共通的统计假设和模型。生信宝典之前的总结在 转录组分析的正确姿势。

样本分析层面：样本相似性度量

作者首先从样本分析层面上讨论了如何利用RNA-seq数据对来自同种或者不同物种的多种细胞类型的样本进行比较，从而研究基因表达机制在物种进化过程中的分化和保守现象。这就需要用到样品间的相似性度量来检测异常样品，进行样品聚类、样品分类和新细胞类型分析。

Pearson和Spearman相关系数是最常用的计算相似性的度量，配合热图展示。但这种通过相关性推测转录组相似性的方式容易受持家基因的影响，并且相关性测量依赖于表达值的精确计算，信噪比低时结果稳定性差。

作者提出了TROM (Transcriptome overlap measure)的计算方式，采用”associated genes“而不是全部基因计算相关性。associated genes是指那些在样品间归一化得到的Z-score的值高于系统检测的阈值的基因。z-score计算见R语言 - 热图美化。作者也提供了R包TROM进行计算, 获得associated genes和进行overlap test，获得格式类似于相关性矩阵的TROM得分矩阵，用于后续分析。

另外还有其他度量方式，如偏相关系数 (partial correlation coefficient)指消除批次效应或其它影响因素后的相关性；

考虑到RNA-seq样品之间并不是线性相关的，信息论中的互信息(mutual information)可以计算非线性相关性。其它类似的还有条件互信息(conditional mutual information)。这些主要用在基因调控网路的构建上。具体见https://rdrr.io/cran/synRNASeqNet/man/parMIEstimate.html和。

直接降维的手段，如PCA、t-SNE和MDS, 见之前的文章：(一文看懂PCA主成分分析, 还在用PCA降维？快学学大牛最爱的t-SNE算法吧, 附Python/R代码)。

基因层面：基因表达动力学

其次，在基因分析层面上，作者着重讨论了两种最常见的基因表达数据分析：基因差异表达分析和基因共表达网络分析。

基因差异表达分析通过统计学的假设检验理论来研究哪些基因在不同条件下(比如实验组和对照组)存在差异表达。常用的数据标准化方法如FPKM/RPKM/TPM，是一种定量方式，一般可比，但也会有protocol特异的偏好。一般不用于差异基因检测。

基于分布的标准化方式：如DESEq, DESeq2, TMM和quantile normalization等，基于不同样品基因表达水平分布相似的原则进行计算。具体计算如下（代码比文字易于理解）。

基于基因的标准化方式是使不同样品间非差异基因或持家基因的表达水平一直，如PoissonSeq。

基于最广泛使用的Negative Binomial模型概括了检验基因差异表达的6个核心步骤：估算基因表达和离散度，估算样品间差异统计量，推导H0下的统计量，计算每个基因的统计了，计算P-value，多重假设检验校正。操作见DESeq2差异基因分析和批次效应移除。

基于log转换后的基因表达符合正太分布的模型，如voom (线性模型分析多因素设计，定量加权解释表达变化，经验贝叶斯估算)和sleuth(log转换的基因表达作为线性模型的响应变量，拆解为3个成分：样品间差异、生物自然噪音、推测噪音，通过bootsrap-test估计零假设和计算p-value)。

注意点：重复次数，FDR，没有方法是最优的，时间序列数据使用maSigPro或multivariate empirical Bayes statistic (MB statistic)。

基因共表达网络分析: 建立基因的共表达关系，推测未知基因的功能。常用的有WGCNA，检测共表达基因簇和模块。具体见WGCNA分析，简单全面的最新教程。

配合Cytoscape使用更佳：

• Cytoscape教程1
• Cytoscape之操作界面介绍
• 新出炉的Cytoscape视频教程
• Cytoscape制作带bar图和pie图节点的网络图

转录本分析层面： 重构和量化全长转录本

在转录本分析层面上，大部分研究的重心在于通过重建和量化全长转录本来研究可变剪切在细胞分化和疾病发展等过程中的作用。这两项任务（重建和量化）对于包括人在内的复杂生物体来说仍然非常困难，因为第二代测序技术一般产生的是不超过300个碱基（ bp）的短片段(reads), 但人体内全长转录本的平均长度在1700 bp左右，而短片段导致的信息损失只能依靠统计建模来推断和弥补。作者将现有的转录本量化方法分为两类：基于极大似然分析或基于回归分析。前者将转录本表达量作为参数建立混合概率模型(mixture model)，通过最大化似然函数来求解最优参数估计。后者将转录本表达量作为系数建立回归模型，通过数据拟合来求解最优系数估计。这里面涉及的软件有Cufflinks, Stringtie, Salmon, Kallisto, RSEM等，其性能评估见39个转录组分析工具，120种组合评估(转录组分析工具大比拼 （完整翻译版）)。

外显子分析层面：计算单个外显子的PSI

在外显子分析层面上，研究重心在于量化单个外显子在可变剪切的过程中存留在全长转录本中的可能性 (percentage spliced in, PSI)，是一个较为可靠的可变剪接分析方式。常用工具有MISO和rMATS。基因组浏览器在此有重要的应用，可视化的reads分布模式对应研究不同可变剪接是必须的，从图谱看到差异，再设计工具寻找差异。

• 测序数据可视化 (一)
• IGV基因组浏览器可视化高通量测序数据
• 高通量数据分析必备-基因组浏览器使用介绍 - 1
• 高通量数据分析必备-基因组浏览器使用介绍 - 2
• 高通量数据分析必备-基因组浏览器使用介绍 - 3

文章的最后作者对RNA-seq统计建模过程中仍然存在的难点进行了总结，并简要讨论了RNA-seq在RNA编辑和非编码RNA等相关问题中的应用。近年来兴起的单细胞RNA测序技术将转录组研究提高到单细胞精度，也给新数据的统计建模带来了新的机遇和挑战，可以参考The Human Cell Atlas White Paper了解更多。

本文主要是帮助师姐推荐宣传Quantitive Biology杂志，如对这篇文章原文有任何疑问欢迎联系QB微信号：13269084698，或直接在文末留言，QB将邀请作者为您解答。

附上两个完整的分析流程

转录组分析培训班第4期

在广大粉丝的期待下，《生信宝典》联合《宏基因组》在2018年11月16-18日北京鼓楼推出《转录组数据分析》专题培训第四期，为大家提供一条走进生信大门的捷径、为同行提供一个转录组实战分析学习和交流的机会、助力学员真正理解分析原理和完成实战分析，独创四段式教学(3天集中授课+自行练习2周+再集中讲解答疑+上课视频回看反复练习)，“教—练—答—用”四个环节统一协调，真正实现独立分析大数据。

关于学习生物信息学分析的重要性，请阅读《生物信息9天速成班—成为团队中不可或缺的人》。

本课程一共3天，每天6节课，共18节课，全部课程均理论与实战结合(只要课上讲的都是可以带你自己实现的分析)。从分析平台搭建、Linux和R基础、图表解读和实战、转录组设计、分析标准流程、差异基因分析、功能富集分析、及各类高级分析(差异剪接、WGCNA分析、通路图绘制等)，和CNS级图片修改排版。3天时间，老司机带您完成自学需要3个月甚至是1年的崎岖之路，助力您真正玩转转录组分析。

课程大纲

每节课1小时一个主题，理论结合实战，学懂原理，实战实操，全是老司机多年经验和代码的无私分享。下面是课程安排，如11代表第一天第一节课，26代表第二天第六节课，41为两周后的线上集中视频答疑。这次课程重点做了调整，侧重后期高级分析，弱化流程性分析。

教程内容简介如下：

转录组的应用、设计和案例分享

1. 转录组学研究技术介绍
2. 转录组学实验设计和测序原则、注意事项
3. 二代、三代测序过程和原理解析
4. 转录组学文章案例分析
5. 在线基因表达资源数据库

转录组分析流程实战

1. 转录组分析流程评估
2. 测序数据质量评估和清洗
3. 不基于比对的差异基因分析
4. 基于比对的差异基因分析
5. 转录本组装和选择性剪接分析
6. 目标基因GSEA/GO富集分析

常见图表解读和图形编辑排版

在培训上，结合发表高水平文章，进一步讲解16种常用分析图的原理和使用范围，让你不仅读懂图，更知道如何应用于自己的研究，并亲自轻松完成绘图。

针对大家使用R语言绘图学习时间成本较高的问题，易生信团队针对常用16种图开发了免费绘图网站，一键出图，更可鼠标点选参数修改图形的个性样式。

成果发表是科研过程中不可缺的一部分，发表成果又少不了图形展示。文章图表排版是否整齐规范、协调一致、重点突出对一篇文章的发表也是有不少贡献的。之前推出的文章发表图的修改和排版讲演了部分图形编辑和排版操作，本次培训也会实践从原始图形、到细节修饰再到排版发表的整个过程和注意事项。

基因组浏览器用于多组学数据的可视化和关联分析，本地有IGV，在线有UCSC genome Browser和Epigenomebrowser，各有特色。

转录组高级分析

1. WGCNA基因共表达分析
2. WGCNA基因、表型关联分析
3. Cytoscape 共表达网络绘制
4. 转录组常见图形在线绘制
5. KEGG/Reactome通路图绘制，表达映射
6. 基因互作的文献挖掘和数据库挖掘展示