SEQprocess: a modularized and customizable pipeline framework for NGS processing in R package

Backgrounds Next-Generation Sequencing (NGS) is now widely used in biomedical research for various applications. Processing of NGS data requires multiple programs and customization of the processing pipelines according to the data platforms. However, rapid progress of the NGS applications and processing methods urgently require prompt update of the pipelines. Recent clinical applications of NGS technology such as cell-free DNA, cancer panel, or exosomal RNA sequencing data also require appropriate customization of the processing pipelines. Here, we developed SEQprocess, a highly extendable framework that can provide standard as well as customized pipelines for NGS data processing. Results SEQprocess was implemented in an R package with fully modularized steps for data processing that can be easily customized. Currently, six pre-customized pipelines are provided that can be easily executed by non-experts such as biomedical scientists, including the National Cancer Institute’s (NCI) Genomic Data Commons (GDC) pipelines as well as the popularly used pipelines for variant calling (e.g., GATK) and estimation of allele frequency, RNA abundance (e.g., TopHat2/Cufflink), or DNA copy numbers (e.g., Sequenza). In addition, optimized pipelines for the clinical sequencing from cell-free DNA or miR-Seq are also provided. The processed data were transformed into R package-compatible data type ‘ExpressionSet’ or ‘SummarizedExperiment’, which could facilitate subsequent data analysis within R environment. Finally, an automated report summarizing the processing steps are also provided to ensure reproducibility of the NGS data analysis. Conclusion SEQprocess provides a highly extendable and R compatible framework that can manage customized and reproducible pipelines for handling multiple legacy NGS processing tools.


