手把手教你从转录组数据中筛选关键基因

转录组分析最重要的目的之一就是筛选候选基因，那怎么从海量数据中筛选目的基因成为大家的困惑。这里小编经过归纳总结得出了下面的三种方法：

筛选关键基因的方法有3类：

1.表达量﹢功能富集

2.表达量﹢实验

3.表达量﹢序列

根据上面的三种方法不难看出，筛选关键基因的核心是表达量，其实转录组的核心就是以表达量为核心展开其他分析的，然后再附加其他一些信息，找出目标基因；最后将分析结果与研究目的巧妙融合，如果再做一些基因功能实验验证（高分文章必备，可以选择不做），那一篇高质量的文章就问世了。下面小编从各方法的详细应用、实际操作、列举实例等方面给您详细说明。

1.表达量﹢功能富集

该方法就是经典的转录组分析流程，帮您筛选关键基因。首先需要做差异表达分析，然后用图形展示差异表达基因数量，具体参见《转录组图形专题之差异基因相关图形介绍》。当然也可以做趋势分析或WGCNA分析，然后拿差异表达基因分析或趋势分析或WGCNA分析的显著基因集做富集分析，富集分析详解及如何挑选显著通路参见《快来领取全套的富集分析相关图形详解吧！》，之后筛选出显著通路，再筛选该通路下的显著基因即可。

如转录组分析与狼尾草叶色相关的关键基因[1]

该文的研究对象狼尾草是一种观赏草植物，在高光环境下会产生紫色的叶子，在低光环境下会产生浅紫色或绿色的叶子，但后者大大降低了其美学吸引力。因此，该文旨在鉴定与叶片着色相关的关键基因，并阐明参与狼尾草叶子的颜色变化的分子机制。

差异表达基因分析总共鉴定了19043个DEGs，与T0（未处理阶段的叶子）阶段的表达相比，在T1（遮阴12天后新叶子完全变绿的阶段）阶段上调和下调的基因分别为10761和8642。KEGG富集分析发现，显著富集的通路主要有类黄酮的生物合成，黄酮和黄酮醇的生物合成以及类胡萝卜素的生物合成。基因筛选发现存在与叶绿素代谢有关的31个差异表达基因，其中21个与叶绿素的生物合成有关，有10个与叶绿素的降解有关，以及3个与叶绿素降解调控有关的转录因子，花青素的合成和积累有31个关键酶基因，4个可能参与花色苷代谢调控的转录因子（图1所示）。

图1 候选基因列表

怎样从表格中快速筛选出感兴趣的基因等信息呢？方法很多，可以写脚本用python去查找关键基因信息，还可以用Excel中的vlookup函数查找我们关心的基因相关信息，如：基因长度、基因在样品中的表达量、差异倍数、注释信息等。

下面用大家熟悉的Excel作以展示。

1.一般我们在富集分析后，就能知道哪些基因被富集到通路下（如图2 所示），第一列为KEGG_A_class的名称，第二列为KEGG_B_class的名称，第三列为Pathway名称，第四列为富集到通路的基因count数，第五列为Pathway ID，第六列为富集到该通路的基因ID，第七列为富集到该通路的K号。

图2 KEGG富集分析结果表

2.然后找到与研究目的相关的显著富集通路的第六列的其中一格，复制粘贴到一个新的excel表中，然后在excel的工具栏点击“数据”→“数据工具”→“分列”→选择“分隔符”后单击“下一步”→在分割符号中选择以“分号”分隔（图3）→单击“完成”。最后全选后转置，将基因ID转换为列（图4）。

图3 以“分号”为分隔符的操作

图4 转置后的结果

3.准备好如下图5的表格，不管是公司分析数据，还是自己做分析，都会有一个汇总表，表里有所有基因ID ，序列长度，NR注释，Swissprot注释，Pfam注释，KEGG注释，FPKM值，差异倍数，FDR值等。

图5 转录组分析结果汇总表范例

4.在excel中用vlookup函数来提取有效信息列。该函数的介绍如下：

Vlookup函数需要输入4个参数：

（1）要查找的值，比如：基因ID；

（2）需要查找的区域；

（3）区域中包含返回值的列号，也就是在查找区域中要匹配的第几列信息；

（4）精确匹配或近似匹配，一般选择精确匹配。其中0/False表示精确匹配，1/Ture表示近似匹配。

5.例如我们想要知道转置后的基因ID长度，在新列表中，插入vlookup函数，按要求把4个参数输入（如图6）。点击确定，第一个基因的长度即被找到。

图6 vlookup函数的4个参数介绍

6.填充所有单元格（图7），同时用同样的方法提取注释信息，差异倍数等。

图7 填充所有单元格

7.最终的提取结果如下表所示。

表1 候选基因信息详细表

2.表达量﹢实验

该方法的“实验”主要指什么呢？其实就是课题组前期的基因功能实验或文献查阅来借助别人的已经验证过的基因功能结果，之后再根据表达量筛选关键基因。

