摘要
为使科研工作者简单高效地分析RNA-seq数据，本研究基于Snakemake工作流程管理系统和Conda环境管理器构建了一个自动化和模块化的工作流程：RNApipe（Github：https：//github.com/ywu019/RNApipe.git），其可对来自任何有参物种的RNA-seq数据自动执行质控、比对、定量、鉴定差异基因，以及GO、KEGG、GSEA等功能注释分析；其中，每一步骤的分析结果均以高质量的可视化图片或报告展示，并保留重要的输出文件。使用RNApipe在多个模式物种中的测试与评估结果表明：RNApipe可以平稳运行，且注释结果准确。与现有的自动化分析流程相比，RNApipe的主要特点包括：工作流程较为完整、默认工具消耗时间与资源较少、适用于任何有参物种、全面的可视化、以及用户友好性（易安装、易使用、易扩展）。研究表明，RNApipe便于研究人员快速地从大型RNA-seq测序数据中获取基本信息。
转录组测序（RNA-seq）是一项对特定组织或细胞中所有的RNA进行高通量测序的技术。近十年来，随着测序技术的发展、测序成本的降低以及分析工具的快速开发，RNA-seq已广泛应用于研究基因表达调控、可变剪切甚至RNA结构等方向，成为生命科学基础研究领域中非常重要的技
目前，RNA-seq的各个分析步骤都已有许多成熟的计算方法。然而，这种分段分析的方法通常会涉及到不同的操作系统和编程语言，这对于科研人员的编程能力提出较高的要求。之前已有多项研究将各个分析步骤集成到自动化分析工作流程中，从而降低对用户编程能力的要求并有效地提高其工作效率。然而，目前已开发的RNA-seq自动化分析工作流程中存在着一些不足之处。例如，VIPE
本研究通过Conda部署环境以确保整个工作流程中的软件快速、平稳地安装，并通过Snakemake工作流程管理系统集成各个分析步骤，最终构建了1个流程更完整且更易使用的RNA-seq自动化分析工具——RNApipe，旨在为研究人员快速、简便地从大型RNA-seq测序数据中获取基本信息提供技术支持。
对于编程经验有限的用户，安装并维护RNA-seq数据分析中必需的工具及其依赖项通常具有挑战性。Conda作为软件包和环境管理系统，可以自动安装、运行和更新软件包及其依赖项，并且支持虚拟环境的创建、切换与移植， 是确保软件管理可持续和可重现的理想工
为减少数据分析过程中手动执行的步骤，传统的工作管道通常使用自定义的脚本将多个工具进行链接，然而这种管道通常与本地计算基础架构高度相关，难以共享与维护，产生的结果也较难重
图1 RNApipe数据分析流程图
Fig.1 The workflow chart of RNApipe analysis
整个RNApipe工作流程（
RNApipe的源代码已经上传至Github供用户下载和使用：https：//github.com/ywu019/RNApipe.git。以下简单介绍RNApipe的安装和使用方法，具体使用方法可见Github中的说明文档（Documentation.pdf）。
1）下载RNApipe工具与安装环境。首先通过Python3.7在Linux系统中安装miniConda3环境，随后通过以下3条命令完成RNApipe工具的安装与环境的创建：①从GitHub下载RNApipe：it clonehttps：//github.com/ywu019/RNApipe.git；②创建环境： conda env create -n RNApipe -f envs/envs.yaml；③激活环境： conda activate RNApipe。
2）用户指定输入文件。在运行RNApipe前需要指定：基因组文件（.fa）、注释文件（.gtf）与压缩的测序文件（.fastq.gz）（
随后修改configs/目录下的样本表（samples.tsv）和配置文件（config.yaml）。用户可以根据实验条件指定样本表中的条件列“condition”与重复列“replicate”的信息，并修改配置文件中的基本设置，配置文件主要包含以下变量：
PROJECT： #项目名称；
READPATH：#测序文件fastq.gz的路径；
END：#测序数据类型是SE（single-end）还是PE（paired-end）；
OUTPUTPATH： #输出文件的路径；
GENOME：#基因组文件路径；
ANNOTATION：#注释文件路径；
CONTROL： #对照组名称，与samples.tsv中condition一致；
TREAT：#处理组名称；
SPEICES_ABBREVIATION：物种简称，在speices_abbreviations中查询；
TRIM： fastp # ［fastp， cutadapt］ 选择工具执行数据预处理；
ADAPTER： AGATCGGAAGAGC # 接头序列；
ALIGN： hisat2 # ［hisat2， STAR］ #比对软件；
BIGWIG： ［"no"］ # ［yes， no］ #是否将bam格式的文件转化为bigwig格式；
DEA： DESeq2 # ［DESeq2， edgeR， limma］ #差异基因分析软件；
QUANTIFY： featureCounts #［featureCounts， htseq］ #定量软件；
THREAD： "10" # 线程数目；
TIME： "5" # 程序延迟时间。
3）运行工作流程。准备好所有的输入文件后，用户只需要通过一条命令（python main.py）即可在本地或集群环境中运行RNApipe整个工作流程。运行结果保存在用户指定的输出目录中，logs/目录下生成各个分析步骤的日志文件，其中记录了程序的运行时间、标准输出和标准错误等信息。
为了让用户了解RNApipe的工作原理和基本功能，本研究利用RNApipe对来自4日龄野生型拟南芥根部（root）和芽部（shoot）的单端RNA-seq完整数据集（PRJNA321304；GEO：GSE81332）进行自动化分析，其中根部设为对照组，芽部设为试验组。
图2 RNApipe的文件与目录结构
Fig.2 The file and directory structures of RNApipe
用户在运行项目时应在project/目录中设置类似于example_data/的项目结构，并将测序数据、参考基因组和注释文件分别置于相应的子目录中（
为了保证分析结果的准确性，RNApipe分别对质控后的数据和比对文件进行质量检测，并将每个样本文件生成的质控报告（html）进行汇总（
图3 数据质量检测报告
Fig.3 Report on data quality inspection
A：读段中碱基质量的平均值 The mean quality value across each base position in the read； B：样本比对情况 A brief mapping summary of all six samples.
在完成质量控制和比对后，RNApipe将表达矩阵（count）进行TPM标准化，以消除测序深度和基因长度对表达量的影响。设置阈值（TPM<1）过滤低表达的基因，得到的表达量矩阵用于可视化图表展示（
图4 样本间基因表达量的可视化
Fig.4 Visualization of gene expression levels among samples
A：热图展示样本间的相关性The heatmap shows the correlation among samples；B：PCA主成分分析PCA principal component analysis；C：箱线图展示样本总体表达量The boxplot shows the overall expression level of samples；D：累积分布曲线图展示样本总体表达量变化The cumulative distribution curve shows the changes of expression levels in samples；E：热图根据样本内所有表达基因的相关性进行聚类The heatmap is clustered according to the correlation of all expressed genes in the samples.
RNApipe使用cor函数与ggplot2包（v3.3.3）将所有样本基因表达谱的成对相关性可视化为热图，相关系数的绝对值越接近1，表明样本之间的相关性越高，这有助于用户进一步识别相关性较低的离群值样本。结果显示拟南芥根部或芽部中各个生物学重复间的相关性较高，而样本间的相关性较差（
RNA-seq下游分析的第一步是鉴定样本间差异表达的基因。RNApipe集成了当下主流的差异表达工具：DESeq2、edgeR、limma。这些工具具有不同的模型和优势（见Documentation.pdf），用户可以根据实验需求在config.yaml中选择其中任一工具进行分析。RNApipe默认使用DESeq2软件在基因水平上执行差异基因的鉴定，DESeq2使用基因表达的原始丰度而不是标准化丰度进行统计检
图5 样本间差异表达基因的可视化
Fig.5 Visualization of differentially expressed genes among samples
A：火山图展示根部和芽部间差异表达的基因 Volcano plot shows differentially expressed genes between root and shoot； B：热图展示样本间变化最大的前30个差异基因 The heatmap shows the top 30 differentially expressed genes.
差异基因的功能注释对于了解特定条件所影响的生物学过程至关重要。基因本体论（GO）分
图6 样本间基因功能注释的可视化
Fig.6 Visualization of gene functional annotations among samples
A：样本间差异基因的GO功能注释 GO functional annotation of differentially expressed genes that among samples；B：样本间差异基因的KEGG通路分析KEGG pathway analysis of differentially expressed genes among samples；C：shoot R-ATH-556833的GSEA富集分析GSEA enrichment analysis of shoot R-ATH-556833；.D：root ATH-00860的GSEA富集分析GSEA enrichment analysis of shoot root ATH-00860.
GO和KEGG富集分析往往侧重于识别两组间差异表达显著的基因，这可能导致遗漏部分差异表达不显著却有重要功能和意义的基因，另外，GO和KEGG分析无法判断差异基因的变化方向是上调还是下调。基因集富集分析（GSEA）不需要指定明确的差异阈值，算法会根据实际数据的整体趋势对基因的表达差异进行排序，随后检验数据库中预先设定的基因集富集在该排序列表的顶端或底端，因此，GSEA检测到的是基因集而不是单个基因的表达变化。GSEA图展示了脂质和脂蛋白通路中的基因集的表达在拟南芥的芽中呈现上调的趋势，而卟啉与叶绿素代谢通路中的基因集的表达在根中呈现下调的趋势（
与现有的自动化分析工具对比（
工作流程 Workflow | 质量控制 Quality control | 物种 Organism | GO/KEGG分析 GO/KEGG analysis | GSEA分析 GSEA analysis | 安装方式 Installation | 硬件需求 Hardware requirement | 编程需求 Programming requirement | 文献 References |
|---|---|---|---|---|---|---|---|---|
| RNApipe | √ | All | √ | √ | 简单Easy | 低Low | 低Low | 本研究This study |
| RASflow | √ | All | × | × | 简单Easy | 低Low | 低Low |
[ |
| UTAP | √ | 5 | × | × | 简单Easy | 高High | 低Low |
[ |
| ARMOR | √ | All | × | × | 简单Easy | 高High | 低Low |
[ |
| VIPER | √ | 2 | √ | √ | 简单Easy | 高High | 低Low |
[ |
| BioJupies | × | 2 | √ | × | 网站Web | 低Low | 低Low |
[ |
| hppRNA | √ | 2 | × | × | 中等Medium | 低Low | 中等Medium |
[ |
| aRNApipe | √ | All | × | × | 困难Hard | 高High | 高High |
[ |
| RNACocktail | × | All | × | × | 困难Hard | 低Low | 高High |
[ |
注Note：√：支持 Support；×：不支持 Unsupported；All：所有的有参物种All participating species.
利用RNA-seq数据鉴定差异表达基因为研究重要的科学问题提供了思路，若将整个过程的分析步骤整合为一个自动化、模块化的流程，将极大地方便科研工作者的使用。为此，本研究开发了一个适用于任何有参物种的、流程更加完整的RNA-seq自动化分析工作流程：RNApipe。RNApipe工作流程旨在分析来自任何有参物种的RNA-seq数据。目前，RNApipe已在以下5个模式物种的完整数据集中完成测试与评估：人（PRJNA390636；GEO：GSE100075）、小鼠（PRJNA630257； GEO：GSE149838）、拟南芥（PRJNA321304；GEO：GSE81332）、粗糙脉孢菌（PRJNA392079；GEO：GSE100539）、酵母（PRJNA318684；GEO：GSE80357）。分析结果表明，RNApipe运行稳定，功能注释结果符合已发表文献的数据。实际数据集的可视化结果存放于https：//github.com/ywu019/RNApipe_pj.git.，与现有的自动化分析工具对比，RNApipe具有安装容易、使用方便、易于扩展、功能丰富等特点。RNApipe将帮助科研工作者以一种简单、高效且可重复的方式从RNA-seq数据中获取有价值的信息。
RNApipe目前主要适用于对二代Illumina测序平台产生的短读长（short-read）数据进行差异基因分析。在RNApipe的基础上，用户可以通过修改rules命令来扩展其他功能，例如可变剪切与单核苷酸变异的检测等。另外，由于常规RNA-seq是对所有细胞的RNA进行测序的技术，因此无法辨别细胞的类型和空间信息。近年来兴起的单细胞转录组测序技术（scRNA-seq）常应用于研究不同组织或器官中（例如正常和癌症）差异基因的表达，为研究单细胞水平的基因表达谱提供了机会。空间转录组测序技术通过将成像技术与scRNA-seq相结合，可以绘制出转录本在组织中表达的位置，进一步提高组织分辨率。scRNA-seq技术与空间转录组学技术的普及同时伴随着研究人员在数据分析方面面临的挑战，RNApipe的结构设计则可以为scRNA-seq自动化数据分析工具的构建提供参考。
参考文献References
STARK R，GRZELAK M，HADFIELD J.RNA sequencing：the teenage years［J］.Nature reviews：genetics，2019，20（11）：631-656. [百度学术]
NAGALAKSHMI U，WANG Z，WAERN K，et al.The transcriptional landscape of the yeast genome defined by RNA sequencing［J］.Science，2008，320（5881）：1344-1349. [百度学术]
CONESA A，MADRIGAL P，TARAZONA S，et al.A survey of best practices for RNA-seq data analysis［J/OL］.Genome biology，2016，17：13［2022-01-13］.https：//doi.org/10.1186/s13059-016-0881-8. [百度学术]
CORNWELL M，VANGALA M，TAING L，et al.VIPER：visualization pipeline for RNA-seq，a Snakemake workflow for efficient and complete RNA-seq analysis［J/OL］.BMC bioinformatics，2018，19（1）：135［2022-01-13］.https：//doi.org/10.1186/s12859-018-2139-9. [百度学术]
WANG D P.hppRNA—a Snakemake-based handy parameter-free pipeline for RNA-Seq analysis of numerous samples［J］.Briefings in bioinformatics，2017，19（4）：622-626. [百度学术]
MONIER B，MCDERMAID A，WANG C，et al.IRIS-EDA：an integrated RNA-Seq interpretation system for gene expression data analysis［J/OL］.PLoS computational biology，2019，15（2）：e1006792［2022-01-13］.https：//doi.org/10.1371/journal.pcbi.1006792. [百度学术]
MARINI F，LINKE J，BINDER H.Ideal：an R/Bioconductor package for interactive differential expression analysis［J/OL］.BMC bioinformatics，2020，21（1）：565［2022-01-13］.https：//doi.org/10.1186/s12859-020-03819-5. [百度学术]
TORRE D，LACHMANN A，MA’AYAN A.BioJupies：automated generation of interactive notebooks for RNA-seq data analysis in the cloud［J］.Cell systems，2018，7（5）：556-561. [百度学术]
SAHRAEIAN S M E，MOHIYUDDIN M，SEBRA R，et al.Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis［J/OL］.Nature communications，2017，8：59［2022-01-13］.https：//doi.org/10.1038/s41467-017-00050-4. [百度学术]
ORJUELA S，HUANG R Z，HEMBACH K M，et al.ARMOR：an automated reproducible MOdular workflow for preprocessing and differential analysis of RNA-seq data［J］.G3 Genes|Genomes|Genetics，2019，9（7）：2089-2096. [百度学术]
ZHANG X，JONASSEN I.RASflow：an RNA-seq analysis workflow with snakemake［J/OL］.BMC bioinformatics，2020，21（1）：110［2022-01-13］.https：//doi.org/10.1186/s12859-020-3433-x. [百度学术]
GRÜNING B，DALE R，SJÖDIN A，et al.Bioconda：sustainable and comprehensive software distribution for the life sciences［J］.Nature methods，2018，15（7）：475-476. [百度学术]
WRATTEN L，WILM A，GÖKE J.Reproducible，scalable，and shareable analysis pipelines with bioinformatics workflow managers［J］.Nature methods，2021，18（10）：1161-1168. [百度学术]
KÖSTER J，RAHMANN S.Snakemake：a scalable bioinformatics workflow engine［J］.Bioinformatics（Oxford，England），2012，28（19）：2520-2522. [百度学术]
CHEN S，ZHOU Y，CHEN Y，et al.Fastp：an ultra-fast all-in-one FASTQ preprocessor［J/OL］.Biochemistry and biophysics reports，2018，34（17）：i884-i890［2022-01-13］.https：//doi.org/10.1093/bioinformatics/bty560. [百度学术]
MARTIN M.Cutadapt removes adapter sequences from high-throughput sequencing reads［J］.EMBnet journal，2011，17（1）：10-12. [百度学术]
BROWN J，PIRRUNG M，MCCUE L A.FQC Dashboard：integrates FastQC results into a web-based，interactive，and extensible FASTQ quality control tool［J］.Bioinformatics （Oxford，England），2017，33（19）：3137-3139. [百度学术]
EWELS P，MAGNUSSON M，LUNDIN S，et al.MultiQC：summarize analysis results for multiple tools and samples in a single report［J］.Bioinformatics，2016，32（19）：3047-3048. [百度学术]
KIM D，LANGMEAD B，SALZBERG S L.HISAT：a fast spliced aligner with low memory requirements［J］.Nature methods，2015，12（4）：357-360. [百度学术]
DOBIN A，DAVIS C A，SCHLESINGER F，et al.STAR：ultrafast universal RNA-seq aligner［J］.Bioinformatics，2012，29（1）：15-21. [百度学术]
OKONECHNIKOV K，CONESA A，GARCÍA-ALCALDE F.Qualimap 2：advanced multi-sample quality control for high-throughput sequencing data［J］.Bioinformatics，2015，32（2）：292-294. [百度学术]
LIAO Y，SMYTH G K，SHI W.featureCounts：an efficient general purpose program for assigning sequence reads to genomic features［J］.Bioinformatics，2013，30（7）：923-930. [百度学术]
ANDERS S，PYL P T，HUBER W.HTSeq：a Python framework to work with high-throughput sequencing data［J］.Bioinformatics，2014，31（2）：166-169. [百度学术]
LOVE M I，HUBER W，ANDERS S.Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2［J/OL］.Genome biology，2014，15（12）：550［2022-01-13］.https：//doi.org/10.1186/s13059-014-0550-8. [百度学术]
ROBINSON M D，MCCARTHY D J，SMYTH G K.edgeR：a Bioconductor package for differential expression analysis of digital gene expression data［J］.Bioinformatics （Oxford，England），2010，26（1）：139-140. [百度学术]
RITCHIE M E，PHIPSON B，WU D，et al.Limma powers differential expression analyses for RNA-sequencing and microarray studies［J/OL］.Nucleic acids research，2015，43（7）：e47［2022-01-13］.https：//doi.org/10.1093/nar/gkv007. [百度学术]
BU D C，LUO H T，HUO P P，et al.KOBAS-i：intelligent prioritization and exploratory visualization of biological functions for gene enrichment analysis［J/OL］.Nucleic acids research，2021，49（W1）：W317-W325［2022-01-13］.https：//doi.org/10.1093/nar/gkab447. [百度学术]
SONESON C，DELORENZI M.A comparison of methods for differential expression analysis of RNA-seq data［J/OL］.eLife，2013，14：91［2022-01-13］.https：//doi.org/10.1186/1471-2105-14-91. [百度学术]
ASHBURNER M，BALL C A，BLAKE J A，et al.Gene Ontology：tool for the unification of biology［J］.Nature genetics，2000，25（1）：25-29. [百度学术]
OGATA H，GOTO S，SATO K，et al.KEGG：Kyoto encyclopedia of genes and genomes［J］.Nucleic acids research，1999，27（1）：29-34. [百度学术]
KOHEN R，BARLEV J，HORNUNG G，et al.UTAP：user-friendly transcriptome analysis pipeline［J］.BMC bioinformatics，2019，20（1）：154［2022-01-13］.https：//doi.org/10.1186/s12859-019-2728-2. [百度学术]
ALONSO A，LASSEIGNE B N，WILLIAMS K，et al.aRNApipe：a balanced，efficient and distributed pipeline for processing RNA-seq data in high-performance computing environments［J］.Bioinformatics，2017，33（11）：1727-1729. [百度学术]