基于Snakemake构建非模式生物转录组分析框架

李岸林 刘长宁

摘要：指出了在生命科学领域，非模式生物具有越来越重要的地位，对非模式生物的研究大幅拓宽了生命科学研究的广度，促进发现更多新的科学问题，弥补模式生物在研究中的不足之处。传统的转录组分析存在流程操作复杂、工作效率低、入门门槛高等问题。基于Snakemake的转录组分析框架，能够运用数据标签批量处理多个数据，无需使用者多次设置参数，大大提高了工作效率。利用Snakemake初步搭建转录组分析框架，在野生型和FUL基因型番茄样品的RNA--Seq转录本数据进行测试。本研究工作为非计算机背景的研究人员进行转录组分析提供了方便，对相关转录组数据分析研究具有一定參考意义。

关键词：分析框架;非模式生物;Snakemake转录组分析

中图分类号：Q75文献标识码：A 文章编号：1674-9944（2019）14-0207-04

1引言

当前模式生物已经广泛地应用于各个研究领域之中，为揭示中心法则以及基因和转录本等方面的一些基本的生物学规律做出了重要贡献。经典模式生物虽然已被广泛应用于各研究领域之中，但其在生物多样性和功能性等方面存在着很多问题。对非模式生物的研究可以拓展生命科学研究的广度，促进发现更多新的科学问题，弥补模式生物在研究中的不足之处。在遗传和进化方面，非模式生物有着重要的研究价值。

在传统的转录组数据分析流程中，一般分为多个步骤进行分析，在每个步骤中会包含一个或者多个分析软件。研究人员需要通过编程脚本来衔接各个步骤，并且在分析的过程中会有大量重复性的数据分析操作，不仅对于研究人员要求有一定计算机技能要求，还会耗费大量的时间和精力在重复性的工作上。当前测序数据量飞速增长，在大数据时代下，需要更加强大的转录组数据分析工具，以Snakemake为基础的分析框架可以大幅提高转录组数据分析效率，且简单易用，为非计算机背景研究者提供了方便。

Snakemake是基于Python3编写的一个流程管理工具，Snakemake将每个步骤简化为一个rule，利用shell命令或者Python代码进行文件输入和输出文件，还可以利用Python语言对软件运行的过程进行管理。利用Snakemake搭建转录组分析可以提高转录组分数据分析效率。本文初步搭建基于Snakemake的转录组分析框架，并利用测试组数据对框架进行测试。

2框架搭建流程

在框架搭建过程中，需要安装许多转录组分析相关的软件，本框架利用Conda自动下载安装需要的软件，避免了手动下载安装软件耗时多、软件版本错误等方面的问题。.Conda还可以配置相关依赖项，完成在本地计算机系统中运行环境的创建、保存、加载和切换等一系列的操作。在框架参数和软件参数、文件路径等设置方面，本框架通过配置文件config.yaml进行统一设置，无需重复对相同的参数进行设置。同时用户也可以通过配置文件对框架进行设置管理，白定义转录组分析模式。在软件选择方面，对转录组分析相关软件进行研究和对比，选择性能较好、应用广泛、稳定性强的主流软件应用于框架中。在框架搭建中将所有软件转化为Snakemake中的rules，然后利用Python语言对rules进行优化和串联并形成Snakefile文件。根据测序数据类型，将单端数据和双端数据进行分别处理，在流程上分为有参转录组分析和无参转录组分析两种类型。共搭建4个分析流程，包括Refer_SE.sk、Refer PE.sk、noRefer_SE.sk和noRefer_PE.sk，分别用于单端有参转录组分析、双端有参转录组分析、单端无参转录组分析和双端无参转录组分析。具体的分析流程如图1所示。构建框搭建过程中所使用的Snakefile文件及其他配置可以在GitHub上获取。

3数据测试与应用

本测的测试数据来自番茄（Lycopersicon esculen-turn MillJ）的转录组测序数据，为有参单端数据测试。番茄在成熟期间，其果实肉质、色素沉着、糖含量、香气和营养品质等方面经历复杂的生理和生化变化。深人探究其果实成熟的机制，有助于我们优化番茄品种的质量。在果安成熟的机制中，MADS基因是果实成熟中心调节的基因，FUL基因是MADS的基因中的一个关键的调节催熟的基因，将探究FUL基因的存在是否会影响不同阶段的番茄基因差异表达。

3.1测试数据

