1.快速筛出候选基因

转录组测序可以帮助研究者全局且快速地对某一生命过程中所有发生变化的基因进行筛查，然而只有从全局的基因变化分析到定位具体基因的讨论，研究者才能有的放矢开展后续的基因表型验证和机制探索。我们对快速筛出候选基因提供以下通用性的建议和方法，您也可以结合其他筛选指标一起筛选出关键基因。

1.1 筛选思路

候选基因应该与您所研究的生物学问题密切相关，并在某一生命过程中发挥核心的生物学功能、参与重要的通路。基因功能的背景信息可以来源于参考基因组注释、基因功能数据库（如GO、KEGG）注释或者相关研究的文献和综述等，对于大且复杂的研究方向，比如“纤维化”、“植物盐胁迫”等，建议通过文献综述（可以在PubMed检索，https://pubmed.ncbi.nlm.nih.gov/）获取可能相关的基因、GO和KEGG通路信息，然后在自己的测序结果中进行筛选。无论通过怎样的方式进行数据筛选，人为的解读和挑选是必不可少的。有时候差异基因过多，还是不得不用肉眼去挖掘每个基因所代表的意义。这个工作不可省略，我们不得不去做这个繁琐但必要的工作！很多大课题组筛选数据也没有其他技巧，都是挨个儿看挨个敲除每个基因来验证。如果暂时没有思路，可以考虑一些调控上游的因素，如转录因子、组蛋白修饰、DNA结合蛋白、RNA结合蛋白等。

1.2 筛选方法

1）基于GO富集分析结果

打开GO富集结果*_GO_enrichment.xlsx文件，在GO_Term列中，用关键词逐一筛选与您研究相关的生物学功能，通常设置pvalue小于0.05。获得相关显著富集的GO Term后建议在基因差异分析表格*_Gene_differential_expression.xlsx中的GO列进行筛选（使用GO ID “GO:0030335”或GO Term “positive regulation of cell migration”进行筛选，由于很多GO Term名称高度相似，建议优先使用GO ID进行筛选），可以得到注释到相关GO Term的所有基因，再筛选significant列显示为“yes”的即可获得注释到相关GO Term且差异显著的基因。由于*_Gene_differential_expression.xlsx包含所有基因的表达量信息和差异分析信息（差异倍数、p值、q值），因此在*_Gene_differential_expression.xlsx中进行筛选可以获得某GO Term中的基因表达全貌，有利于捕获那些不满足默认差异阈值条件（|log2FC|>=1 & q<0.05）但同样具有统计学差异（q<0.05或p<0.05）的基因。示例如下：

（使用Excel或WPS的筛选功能，从GO富集结果文件中筛选“positive regulation of cell migration”功能，找到11个差异表达的基因，随后可以在*_Gene_differential_expression.xlsx中的GO列筛选获得这11个差异表达基因的信息。

详情请点击：

1. 玩转GO和KEGG富集因子图的N种姿势：3种数据处理（含在线筛选条目），3种排序方式，本地交互图片  

2. 跟着小姐姐一起学习如何做GO和KEGG基因功能富集分析 | 生信快闪

3.【疯狂7月】⑬—KEGG通路标色：表达量/上下调、单组/多组 | 云平台

（使用Excel或WPS的筛选功能，从KEGG通路富集结果文件中筛选“PI3K-Akt signaling pathway”通路，有17个差异表达基因与“PI3K-Akt signaling pathway”相关，随后可以在*_Gene_differential_expression.xlsx中的KEGG列筛选获得这17个差异表达基因的信息。）

请注意：在某些实验处理条件下，差异表达基因可能较少，因而GO和KEGG通路富集结果较少，如果已经有关注的GO或KEGG通路，可直接从*_Gene_differential_expression.xlsx中的GO列或KEGG列进行筛选，获取候选基因信息。

详情请点击：

1. 玩转GO和KEGG富集因子图的N种姿势：3种数据处理（含在线筛选条目），3种排序方式，本地交互图片  

2. 跟着小姐姐一起学习如何做GO和KEGG基因功能富集分析 | 生信快闪

3.【疯狂7月】⑬—KEGG通路标色：表达量/上下调、单组/多组 | 云平台

GSEA从基因集的富集角度出发，理论上更容易囊括细微但协调性的变化对生物通路的影响。我们可以基于GSEA分析结果筛选在实验组或对照组中富集的基因集（GO条目、KEGG通路等），然后可以关注CORE ENRICHMENT（对ES值有主要贡献的基因），即Leading edge subset。

详情请点击：

1.三分钟绘制一张优美的GSEA图 | 云平台

2.常规分析找不到差异基因？GSEA分析了解一下 | 学习专栏

综合上述两种方法，初步确定候选的基因范围。接下来通过以下条件，进一步缩小候选基因范围（仅供参考）。