转录组研究的目的就是寻找与实验设计相关的关键基因，一般来说研究某生理现象都先要阅读大量文献来判断该实验的可行性。如转录组分析揭示了雌性鲁波山羊下颌腺在不同发育阶段的潜在免疫功能相关调节基因或途径[2]

该文研究的目的是通过转录组测序来定位差异表达基因（DEG）在三个不同阶段的表达谱，预测在不同发育阶段的下颌下腺的免疫功能。由于人、小鼠、大鼠、牛的下颌腺都检测到了相关抗体，并且研究发现了下颌腺中的类红细胞分化因子、内皮素、肝细胞生长因子（HGF）、转移性生长因子（MGF）、转化生长因子-α（TGF-α）和其他因子。所以在后续的基因筛选中会重点关注相关因子的表达基因。

从上面例子中可以发现一个套路，先查阅模式生物中该方面的研究结果，然后再结合自己的研究项目筛选出基因，并且重点关注模式生物中已知的功能基因有没有较大的差异倍数。

下面小编用一个例子在Excel中作以展示。

转录组研究大豆与共生真菌的共生过程的响应机理。植物与内生菌根真菌AMF的共生大约发生在3.7亿年以前，这种稳定共生体的形成依赖于宿主与菌根真菌之间持续发生信号关联，并在植物宿主的根系皮层中表达大量新的转录本。目前在百脉根、蒺藜苜蓿等模式植物中研究广泛，从没有真菌到植物与真菌形成共生关系会发生很多分子响应过程，如：植物产生类黄酮、独角金内酯等，植物钙信号传导，真菌与植物的营养物质交换等。所以在本研究的关键基因筛选时，我们会重点关注独角金内酯的化学分子合成的相关基因，钙离子信号传导基因，糖转运或脂质转运的相关基因。

这时我们就可以从图5的总表中提取上述所描述的关键基因了。具体操作如下：

1.先在Swissprot注释列中筛选糖转运基因，这时从110843个基因中找到了30个双向糖转运基因；

图8 筛选步骤

2.然后精细筛选，手动看各数据库注释信息等；

图9 筛选结果

3.最终筛选的基因长度在600-1000 bp左右，NR和Swissprot的注释结果为双向糖转运（bidirectional sugar transporter），Pfam注释信息为（PF03083， Sugar efflux transporter for intercellular exchange），KEGG通路注释信息为Transporters（ko02000）及基因的注释信息为SWEET/ solute carrier family 50（K15382），差异表达倍数较大。筛选结果如图4所示。

图10 筛选结果

3.表达量﹢序列

该方法的“序列”指什么呢？主要指碱基序列或氨基酸序列信息，一般包括可变剪切、碱基突变、融合基因、结构特征（也就是结构预测）、RNA编辑等。实质上该方法就是根据特殊基因结构和表达量差异筛选候选基因。

例如根据预测表达基因的蛋白结构来筛选基因。在筛选真菌的效应蛋白（外泌体）时，候选分泌效应蛋白的标准：

a. 氨基端具有一段信号肽；

b. 无跨膜结构域（transmembrane domain ,TMD）（不排除预测的TMD与信号肽重叠的蛋白质）；

c. 无糖基磷脂酰肌醇锚定位点；

d. 没有将蛋白输送至线粒体或其他胞内细胞器的预测定位信号；

e. 氨基酸数量大约为50~300个氨基酸；

f. 富含半胱氨酸并且包含二硫键；

g. 含有重复序列；h. 具有核定位信号；i. 特异的motif或者domain。

通过上面一些蛋白序列标准的预测可得到候选效应蛋白库，随后再根据实验设计样品的差异表达倍数寻找关键的基因即可。知道了基因ID后根据方法1表达量﹢功能富集的vlookup函数提取数据即可。

例如下面的例子：转录组和蛋白组研究印度梨形孢调节植物质外体中的胞外核苷酸水平并影响真菌定植[3]

目的：细胞外腺苷5'-三磷酸（eATP）是必需的信号分子，通过与动植物中膜相关受体蛋白的相互作用介导不同的细胞过程。eATP调节植物的生长、发育以及对生物和非生物胁迫的响应，它在质外体中的积累诱导ROS的产生，细胞质钙的增加介导对入侵微生物的防御反应。植物-真菌相互作用过程中，生物活性核苷酸触发信号和真菌细胞外E5'NT活性的重要性尚不清楚。因此作者利用转录组和蛋白组的手段研究印度梨形孢的E5'NT能否将ATP、ADP和AMP水解为腺苷和磷酸盐，并研究eATP含量和植物对真菌定植的反应。

结果：蛋白组检测出了102种蛋白，经过结构预测得到了48种蛋白是质外体分泌蛋白，再根据转录组表达数据的差异表达分析，最终筛选出关键基因（图11）。其中效应因子SiE5’NT（PIIN_01005）是一种核苷酸水解酶，它可以在质外体空间产生活性氧物质（ROS），并将核苷酸（eATP）水解为腺苷和磷酸盐，从而改变eATP含量和植物对真菌定植的反应，达到定植的目的。