本次测试数据的参考基因组来自于Sol GenomicNetwork网站中的International Tomato Genome Se-quencing Project中的最新的番茄全基因组，其版本号为SL3.0番茄基因组（GCA_000188115.3）。RNA-seq数据下载于National Center for Biotechnology In-formation（NCBI）数据库，测序设备为Illumina HiSeq2000，数据量大约为20G。样本包括野生型和FUL基因型两种番茄样品，同时前期和后期番茄果实样品各3个，共12个样品。样本的具体信息如表1所示。

3.2数据质量控制结果

在数据质量控制部分，本框架利用FastQC软件进行质量查看，Trimmomatic进行序列的修正和处理。通过Snakemake转录组分析框架，成功获得12个测试数据的转录组质量分析结果，如表2所示。表中显示数据的整体质量较好，质量分数中值（Mean Quality Scores）均大于20;不能分辨的碱基N几乎为O%;读段长度均为98到100bp，长度较为一致。

3.3基因的表达量计算

基因表达量分析最直接的方法是计算比对到每个基因上的读段的数量，在转录组测序中通常称其为量化（count）。通过使用基因组比对的方法，本框架的比对分析过程是利用HISAT2软件进行的，比对分析获得每条reads比对到基因组上的位置，以及每个基因比对的reads数量，这些信息通常以GTF格式文件进行储存。HISAT2软件在比对分析过程中，结合GTF文件中提供的转录本剪接信息，在一定程度上提高了读段比对的准确性。通过比对，总共大约有147百万读段比对到了基因组上，对每个基因比对的reads数量进行统计，如表3所示。

3.4基因的聚类分析

聚类分析热图可以根据转录本的表达量高低，对转录本和样本进行聚类分组，框架利用DEseq2进行聚类分析，R包来完成结果的可视化如图2所示。在聚类分析热图中，通过颜色的深浅变化代表每个基因或转录本在不同样品中的表达情况，红色越深表达量越高，蓝色反之。通过颜色的深浅就可以直观的看到不同样品或分组之间表达谱的差异。此外，将聚类（Clustering）的结果表示在热图的上方，以进化树的结构展示。若基因或样品的表达趋势越接近，则在聚类结果中也更为接近，通常认为能聚到一起的基因簇或样品簇有着更为接近的生物学特征，以此可以推断未知基因的功能。

3.5基因的差异性表达

差异表达分析的目的是分析在不同的环境或者不同组织中的基因表达差异。通过对基因的表达量来发现是否存在基因的差异性表达。衡量基因差异性表达首先需要将基因的表达量进行标准化，接着利用统计模型进行差异性表达检验。差异性表达分析主要指标是Fold-change，Fold-change指的是两个样本之间表达量差异的倍数。p值（p-value）是反映了结果的可信度，当p值越小时，基因在不同条件下的差异性表达越可信。而adjusted P或q值（q-value）或是False dis-covery rate（FDR）也叫假阳性率（False Positive Rate）修正下的p值，通常使用q值或FDR来衡量表达差异的可信度。

通过框架DEseq2软件计算得到Fold-change以及p值或q值后，R包绘制得到火山图（Volcano Plot），如图3所示。在生成的火山图中，横轴logFC代表Foldchange的log值，纵轴代表adjusted p值的loglO的值，通过对横纵坐标阈值的限定，可以非常直观筛选出在两样本间发生差异表达的基因。本框架配置文件中设置p值默认值为0.05。而Fold Change值只有大于2时才认为是显著差异表达，所以横坐标值logFC应大于1。通过这样的筛选条件，可以筛选出显著差异表达的相关基因作为研究目标。

4总结与展望

基于Snakemake搭建的转录组框架降低了流程对于数据格式转化、脚本编写的要求，可以简便、高效、方便地处理转录本数据。本研究初步搭建了非模式生物分析的框架，利用本框架可以高效完成一整套转录组数据分析的流程，获得相应的数据和可视化结果。基于此框架可以使转录组分析更容易进行操作，为非计算机背景的研究人员进行转录组分析提供了方便。基于Snakemake进行转录本数据分析框架的搭建，能够为后续相关分析框架的构建提供参考。

在Snakemake框架的搭建过程中，提高相关参数的普适性是一个难题，特别是对于没有参考基因序列的非模式生物，无法确定固定的参数，需要用户不断地进行尝试。在框架的功能方面，现在只有转录组分析的基本功能，框架的功能的丰富性有所欠缺，在后期需要进一步丰富框架的功能。另一方面，需要增强框架应用的易操作性，尽量优化用户操作步骤。在框架的软件选择方面，在后期会添加备选软件供用戶选择。在后续工作中，会不断完善框架内容和参数调配，并且针对不同的非模式生物匹配相应的配套参数。对于后续更新的一些框架内容，会定期上传到Github上以供下载。