表达丰度：受系统噪音影响，低表达丰度基因（FPKM值<10）的差异倍数并不可靠，建议选择表达丰度中等水平的基因进行后续研究。

基因新旧：在PubMed上查询候选基因的相关文章数量，建议选择较新的基因（文献数量<100）开展深入研究，请注意同一个基因常存在多个别名。

分子大小：考虑到后续功能验证实验是否能开展，建议选择0.5~2.5K大小的基因。

除了上述建议，您还可以关注测序结果里面的特殊基因，比如相关修饰酶、转录因子信息，其中动物转录因子可以前往AnimalTFDB（http://bioinfo.life.hust.edu.cn/AnimalTFDB#!/）查询，植物转录因子可以前往PlantTFDB（http://planttfdb.gao-lab.org/）查询。如果您研究的物种在AnimalTFDB和PlantTFDB中没有收录，可以通过序列比对近源物种或模式物种转录因子的方式获得可能的转录因子信息。如果有类似需求，可以将需求（附上项目编号）和需要参考的物种信息反馈至售后技术支持邮箱（support@lc-bio.com）。

通过蛋白质-蛋白质相互作用（PPI）网络分析可以在STRING数据库中查询关注基因或基因集的相关性蛋白，以及哪些蛋白处于网络的核心。

详情请点击：

云课堂(17) | Cytoscape内插stringAPP完成PPI蛋白互作分析指南

韦恩图（Venn），也叫温氏图、维恩图，是用于显示集合重叠区域的关系图。我们可以用韦恩图筛选不同比较组（不同基因集）共同的元素和特有元素，比如共同差异表达基因、特异差异表达基因等，并将这种筛选过程和结果可视化。灵活利用韦恩图可以通过一次筛选或多次筛选获取我们关注的核心基因集，然后可以对核心基因集进行进一步分析。比如利用韦恩图筛选某基因敲减组差异上调基因和过表达组差异下调基因的交集，韦恩图详细方法介绍参见：

1.小姐姐告诉你如何画韦恩图 | 生信快闪

2. OmicStudio开挂绘图 | UpSet图原理讲解+实操演示（扫描下方二维码）

       当样本为时间节点样本、浓度梯度样本、治疗或用药前后设置时，若进行两两比较筛选，当样本组较多时可能需要多次的相交筛选才能获得目标变化趋势的基因集。基于STEM（http://www.cs.cmu.edu/~jernst/stem/）的基因表达模式聚类分析（趋势分析）可以快速获得我们关注趋势的基因集或者哪种变化趋势是最显著的变化趋势。OmicStudio已经上线表达模式分析工具（https://www.omicstudio.cn/tool/37），输入基因表达量数据即可进行表达模式分析。表达模式分析需要将所有组（或自己关注的组）按照时间节点排序、浓度梯度排序或对照、模型、治疗的顺序排序，输入的是基因表达量（比如FPKM值）信息，每组需要使用一个表达量表征该组的基因表达量，通常使用表达量平均值（对该组所有样本取表达量平均值）。关于OmicStudio上的表达模式分析可以参考：

小姐姐教你使用STEM进行趋势分析

加权共表达网络分析，即WGCNA，旨在寻找协同表达的基因模块（module），并探索基因网络（模块）与关注的表型之间的关联关系以及网络中的核心基因（Hub基因）。WGCNA适用于复杂的数据模式，如果样本数大于15（建议5组以上），可以考虑开展WGCNA分析。关于WGCNA介绍可以参考：

高分文章中都在用的找核心基因的神器—WGCNA你确定不了解下吗？

我们默认是对差异显著基因进行富集分析，在我们更改差异显著阈值、韦恩图筛选、趋势分析和WGCNA分析后可以获得相关关注基因集信息，随后可以对关注基因集进行富集分析，以获悉与这些筛选基因相关的GO或通路。OmicStudio已上线通用版富集分析工具，可以适用于几乎所有物种。

详情请点击：

1. 玩转GO和KEGG富集因子图的N种姿势：3种数据处理（含在线筛选条目），3种排序方式，本地交互图片  

2. 跟着小姐姐一起学习如何做GO和KEGG基因功能富集分析 | 生信快闪

3.【疯狂7月】⑬—KEGG通路标色：表达量/上下调、单组/多组 | 云平台

