转录是遗传信息从DNA传递给RNA的过程,是中心法则的第一步。转录组在组学发展史上一直都有着举足轻重的作用,也是目前应用最广的高通量测序分析技术之一。
作为转录组学研究的基础,基于表达量的分析是理解生物学过程的关键,也是最为经典的一种分析手段。在这个专栏中,我们将深入探讨转录组学的三种基于表达量的分析方法,分别是差异分析、趋势分析和WGCNA分析。今天我们的首篇文章将聚焦于转录组的差异分析,这一经典方法通过两两比较样本之间的基因表达水平,来识别关键的生物学变化,帮助我们缩小范围,寻找目标基因及通路。
基本概念
在转录组数据分析中,差异分析是最常见的分析方法。差异分析的核心在于评估不同组别(处理组)之间的表达差异是否显著大于组内(随机误差)的差异。
通过分析表达量的变化倍数(Fold change)和差异检验值(P/FDR值),可以确定哪些基因在样本间或组间表现出显著的差异,并将这些结果与后续的富集分析等方法相结合,以探究相关的生物学过程和分子机制。
转录组差异分析主要流程
基因表达量统计:
计算每个基因或转录本的表达量,通常采用原始 reads count 、FPKM(每千碱基每百万读段)、TPM(每百万转录本数)等方法进行展示。
原始 reads count 表示该转录本所包含reads数目,但受测序量和基因长度的影响。为保证后续分析准确性,我们先对测序深度进行校正,再对基因或转录本的长度进行校正,获得标准化后的FPKM值或TPM值,再进行后续分析。
差异表达基因筛选:
使用统计方法(如DESeq2、edgeR等)比较不同样本间的基因表达水平,筛选出差异表达的基因。常用的筛选指标包括差异倍数(FC)和差异检验值p值/FDR。
FC(Fold Change):
FC表示两个样品中同一个基因表达水平的变化倍数。它反映了实验组和对照组之间的表达量差异。FC值越大,表示表达量的变化越显著。在差异表达分析中,通常会设定一个log2FC阈值,例如大于1或2,来筛选差异表达基因。
p值:
p值是统计学中用来表示观察到的数据与假设的理论数据差异是否显著的一个指标。p值越小,表示基因的表达差异越显著。然而,p值本身容易受到多重假设检验的影响,因此在实际应用中需要结合其他指标进行综合判断。
FDR(False Discovery Rate):
FDR用于控制假阳性结果的比例,即在所有被错误认为有显著差异的基因中,真正有差异的基因所占的比例。在转录组分析中,由于需要对大量基因进行假设检验,如果不进行FDR校正,即使单次检验的假阳性率很低,累积的假阳性率也可能很高。FDR的计算通常采用Benjamini-Hochberg方法,通过对每个基因的p值进行排序和校正,从而控制最终的假阳性率。
一般情况下筛选 FDR<0.05 且 |log2FC|>1 的基因为显著差异基因。也可以根据具体情况适当放宽或缩紧筛选标准。
常用差异分析可视化图形
经过以上的一系列操作,我们终于筛选出来一部分真正具有显著差异的基因,但这些基因的数量可能是巨大的,因此我们需要利用直观的图形将数据可视化,以便更好地阐述生物学问题。
常用的图形有柱状图,火山图,热图、韦恩图等,更加高端一点的图形有小提琴图、多组差异散点图、雷达图等等,他们所能展示的细节会更加丰富。在文章中常常采用以上多种图形组合的方式来多角度展示差异分析结果。
01
柱状图
柱状图,也称条形图,是转录组类型文章中最基础,也是高频出现的经典图形。它适用于比较两个或多个组别的数据(如不同时间或条件),通常只有一个变量,常用于小数据集分析。
通过柱状图,我们可以直观地看到组内组间上调和下调的差异基因数量,从而反映出生物过程的大致趋势,也可以帮助我们验证一些预期的实验结果。
02
火山图和多组火山图
火山图是一类用来展示比较组间差异基因表达情况的图形。该图形本质是散点图的一种,它能将统计显著性数值(-log10(FDR))和差异倍数(log2FC)相结合,从而能够帮助人们识别那些变化幅度较大且具有统计学意义的基因。
虽然火山图所展示的细节很丰富,但它的局限性在于只能比较两个组的数据。而下面要介绍的这个图形弥补了这个缺陷。
多组差异散点图(也叫多组火山图)可以一次性展示多个比较组的差异基因分布情况。它的横坐标为比较组的名称,纵坐标为log2FC,绘制时会自动筛选Pvalue/FDR<0.05的差异显著的基因进行展示。
03
MA图
MA图(Minus-versus-Add)是另一种用来展示差异分析的散点图,在早期microarray(芯片转录组测序)数据分析中会被使用。
有时,研究人员需要利用点图比较两组表达量的差异,X轴是两组数值的均值,即(a+b)2,Y轴是两组的差值,即b-a,一个是add,一个是minus,这就是MA(Minus-versus-Add)plot了。
根据上述公式,M就变成了log2(FC),A是组间表达量的平均强度,给差异显著(一般为差异倍数为2倍以上,且p值小于0.05)的点不同的颜色,显示两组样品差异显著性。
当一个点(即一个基因)的Y值是0,说明它在两组间没有差别,当它X轴数值越大,说明它在两组的均值越大。那么,当一个点的X轴数值很大,Y轴绝对值也很大的时候,就说明它是那种平均表达量高,组间差别还很大的基因,意味着这个基因一定在其中一组有着惊人的表达量。反过来,如果它Y轴绝对值很大,但是X轴数值很小,说明它很有可能是小量表达的基因,微量的变化带来较大的倍数波动。
04
聚类热图
在差异基因分析完成后,我们可以将感兴趣地基因用聚类热图来展示。每个小方格表示每个基因在不同样本中的表达量情况,表达量越大颜色越深。上方和左侧的树形图表示对不同样品和基因的聚类分析结果。
热图可以直观呈现多样本多个基因的全局表达量变化,并展示出不同样本的基因表达趋势。除此之外,我们还可以根据聚类的结果,观察样本的重复性,看是否有不同组的样本被聚类到一起,还可以从基因聚类的角度,看哪些基因具有比较一致的表达变化,或许能够在功能预测上给我们一些启发。
在绘制聚类热图时,我们不可避免地要思考一个问题,就是归一化的问题,归一化简单来说就是将一组数据通过归一化处理,使其符合均值为0,方差为1的标准正态分布。那么我们是选择归一化的形式呢?
a)按行均一化:
将每一行数值分别单独处理,使其符合标准正态分布。通常我们是以基因为单位来观测这些表达量数值的变化。例如,A基因表达量从10变化到20,B基因表达量从100变化到200,我们更关心它们变化的倍数,这时我们应该选择以基因为单位进行归一化。
因此,按基因归一化处理,可以最大程度地呈现每个基因的变化信息,避免一个超高表达的基因掩盖了其他基因的变化。在绘制热图时,它是最常用的一种归一化策略。
b)按列均一化:
将每一列数值分别单独处理,使其符合标准正态分布。如果我们想通过基因的表达量来观测样本的重复性好坏或样本分类,这时可以通过按列归一化处理。
c)行和列都均一化:
将所有的行列数据一起处理,使其符合标准正态分布。如果想让高表达的基因对样本的分类起到更大的作用,这时可以选择对所有数值归一化。
d)不做归一化处理:
表示只想通过颜色的渐变来观察基因的表达量变化,可以选择不做归一化处理。
05
韦恩图和Upset图
维恩图和upset图是都基于组与组之间比较得到的差异基因进行的集合,重叠部分即是不同处理下各组样本中被共同调控的基因集,单独的部分则是某种处理下特定调控的基因集。基于维恩图或upset图我们可以分别对共有或特有基因集进行深度挖掘。
韦恩图相信大家已经不陌生了,也比较好理解。但是韦恩图一般只能基于2-5个比较组来展示,更多的比较组所做的韦恩图就比较繁杂不美观了。而Upset图是升级版的韦恩图,它可以更好地展示6个及以上的比较组的交并集情况。
06
雷达图
雷达图是以从同一点开始的轴上表示的三个或更多个定量变量的二维图表的形式显示多变量数据的图形方法。图上可以一次性展示信息,包括log2(FC)值、平均表达量和上下调情况等。我们可以展示两组样本间根据P值或者Q值整体展示差异程度最大的TOP n基因的信息。
转录组差异分析应用
讲了这么多理论,接下来我们就要开始“实战”了,也就是学习如何将差异分析运用到我们自己的研究中。在看完这篇文章后,你定会豁然开朗,意识到看似简单的差异分析也能被应用地如此巧妙!
这篇来自江苏科技大学的基迪奥客户的文章,在2024年2月发表在Plant Physiology and Biochemistry(IF:6.5)杂志上。作者通过RNA测序技术对桑树(Morus alba)叶片在不同浓度硼处理条件下的转录组进行了差异分析,旨在揭示桑树应对硼胁迫的分子机制。实验一共设置了5个处理,包括缺硼(T1)、中度缺硼(T2)、适量硼(CK)、硼中毒(T3)和严重硼中毒(T4),其中CK为对照。
作者在第一部分首先展示了转录组数据的质量以及样本的序列信息,包括reads过滤情况、主成分分析(PCA)、样本相关性热图、FPKM 密度图、样本小提琴表达图、样本表达韦恩图。
图1 转录组数据的序列统计
在第二部分,作者开始对转录组数据进行整体上的差异分析。将4个实验组与CK进行两两差异分析,并筛选FDR<0.05且|log2FC|≥1的基因为差异表达基因(DEGs),通过柱状图、韦恩图以及火山图对这几千个差异基因进行展示。
图2 不同硼水平下桑树的上下调差异表达基因(DEGs)
紧接着,作者对筛选出来的DEGs进行了GO和KEGG富集分析,分别依据p值和富集因子(enrichment factor)来绘制热图,p值代表了富集的显著程度,p<0.05(红色)代表显著富集。而富集因子代表了该term/pathway对应的差异表达基因数目与基因组背景中富集到这一term/pathway的数目的比值,比值越高富集程度越显著。
图3 DEGs的GO和KEGG聚类热图
在第三部分,作者开始聚焦于那些与植物对硼胁迫有着密切响应的功能基因,包括转录因子家族、转运蛋白家族、光合作用相关基因、植物激素和信号传导相关的差异表达基因、与细胞壁结构相关的差异表达基因和抗氧化相关基因等。
对于筛选出来的每一类基因,作者都会用柱状图与聚类热图统一而直观地展示了它们在不同处理条件下的表达模式。使用差异倍数(log2FC)绘制的柱状图可以直观地观察每种基因的差异上下调情况,而使用聚类热图则可以方便展示基因在不同硼处理下的表达模式差异。
图4 部分功能基因对不同硼处理的响应展示
在第四部分,作者利用散点图展示了16个基因进行RT-qPCR验证的结果。结果显示RNA-seq与RT-qPCR数据之间存在显著的相关性(R2=0.88),证明了RNA-Seq结果的高度可靠性。
图5 采用qRT-PCR方法对桑树16个基因表达的RNA-seq数据进行回归分析验证。
最后,作者结合之前对不同硼处理下的生理生化以及代谢组的研究,总结出来以下结论并绘制出了一些硼处理后的关键代谢途径并梳理出完整的基因网络图。
小结:
这是一篇经典的差异分析+富集分析来完成基因表达量分析的文章,这些分析方法和图形的运用,为我们提供了一个非常全面的基因表达差异分析范例,同时也给我们带来了如下启示:
※ 多种差异分析结合使用:通过结合多种差异分析方法和图形,灵活调节筛选阈值,可以更全面地揭示基因表达的差异模式。
※ 重点关注基因家族和通路:根据之前的表型实验以及文献调研,选择重点关注的基因家族和通路进行分析,可以更准确地揭示分子机制。
※ 系统性的分析思路:采用系统性的分析思路,从基因表达差异到转录因子、信号通路等多个层面进行深入分析。
掌握了一大堆理论知识,但是对不会R语言的生信小白来说,一上手还是困难重重。没关系,登录我们的OmicsShare,使用里面的差异分析小工具,不会编程也可以分分钟搞定!点击下面的文章链接领取教程吧:
如何轻松完成差异基因分析?(qq.com)
工具链接:
https://www.omicshare.com/tools/Home/Soft/diffanalysis
如果是在我们基迪奥做转录组项目的小伙伴,我们的结题报告里包含差异分析结果。如果还想自己对数据进行挖掘,还可以在Omicsmart平台上,不用自己导入数据,直接生成个性化差异分析结果及可视化图形,随意调节筛选参数,可以说是超级无敌方便了!
下一期我们将重点介绍转录组基于表达量的分析方法二——趋势分析,这种分析方式不同于常见的差异分析,更加“高大上”,对于存在着浓度梯度或者时间梯度的样品,差异分析可以让你的文章锦上添花!大家记得关注我们的公众号,千万不要错过哦~