Background
Next-Generation Sequencing (NGS) technology is now widely used in biomedical research fields, and is extensively being used in the clinic [9]. Applications with NGS technology include identification of DNA or RNA sequence variants, and the quantitation of RNA abundances or DNA copy numbers. However, processing and analysis of NGS data remain difficult as data are generally processed through by multiple processing steps, and each step requires different legacy programs. To handle these complex processing steps, several pipeline programs have been released. For example, 'NGS-pipe' [18] and 'NEAT' [17] provide automated pipelines for NGS data analysis. Another tool 'systemPi-peR' provides an NGS analysis workflow in R program that can be customized according to the various NGS applications such as whole-exome sequencing (WES), whole-genome sequencing (WGS) and transcriptome sequencing (RNA-seq) data [2]. However, these tools do not handle the recently updated NCI Genomic Data Commons (GDC) pipelines, which have been used as standard pipelines to process The Cancer Genome Atlas (TCGA, https://cancergenome.nih.gov) data. Moreover, recent progress in clinical applications of the NGS data has generated new platform data, such as cell free DNAs, exosomes, and cancer panels. These applications require customized analysis for data quality control and processing.
With this concern, we developed a SEQprocess that provides fully customizable NGS processing pipelines covering the GDC pipelines and new data for clinical applications. SEQprocess is implemented in an R program, providing six pre-customized pipelines that are widely used as standards in NGS data processing and can be executed easily by non-experts such as biomedical scientists.

Implementation
SEQprocess is a framework implemented in R package, providing pipelines for NGS data processing operated by multiple programs. It can be run from start-to-end with a single command in the R console, or through stepwise customization with an interactive mode. The pipelines are designed to support processing pipelines for DNA and RNA sequencing data, including the data processing steps for quality control of raw sequencing data, trimming, alignment, variant calling, annotation, DNA copy number estimation and RNA quantitation. Each pipeline is modularized to run sequentially or separately. The following programs are supported by the pipelines. Quality control of raw data is assessed by FastQC (https:// www.bioinformatics.babraham.ac.uk). Sequence trimming is performed by TrimGalore (https://github.com/ broadinstitute/picard) or Cutadapt [14]. Sequence alignment is supported by BWA [12], STAR [3], TopHat2 [7], bowtie2 [10], or samtools [13]. Removal of duplicates is performed by Picard (https://github.com/broadinstitute/ picard) and re-alignment by GATK [15]. Variants calling is supported by GATK, VarScan2 [8], MuSE [4], or SomaticSniper [11]. Variant annotation is supported by VEP [16] or ANNOVAR [20]. For RNA-seq data, SEQprocess performs RNA quantitation by HTSeq [1] or Cufflinks [19], and DNA copy number estimation is conducted by Sequenza [5]. These programs are implemented as modularized functions with optimized default parameters. These external programs can be installed easily using Conda package manager (https://conda.io/ en/latest). Subsequent steps for NGS data processing can be easily included or excluded in the pipeline. This modular framework provides a highly flexible and extendable platform; thus, new pipelines for upcoming data types such as single cell RNA-Seq data can be implemented.

Results
The current version of SEQprocess provided six different pre-customized standard pipelines, including the pipelines for GDC processing and the newly adapted clinical applications for cell-free DNAs (cfDNA) or exosomal miRNAs (Fig. 1). These pipelines ran by a one-step command that could be executed easily by non-expert users. For WGS/ WES, a GDC compatible pipeline of TrimGalore-BWA--Picard-VarScan2-VEP was implemented. We also implemented a popularly used standard Custom pipeline of TrimGalore-BWA-Picard-GATK-ANNOVAR.
In addition, SEQprocess could estimate allele frequencies for each variant by calculating the sequence read depths of the mutated and wild-type sequences with a GATK function 'DepthOfCoverage'. For liquid-biopsied cfDNA or targeted sequencing data, such as a cancer panel, an optimized pipeline excluding the duplicate removal step was provided, because cfDNA sequence reads usually have the same sequences. For barcoded data (BarSeq), the duplicate removal step was performed using the barcodes. For RNA-Seq data, a GDC pipeline (STAR-Samtools-HT-Seq) was implemented. A popularly used standard pipeline Tuxedo (i.e., Tophat2-Cufflinks) was also implemented. For miR-Seq data from exosomes, cells, or tissues, the Cutadapt-BWA/bowtie2-HTSeq pipeline was implemented with optimized parameters. SEQprocess operates multiple legacy programs and reference data, which might require installation in the system. Configuration of the installed programs and data could be managed simply by editing the 'data/config.R' file ( Table 1). The current version of SEQprocess supported the Linux-operating system because some of the required programs only support the Linux-operating system. Parallel computation on multi-core machines was also supported by using the 'parallel' R package. In addition, multi-threading support in each program of GATK, TopHat2, BWA, STAR, and Cufflinks could be controlled by the program arguments.
Each step of these pipelines are modularized as a wrapper function in R package to provide an easy customization platform.
Step-by-step pipelines could be conducted by a single command 'SEQprocess' , and which  could be readily customized by setting the function parameters ( Table 2). The processed data were transformed into an R/Bioconductor compatible data type (i.e. 'ExpressionSet'), which is popularly used for the subsequent NGS data analysis for biological interpretation [6]-. Each data object for RNA expression, variant, and DNA copy number was provided with the filename extensions of '.eSet' , '.vSet' , or 'cSet' , respectively. These ExpressionSet data types could be transformed into another data type 'SummariazedExperiment' , i.e. a modified data type of 'ExpressionSet' containing 'GenomicRanges' data type (Fig. 2). These will serve as a framework facilitating the subsequent analyses in the R environment. In addition, SEQprocess provided a report summarizing the processing steps and visualized tables and plots for the processed results (Fig. 3). The report file is automatically generated workflow records for data processing steps, arguments, and outcome results. Moreover, users can find error and processing messages from the log file in each program. These reporting systems will ensure the reproducibility of the data analysis. We have also provided an example data ('inst/example') and a script ('example/example.R').

Conclusions
In summary, SEQprocess provides a highly extendable and R-compatible framework that can be manage customized and reproducible pipelines for handling multiple legacy NGS processing tools.

Availability and requirements
Project name: SEQprocess.