2. qPCR验证候选基因

 开展后续实验前，通常会使用qPCR对4.2中的候选基因进行定量检测，以验证是否与测序结果（变化趋势）一致。由于不同检测技术具有不同的检测误差，建议选择同时满足以下条件的差异表达基因进行qPCR验证。

 1）兴趣是第一位的。选择自己感兴趣的基因进行验证，即使其表达倍数可能没有到达2倍（比如1.5倍）。如果没有明确感兴趣的基因，可以参考以下几条；

 2）优先选择表达量变化幅度较大的基因，比如差异倍数 ≥ 2，即fc ≥ 2 or ≤ 0.5，或是log2(fc) ≥ 1 or ≤ -1；

 3）优先选择具有统计学差异的基因，比如qvalue < 0.05；如果满足此q值条件的基因较少，可以使用pvalue < 0.05进行筛选，对于基因群体而言，使用q值可以降低假阳性；对于某基因个体而言，如果p值小于0.05，也可以考虑进行后续验证与功能机制探索。如果自己的样本类型异质性较高且在进行测序时生物学重复数偏低，即使p值大于0.05，在其表达量尚可且变化幅度尚可的情况下也是可以考虑的；

 4）优先选择表达量比较高的基因，在比较的样本组中，至少有一组样本的基因表达丰度（FPKM值）平均值大于10（不是绝对标准）；

 5）建议优先选择生物学重复样本间一致性较高的基因（对于异质性较高的样本，基因在生物学重复样本间可能也有明显差异，但差异倍数是实验组和对照组表达量平均值相除，一些异常表达量可能对差异倍数产生影响）。

注意：

 1）由于测序结果只对测序样本负责，因此qPCR验证测序结果，请使用测序同批样本；

 2）您只需对候选基因进行qPCR验证，我们不建议您对非候选的差异表达基因或非差异表达基因进行qPCR验证，因为测序技术已较为成熟，不必要对技术可靠性进行验证；

 3）考虑到qPCR验证测序结果的吻合率通常为60~70%，如可能，建议选5~10个或更多的候选基因进行验证；

 4）如果样本类型异质性比较高，比如临床样本、动物模型、组织液等，那么注意验证样本的选择或者使用更多的样本（远大于3）进行qPCR验证，以判断某基因是否在两组中存在统计学差异。异质性较高的样本的组内表达平均值容易受偏离值或异常值的影响，从而影响差异倍数的计算，因此对于此种情况，不要只看差异倍数和p值，也应该重点关注组内表达量的一致性或两组间的整体差异；

5）低表达基因（比如所有样本中FPKM值小于1）无论是测序定量和qPCR定量，结果可靠性都会降低，因此不建议选择低表达基因；如果您对其感兴趣、并且使用多对引物（针对不同外显子设计引物）在多于三次重复实验的基础上得到一致的结果，以您的qPCR结果为准。

6）由于可变剪切的存在，一个基因可能存在多个转录本，针对于不同外显子设计引物的qPCR结果可能不同（甚至相反）。基因定量是基因座所有转录产物的综合定量结果，在设计引物时可以考虑引物设计在大部分转录本的共有序列区间。当然如果对某条mRNA（比如经典转录形式）或lncRNA感兴趣，引物尽可能设计在其特有序列区间。

3.常见问题

A：fastq是一个文本格式，用于储存生物学序列及其相应质量值（通常是核酸序列的）。为了方便储存及可读这些信息，这些序列以及质量信息使用ASCII字符标示。该格式最初由Sanger开发，目的是将FASTA序列与质量数据放到一起，目前已经成为高通量测序结果的事实标准。通常fastq文件中每一个序列含有4行信息，第一行：以‘@’开头，是这一条read的名字，这个字符串是根据测序时的状态信息转换过来的，中间不会有空格，它是每一条read的唯一标识符，同一份FASTQ文件中不会重复出现，甚至不同的FASTQ文件里也不会有重复；第二行：表示序列信息，制表符或者空格不允许出现。一般是明确的DNA或者RNA字符，由A，C，G，T和N这五种字母构成，N代表的是测序时那些无法被识别出来的碱基；第三行：用于将测序序列和质量值内容分离开来。以‘+’开头，后面是描述信息等，或者什么也不加。在旧版的FASTQ文件中会直接重复第一行的信息，但现在一般什么也不加（节省存储空间）；第四行：测序read的质量值，每个字符与第二行的碱基一一对应，按照一定规则转换为碱基质量得分，进而反映该碱基的错误率，因此字符数必须和第二行保持一致，它描述的是每个测序碱基的可靠程度，用ASCII码表示。第四行中每个字符对应的ASCII值减去64，即为对应第二行碱基的测序质量值。如果测序错误率用E表示，Illumina的碱基质量值用Qphred表示，则有下列关系：Qphred=-10log10E。

A：FPKM值和物种、不同的发育生长时期都是有关系的，没有一个绝对和统一的标准。根据文献中的报道，一般认为0.5以下的受测序误差影响比较大，可以认为是没有表达，0.5-1属于低丰度表达基因，1-10属于中低丰度表达基因，10-50以上属于中高丰度表达基因，50以上属于高丰度表达基因，该标准仅供参考。

A：MD5全称是Message-Digest Algorithm 5，是一种被广泛使用的密码散列函数，可以产生出一个128位（16字节）的散列值（hash value），用于确保信息传输完整一致。MD5值等同于文件的ID，它的值是唯一的。如果文件内容（不是文件名）发生改变，那么MD5值会发生变化。由于二代测序原始数据一般较大，为了保证数据在传输过程没有发生损坏，可以使用MD5值进行文件校验。Windows平台可以使用notepad++生成文件的MD5值（工具>>>MD5>>>从文件生成...），Mac平台可以打开终端，定位到文件的位置（比如如果文件在桌面，可以使用指令cd Desktop定位到桌面），然后使用指令MD5 Control1_Data1.fq.gz生成Control1_Data1.fq.gz文件的MD5值，而使用MD5 *可以生成此目录下所有文本文件的MD5值。

A：分析结果中的基因ID、基因名（symbol）和转录本ID都来源于参考基因组，一般只有部分物种（如小鼠、大鼠等）有类似人中的基因名（如GAPDH、Gapdh），当然即使是大鼠目前也有部分基因没有基因名而以基因ID代替。其他研究较多的物种可能有自己专门的网站和基因命名方式，比如拟南芥以AT1G开头的基因名，番茄以Solyc01g开头的基因名等。另外部分物种可能有多个数据库，比如水稻、茶叶等，不同数据库间的基因ID、基因名和转录本ID形式可能都不同。此外对于研究较少物种或新组装的基因组可能基因ID、转录本ID和基因名三者相同。

A：所有转录本的序列信息已放入3_transcripts_sequence文件夹中，您可以使用基因对应的转录本ID在merged.fa中查询其序列。OmicStudio已上线序列批量提取工具，提交merged.fa文件和转录本ID列表，即可批量提取关注转录本的序列。另外对于Windows平台，可以使用notepad++软件（https://notepad-plus.en.softonic.com/）打开此文件（merged.fa），使用快捷键（Ctrl+F）检索转录本所在位置，将序列复制出来即可；对于Mac平台，可以使用CotEditor打开此文件，使用快捷键（command+F）检索转录本所在位置，将序列复制出来即可。

A：联川生物针对定量和差异分析提供了基因层面和转录本层面两个维度的信息。summary/4_transcript_expression/5_isoforms_fpkm_expression.xlsx为转录本定量结果，summary/6_differential_expression/*/*_Transcript_differential_expression.xlsx（*为差异比较组）为转录本差异分析结果。由于可变剪切，一个基因可能转录多条转录本，一般情况下，来自于同一基因的转录本具有相同的基因名和基因ID，如果我们在转录本分析结果中查询基因名，查到的是基因对应转录本的表达量信息和差异信息，不同转录本的表达量和差异信息往往是不同的。

当然在基因定量结果（summary/4_transcript_expression/4_genes_fpkm_expression.xlsx、summary/6_differential_expression/*/*_Gene_differential_expression.xlsx）中可能检索基因名（极少情况）也能检索到多个结果，原因是不同的基因座具有相同的名称。

A：一般基因是多功能性的，其在不同发育时期、不同组织和不同状态下可能具有不同的功能。基因富集分析是将所有显著性差异基因向Gene Ontology数据库（或其他数据库）的各Term映射，也就是会统计显著性差异基因在数据库中收录的所有功能（或涉及条目、通路）。如果一个基因差异显著，那么其涉及的所有通路都会被富集，只是通路的富集程度不同。因此如果出现显著富集通路与自己研究不相关或相关性不大的情况，可以直接忽略；或者，如果有兴趣的话，也可以探索其他组织中有明确功能报道的基因是否在我们研究的模型中发挥作用以及发挥怎样的作用。

除了上述情况，极少情况我们会遇到富集分析结果中出现不属于本物种的通路，比如动物样本富集分析结果中出现光合作用（Photosynthesis），原因可能是在进行基因的注释时没有区分物种。一般如果参考基因组中没有基因的GO、KEGG注释，同时GO、KEGG官网又没有收录此物种的情况下，通常策略是将基因序列与GO、KEGG收录物种的蛋白序列进行比对，以获得相关注释，然后基于此进行富集分析。如果在比对时没有区分物种或大类，可能会出现动物基因获得了植物的通路注释，从而在富集分析结果中体现出来。如果遇到此类问题，可以重新进行基因注释和富集分析。

A：比较组格式为：<实验组>VS<对照组>，即分析流程判断VS前为实验组，VS后为对照组。实验组和对照组的顺序只是影响差异基因上下调的判定（上调和下调刚好对调），对p值和q值无影响，对差异基因的判定和富集分析结果无影响，只需要将差异倍数（fc）变为倒数(1/fc)，log2fc变为相反数（-log2fc），up和down对调即可。上下调是相对而言的（比较FPKM的大小可以判断在哪组或哪个样本中的表达量高），A相对于B下调也可以描述为B相对于A上调。

A：summary/6_differential_expression/*/*_Gene_differential_expression.xlsx是某比较组的差异分析结果，在表格中的KEGG列使用Excel或WPS筛选”NF-kappa B signaling pathway“即可筛选与NF-kappa B signaling pathway相关的基因，结合“significant”为“yes”的筛选可以获得与NF-kappa B signaling pathway相关的差异显著基因，然后复制基因名（gene_name列）和样本表达量列（FPKM值列）到新的Excel，将此Excel上传至OmicStudio的热图绘制模块，即可进行热图绘制（对于有生物学重复的情况，在“绘图前数据处理”部分选择“按行中心化”和“按行标准化”，数据处理可以选择“不处理”或其他两个选项）。

其他常见问题FAQ

1. 为什么差异分析表格中只有 100 行，99 个基因的信息？

2. 样本间的相关性有何意义，是如何计算的？

3. 差异基因数目是m 个，为什么 GO、KEGG 富集分析中的 TS 数目小于 m?

4. 热图使用什么数据绘制的，从哪个文件中可以找到绘图数据?

4_genes_fpkm_expression.xlsx 是所有样本基因表达量文件，其中基因名列和所有 FPKM 列（表达量，其他组学也是上传表达量矩阵），将其复制到新 Excel 即可得到基因表达量矩阵。从中筛选我们关注的基因即可针对我们关注的基因进行热图绘制。

summary\6_differential_expression\*\*_Gene_differential_expression.xlsx 是*（差异比较组）的差异分析文件，取关注基因的基因名列和所有 FPKM 列（表达量）至新的Excel 中即可获得这些基因的热图绘制文件。热图绘制可以在 OmicStudio（https://www.omicstudio.cn/tool/4）进行绘制。

Zsample-i =[(FPKMsample-i)-Mean(FPKM of all samples)] / [Standard deviation(FPKM of all samples)]

当对数据进行如下处理时等同于 Z 值，建议在绘图时参照如下勾选。其他的颜色、图片比例可以基于个性化选择进行调整。

1） 如果绘制热图时基因数目太多（比如超过 100 个），不建议显示基因名称，因为如果需要显示清楚每个基因的名称，图形可能很长。如果是此情况，可以不显示行名，并通过图片高度和宽度调整比例以体现选择基因的表达量变化趋势；

2） 一般聚类可以将表达模式相似的基因聚在一起，从而更好的通过热图直观看出表达趋势。如果需要输出的热图和输入的基因顺序或样本顺序一致，可以不对行（一般是基因）或列（一般是样本）进行举例。

5. 差异基因GO/KEGG 富集分析散点图怎么看？数据筛选的依据是什么？

6.差异基因GO 富集分析柱状图为什么始终都是 biological_process（25 个），cellular_componen（15个）molecular_function（10 个）？数据筛选的依据是什么？

7. 请问为什么差异表达分析结果中，最后significant 都是no?

8. 请问为什么差异表达分析结果中，有的基因q 值小于 0.05，定义为差异显著基因，有的 q 值小于 0.05，定义为非差异显著

基因?

9. Coverage 怎么计算的？0-1,2-5 数值表示什么意思?

10. 你们云平台上的通用版富集分析结果中的 GO 富集分析柱状图的纵坐标百分比是怎么计算的?

11. KEGG 富集的表格中那个 p 值有什么用？是越小越好？大于一定值这个结果就不能用了吗？

12. RNA-Seq 测序中，增加生物学重复个数和单个样本数据量对差异分析的影响是什么？

13. 链特异性建库是什么，有什么优势？

链特异性建库的关键就在于合成 cDNA 的第二链时，由 dUTP 代替dTTP，然后用 UDG 处理，第二链就会降解，而第一链保留下来, 继而测序。因此，测序得到的转录本序列信息，只是来源于第一链的。

链特异性建库有如下优势：

1） 定量更准确

由于链特异性测序方法可以区分转录本的来源，因此在计算某些转录本的表达量时，可以排除来自其互补链的转录本。

2） 可变剪切事件的检测更准确

因为链特异性文库可以排除反义链上 antisense 转录本的影响，可变剪接事件的检测假阳性更低。

3） Non-coding transcript 的检测

链特异性文库可以显著提高 non-coding transcript 的检出效率。对于 antisense 的 non-coding 转录本，如果用普通文库，是无法区分的；如果是基因间的 non-coding 转录本，普通文库无法确定转录本的方向。

4） 原核生物操纵子（operon）的预测

原核生物的基因是多顺反子的结构，反义转录本上的基因，如果不加区分，那么对应位置的基因表达量会计算不准确，并且预测 operon 以及基因结构也更不准确。

5） 组装结果更真实

一般的转录组组装出来的 unigene 既包括编码转录本，也包括一些非编码转录本（比如lncRNA），但是如果不区分正反链， 那么有互补配对关系的编码与非编码转录本会被组装成一条转录本。

链特异性文库详细情况，请点击：

你可能做了假转录组！揭秘yyds“真”转录组秘密——链特异性文库

 14.GO、Pathway 富集分析中是否一定需要选择显著富集的通路?

1） 基因调控是个系统，不要仅仅看成一个一个孤立的 Pathway

基因调控是个系统，可以从两个层面进行解读：

例如。基因 X 理论上在不同情况下，有可能参与 A、B、C 通路。在某个生物处理下，或许基因 X 只在A 通路里起作用。但如果进行基因注释的话，X 同样也会被注释到 B、C。所以，富集分析的结果总是会涉及特别多的通路。例如，研究人的项目，无论什么研究背景，常常会富集到帕金森综合症通路。不是你的材料真的得了帕金森综合症，只是那些与你实验处理相关的基因，在一定条件下也可以参与到帕金森综合症的过程，所以被注释到了这个通路里。

2） Pathway 富集分析的统计假设，并非在任何情况下都适用

例如，代谢物 X 的合成修饰。基因 A、B、C 分步骤参与合成的 3 个步骤。基因 A 给X 前体加了羟基，然后传递到下游；基因 B 又给 X 前体加了苯环，再传递到下游；基因 C 又给X 的前体加了个乙酰基 ，完成 X 的合成。那么，基因 A、B、C 是参与了的相同的通路。如果基因 A 发生表达量变化，会直接调控影响 B、C 的表达量变化吗？看来很有可能不会，所以从RNA-seq 差异分析的富集分析结果中，这个通路是不显著的。那么基因A 的表达变化是否有生物学意义？当然有，因为代谢物 X 的合成的确受影响了。

类似的例子，理论上 DNA 差异甲基化的结果，就不能看Pathway 富集分析的结果。1 个Pathway 中的 1 个基因的 DNA 甲基化变化，就足以改变这个通路的基因表达，而不需要整个通路的甲基化都发生变化。DNA 甲基化、组蛋白 CHIP-seq 的结果，其实只看功能注释或通路注释就足够了，不需要考虑富集。

所以，我们还是要观察、理解某个核心 Pathway 中基因的相互作用，才能判断其中的基因变化是否有生物学意义，而不

仅仅看富集分析的 p 值或 q 值。

3） 目前的 Pathway 是不完整的

总结：Pathway 和 GO 富集分析结果的解读，应该从生物学意义的角度出发，p 值和 q 值只是个参考而已，那些不显著的通路也值得解读（从功能注释的角度解读，而不是从富集分析的角度解读）。只要结果可以解释，有意义，不用太迷信 p 值。

 15.1. TPM 是什么，标准化方法比 FPKM 更好？

一．RPKM、FPKM、TPM 的影响因素

RPKM、FPKM 和 TPM 都是基于基因的 reads count 标准化得到的，而影响 reads count 除了基因的客观表达水平外（我们关注的），还包括基因的长度和样本的测序深度（简单理解为测序数据量）。

在同一个样本中，基因越长，随机打断得到的片段就越多，该基因被测到的概率就越大，比对到该基因的reads 就越多，从而

reads count 越大；不同样本里，样本的测序深度越高，同一基因被测到的次数越多，比对到该基因的 reads 数就越多，从而 reads

count 越大。

如果我们想比较同一样本不同基因的表达量，只需要标准化基因长度影响即可；如果我们想比较同一基因不同样本中的表达量，只需要标准化测序深度的影响即可；如果我们既想比较同一基因在不同样本中的表达量，又想比较同一样本不同基因的表达量，那么需要同时对基因长度和测序深度进行标准化。

二．RPKM、FPKM、TPM 的区别

早期的转录组测序为单端测序，以 RPKM（Reads Per Kilobase Million）值作为基因的标准化方法，其先将测序深度标准化， 然后将基因长度标准化。而随着双端测序技术（现在的转录组测序一般均为双端测序）的出现与普及，逐渐使用FPKM 替代

RPKM 进行基因的标准化。FPKM 使用 Fragments 概念代替 Reads，是目前 RNA 测序最广泛使用的标准化方法之一。

对于双端测序，paired-reads 表示一个 Fragments 从两端测序产生的 reads，通常以 read1 和 read2 表示。如果一对 paired-reads 都比对上了，那么这一对 pair-reads 称为一个 fragment；如果一个比对上了，另一个没比对上，那么这个比对上的 reads 就称为一个 fragment。

FPKM 与 RPKM 类似，同样是先将测序深度标准化，然后将基因长度标准化。当为单端测序时，FPKM 和 RPKM 是相同的。但是按照 FPKM 值进行标准化的方式有个特点，即每个样本的FPKM 值总和并不一致。基于此有研究者提出了 TPM 的标准化方法，以解决 FPKM 总和不一致的情况。

TPM 表示 Transcripts Per Million，其是先对基因长度标准化，再对测序深度标准化，与FPKM 正好相反。基因的 TPM 值等于基因 FPKM 值在样本中的占比（百分数）乘以 10⁶，每个样本的基因的 TPM 值总和均为10⁶。

三．一定需要使用 TPM 替换FPKM？

那么我们是否一定需要以 TPM 代替FPKM 进行标准化呢？答案是否定的，原因如下：

1） 当两个样品中的绝大多数基因都发生相同方向变化时（概率比较小），TPM 和 FPKM 均不可靠。特别是当超高表达量的某基因 A 表达量发生变化时，在测序结果中会测得基因 A 更多的 reads，其他基因的 reads 就会相对较少，使用 FPKM 值进行标准化时，会得到其他基因的表达量都降低的假象；此时如果使用 TPM 的标准化方法，其标准化表达量会进一步被影响。

2） TPM 的标准化方式影响的是基因的表达量数值，但是对于差异分析的 p 值和 q 值无影响，因为 p 值和 q 值的计算是基于的基因的 reads count，不是基于 FPKM 或TPM。

3） TPM 影响的是差异分析中基因的差异倍数。TPM 是 FPKM 的百分数，虽然各样本的 FPKM 值总和不一致，但 TPM 和 FPKM值是线性相关的，基于 FPKM 和TPM 值对于大部分基因的变化趋势判断是一致的。

4） 研究者基于 FPKM 值的定量结果能够验证出差异并基于此进行下游机制研究（这也是一直沿用的重要原因），而且目前大量文章依然采用 FPKM 作为基因表达量的标准化方式，特别是顶级期刊依然不断有采用FPKM 值方式的文章发表，采用FPKM 或 TPM 不是制约文章发表的因素。

5） 目前尚没有生信研究文章从整体上表明 TPM 在定量上比 FPKM 有什么优势，TPM 只是解决了 FPKM 总和不一致的问题。如果科学界内有广泛共识认为 TPM 比FPKM 更有优势，结果更准确、更可靠，那么提供相关测序和分析服务的公司肯定会快速跟上的，但是目前并没有广泛共识。

6） TPM 和 FPKM 可能对那些变化微弱的基因影响较大，因为两者定量后可能得出的结论有差异，比如上下调（差异倍数和

1 的比较；0.99 和 1.01 的结论是不同的）以及显著性（在满足统计学阈值的情况下，差异倍数与 2 或 0.5 的比较；1.99 和 2.01 的结论是不同的）。但是基于某个基因的验证结果，我们不能得出 TPM 和FPKM 谁更好的结论，因为换个基因您可能会得出相反的结论。比如在p < 0.05 或 q < 0.05 的基础上TPM 计算出的差异倍数为 2.01，FPKM 计算出的差异倍数为 1.99，以默认的差异显著阈值来看，TPM 分析出是差异显著的，FPKM 分析出是差异不显著的。但是 1.99 和 2.01 没有本质差别。又比如 TPM 算出某基因变化倍数为 1.8 倍，FPKM 算出变化倍数为 1.5 倍，不能因为 TPM 算出的差异倍数可能与qPCR 算出的差异倍数更接近，说 TPM 更准确。可能产生差异的根本原因是二代测序和 qPCR 这两个技术本身的差异。

测序是基于比对到基因的 reads 去衡量基因的表达水平，也就是基于基因的所有外显子；而 qPCR 只是基于基因的一两个外显子去衡量基因的表达水平。

7） 不用过度纠结于 FPKM 还是 TPM，应该将精力放在数据的深度挖掘和实验验证上。如果您依然想采用 TPM 值，那么可以基于 FPKM 值进行TPM 值的计算。某样本 S（Sample S）中基因 A（geneA）的TPM 值为：TPMgeneA = （FPKMgeneA / Total FPKMSample S） X10⁶（使用Excel 或 WPS 求取样本S 的 FPKM 值总和并计算TPM）。

附上两篇使用 FPKM 进行标准化的参考文献：

16. 构建基因 A 敲除的细胞系或动物模型，为什么测序结果中还能够检测到基因 A 或基因 A 的表达水平没有明显变化？

如果我们针对的是目标基因的启动子区域（一般针对 RNA 基因），那么理论上，无论是转录水平还是蛋白水平都是下调的。但是如果针对的是外显子，那么比较转录水平是没有意义的，有 WB 结果支持即可。

那么为什么在转录组中检测到敲除基因的表达，但是蛋白却又明显下调呢？原因是虽然移码突变不一定会影响转录，但是转录出来的 RNA 不能翻译（相当于 mRNA 变成了 lncRNA），从而影响了蛋白水平。

如果我们研究的基因敲除模型并不是自己构建的，可以先咨询敲除是针对基因的哪个区域，再去分析可能的原因。

17. GSEA 分析结果和 GO、KEGG 富集分析结果结论不一致，我该相信哪个？

研究背景

非酒精性脂肪肝（NAFL）及其后遗症是日益严重的全球健康问题。非酒精性脂肪肝是非酒精性脂肪肝疾病（NAFLD）的第一阶段，指在没有明确原因的情况下，大于5%的肝细胞内脂肪过度沉积的病理状态。NAFL可进展为非酒精性脂肪肝炎 (NASH)，甚至进一步发展为肝硬化和肝细胞癌 (HCC) 。NAFLD难以诊断、分期和监测，且目前没有被批准的治疗药物。因此，发现并鉴定潜在药物靶点和生物标志物用于NAFLD的诊断，检测及治疗至关重要。

新文速递

安进制药公司旗下deCODE genetics研究团队在Nature Genetics 杂志（IF=38.33）发表题目为“ Multiomics study of nonalcoholic fatty liver disease”的研究论文。研究团队对非酒精性脂肪肝、肝硬化和肝癌进行了全基因组关联分析（GWAS），并将基因组研究结果、RNA转录组学和蛋白质组数据进行了多组学结合分析。

研究解读

来自英国生物样本（UK Biobank）的36,116例样本和冰岛等国家的9,491个临床病例，在该项目中进行了系列研究，其中血浆蛋白组研究部分，研究团队采用Olink Explore 1536 panel及其他方法，分析了47,151名英国生物样本库个体的1,459种血浆蛋白和35,559名冰岛人的4,907种血浆蛋白丰度，以深入了解影响非酒精性脂肪肝或肝硬化发病机制的蛋白质，并寻找相关的生物标志物。研究发现了2741种蛋白质与NAFL（冰岛=234例，UKB=572例）相关，948种与肝硬化（冰岛=111例，UKB=303例）相关。结果表明，蛋白质组学可以区分NAFL和肝硬化，并为开发NAFL的无创评估方法和新的治疗方案提供见解。

为探索在血浆中测量的蛋白质是否可能在疾病中具有因果效应，作者利用每个蛋白质的trans-pQTLs进行孟德尔随机化计算。转铁蛋白受体蛋白1（TFRC）的pQTLs的效应与它们跟PDFF的效应成正比，表明TFRC可能在NAFL中有因果作用。除TFRC外，分析表明，血浆中蛋白质水平的改变不是疾病的原因，而是疾病的结果，因为NAFL的位点对PDFF的效应与它们对蛋白质水平的效应成正比。作者使用PANTHER v.16.0工具进行了路径分析，以研究对变异位点-肝病关联的理解。与NAFL变异位点相关的273个蛋白质富集于多种代谢和分解过程，包括激素、脂类、酒精、维生素、类固醇和单羧酸的代谢以及其他途径和生物过程。

作者进一步比较了被诊断为肝硬化病人和非酒精性脂肪肝病人的血浆蛋白组。在冰岛队列中，非酒精性脂肪肝和肝硬化之间最显著的差异蛋白，首先是钙质蛋白2 （CSTN2），其次是胰岛素样生长因子结合蛋白2（IGFBP2）。与人群对照组相比，IGFBP2的水平在肝硬化患者中升高，但在非酒精性脂肪肝患者中降低。在UKB中，非酒精性脂肪肝和肝硬化之间最显著的差异蛋白首先是血栓软骨素2，其次是IGFBP7。而NAFL与大人群相比，氨基酸酶-1（ACY1）显示出最强的关联。

为研究血浆蛋白是否能有效区分非酒精性脂肪肝和肝硬化，该项目使用肝酶、年龄、性别和BMI作为基线，以及使用血浆蛋白组和遗传风险评分（GRSs）来训练惩罚性逻辑回归模型。结果表明，用血浆蛋白组训练的模型在区分NAFL和肝硬化、NAFL和大人群以及肝硬化和大人群方面的表现优于其他模型。

写在最后

deCODE Genetics研究团队利用多组学数据阐明了非酒精性脂肪肝的遗传机制，对pLOF变异、血液RNA表达和血浆蛋白质组学数据的整合分析找到了致病基因以及其功能变化是否有影响发病机制。同时，进一步证明了非酒精性脂肪肝风险等位基因对其他疾病和性状（血脂，蛋白质等）的不同影响，并表明血浆蛋白质组学有可能对非酒精性脂肪肝进行诊断和分期，成为诊断和监测疾病的新型无创工具。

参考资料：

1. Multiomics study of nonalcoholic fatty liver disease. Nature Genetics volume 54, pages1652–1663 (2022).

    相关阅读