SeqAcademy: an educational pipeline for RNA-Seq and ChIP-Seq analysis

Quantification of gene expression and characterization of gene transcript structures are central problems in molecular biology. RNA sequencing (RNA-Seq) and chromatin immunoprecipitation sequencing (ChIP-Seq) are important methods, but can be cumbersome and difficult for beginners to learn. To teach interested students and scientists how to analyze RNA-Seq and ChIP-Seq data, we present a start-to-finish tutorial for analyzing RNA-Seq and ChIP-Seq data: SeqAcademy ( source code: https://github.com/NCBI-Hackathons/seqacademy, webpage: http://www.seqacademy.org/). This user-friendly pipeline, fully written in markdown language, emphasizes the use of publicly available RNA-Seq and ChIP-Seq data and strings together popular tools that bridge that gap between raw sequencing reads and biological insight. We demonstrate practical and conceptual considerations for various RNA-Seq and ChIP-Seq analysis steps with a biological use case - a previously published yeast experiment. This work complements existing sophisticated RNA-Seq and ChIP-Seq pipelines designed for advanced users by gently introducing the critical components of RNA-Seq and ChIP-Seq analysis to the novice bioinformatician. In conclusion, this well-documented pipeline will introduce state-of-the-art RNA-Seq and ChIP-Seq analysis tools to beginning bioinformaticians and help facilitate the analysis of the burgeoning amounts of public RNA-Seq and ChIP-Seq data.

The expanding importance of RNA-seq and ChIP-seq data is reflected by its explosive growth in terabytes in the primary public repository for storing this data -the Sequence Read Archive (SRA) (Wheeler et al., 2008). This incredible increase in the amount of public data has not been met with an equal increase in the number of scientists who can skillfully and thoughtfully analyze this important resource. Given the fundamental role that RNA-seq and ChIP-seq data, among other next-generation sequencing data types, are likely to play in the coming decades, there is a critical need to teach RNA-seq and ChIPseq analysis to life scientists with diverse interests and backgrounds.
The goal of analyzing RNA-seq data is often to identify and characterize quantitative differences in gene expression between biological samples from two or more groups. For ChIP-Seq, the goal is to characterize DNA-protein interactions. Biological samples may originate from several different study designs including: different tissue types from the same individual (e.g. cancerous tissue vs. non-cancerous tissue), the same strain of cells under different environmental conditions, or the same tissue under a time-course experiment.
There are major barriers to the novice bioinformatician who is interested in learning how to analyze RNA-Seq and ChIP-Seq data. RNA-Seq and ChIP-Seq data are costly to generate (>$1,000/sample) and cumbersome to store; with data from a single sample often occupying several gigabytes of storage space. However, recent advances, including a greater impetus to deposit sequencing data in SRA ("Principles and Guidelines for Reporting Preclinical Research," 2015) and the innovative alignment of streamed sequencing data (Kim et al., 2015), offer new opportunities to overcome these long-standing problems.
The second barrier to entry is inherent to RNA-Seq and ChIP-Seq data. These datasets are large and complex: there are over 20,000 known genes in the human genome (Naidoo et al., 2011) and the transcriptional diversity of the human genome is not yet fully characterized (Yamashita et al., 2011).
Furthermore, RNA-Seq data is susceptible to "batch effects" and other confounders that can occlude real biological effects or, worse, mislead the un-skeptical researcher. Thus, appropriate analysis of these data requires advanced algorithms and sophisticated statistical methods, coupled with traditional scientific skepticism, to uncover biological insight buried in the data.
These difficulties dissuade many from attempting RNA-Seq and ChIP-Seq analysis, particularly those lacking previous data analysis experience, but the genomics community needs more scientists who can adeptly analyze RNA-Seq and ChIP-Seq data. Moreover, shared understanding of RNA-Seq and ChIP-Seq analysis will produce higher quality discourse between the biologists who are responsible for conducting RNA-Seq and ChIP-Seq experiments and the bioinformaticians who are experts at analyzing the resulting data produced from these experiments. Several well-developed pipelines currently exist for processing RNA-Seq and ChIP-seq data from start to finish ( ; however, these pipelines are generally designed for advanced bioinformaticians who often have existing practical experience in analyzing high-throughput data. A pipeline designed to teach those with little experience how to analyze high-throughput sequencing data is therefore needed. Thus, we developed a proof-of-concept, well-documented "tutorial pipeline" over the course of a three-day NCBI-sponsored hackathon intended to teach RNA-seq and ChIP-seq analysis to beginners. This tutorial pipeline, "SeqAcademy," incorporates state-of-the-art RNA-Seq and ChIP-seq analysis tools into a simple, easy to use workflow tutorial and we demonstrate its use with publicly available data.

Implementation
SeqAcademy uses self-contained tutorials, which runs Python, R, and Bash scripts among others, all from the document itself. It requires about 16 GB of memory storage. The tutorial files facilitate open science and reproducible code by mixing code chunks with notes and markup. This format, known as "literate programming," is particularly amenable to teaching bioinformatics because it allows learners to follow along in the document while running each code step directly within the notebook.

Operation
The tutorial begins with an explanation of how to install necessary dependencies and select interesting data from the BioProjects browser. Alignment while streaming the data is done with HISAT2 version 0.1.6 and subsequent quality control with MultiQC version 1.5. The tutorial then splits into two separate protocols: one for RNA-seq, the other for ChIP-seq analysis.
The workflow involved setup, alignment, quality control, analysis, and visualization steps for publicly available RNA-Seq and

Amendments from Version 3
We have fixed the wording of "protein interactions" to "protein-DNA interactions," clarified the motivation of why many of the parts of the protocol are the way that they are, and added gene annotation to the pipeline. Much of the writing has been updated with more motivation in general. We have placed "Quality control" after "Alignment" in the "Setup" section.
Any further responses from the reviewers can be found at the end of the article

REVISED
ChIP-seq data sets. There are many appropriate tools available for each step of RNA-seq and ChIP-seq analysis. Our goal is to present an easy to use and understandable pipeline rather than an exhaustive list of analysis tools. For each step below, we will explain the role of the bioinformatic tool, as well as our rationale for including it in this tutorial pipeline ( Figure 1). Here, we present an overview of the steps; further details for each subsection can be found on the project's Github page.

Setup
The We note that the SraRunTables file can be adjusted to specific user data, different from the RNA-seq or ChIP-seq data sets used in this project. Thus, this lightweight, portable educational pipeline can be adapted to meet the usage needs and interests of a broad base of bioinformatics beginners and teachers.

Alignment
The purpose of alignment, map raw sequence reads to a reference genome, thereby allowing quantification of a genomic property (e.g. gene transcription in the case of RNA-seq). HISAT2 is a software program used for the alignment of raw sequence data, consisting of FASTQ files (Kim et al., 2015. We chose to use HISAT2 because it allows users to stream raw sequence data rather than downloading it to the local machine, reducing disk space and time requirements for users of the SeqAcademy educational tool -an exemplary use of "edge-computing" in bioinformatics. One disadvantage of this approach is that it requires a stable internet connection, as the aligned raw sequence files are downloaded as SAM (sequence alignment mapping) files along with the log files. Nevertheless, by choosing to use HISAT2 for alignment, we reduced required disk space and broadened the potential user base of this pipeline.

Quality control
Quality control is a critical step given that sequencing data is often of heterogeneous quality, and is a way of i) identifying outliers ii) assessing whether sequencing data is a valid measure of a genomic property To generate a quality control report about the success of the alignment, we used MultiQC (Ewels et al., 2016). MultiQC reports the number of reads mapped to one unique location, reads mapped to multiple unique locations, and reads not mapped to any location in the reference genome. MultiQC can provide reports for both RNA-Seq and ChIP-seq data. Reads mapped to one unique location have a higher confidence level of being correctly mapped, as reads mapped to multiple unique locations cannot be localized to the reference with a high degree of probability. While MultiQC is not strictly necessary for this pipeline--the plots and statistics it produces are based off of the HISAT2 alignment summary files -we chose to include it to introduce users to a useful tool that is built for quality control.

RNA-Seq
RNA-sequencing is a high throughput method for studying gene transcription. After alignment and quality control, users convert the SAM files to BAM files with the samtools package version Afterwards, we demonstrate how to extract biological significance from these various analyses, by showing students how to visualize gene expression patterns and undertake exploratory data analysis with principal component analysis. Principal component analysis (PCA) is an unsupervised clustering method best suited for studies including multiple samples. If only one RNA-seq sample is present, PCA is not an appropriate analysis as no dimensional reduction can be performed. For multiple samples with a single condition, PCA is a valuable tool for identifying and quantifying potential batch effects. When batch effects are successfully isolated by PCA, the corresponding batch PCs may be valuable as adjustment variables (i.e. covariates) in downstream analysis. For example, including batch PCs as covariates in differential gene expression analysis can help reduce confounding by batch. For multiple samples with multiple conditions, PCA can potentially distinguish groups and determine how much transcriptome-wide variance the condition explains.
Notably, PCA offers a global picture of transcription and cannot determine which specific genes are different between conditions--individual genes are best identified via differential gene expression analysis (see DESeq2 below). Likewise, PCA may again be useful in this scenario for quantifying batch effects. Lastly, PCA of multiple samples can be used for as an additional quality control step with visual identification of outliers.
Finally, we show how to undertake differential expression analysis using DESeq2 version 1.21.0 (Love et al., 2014) and how to visualize these differences with volcano plots and experiment-specific visualizations in the R package ggplot2 version 2.2.1 (Wickham, 2009). Thus, students can learn how to quantify gene expression, answer biologically relevant questions through differential gene expression analysis, and visualize gene expression patterns.

ChIP-Seq
After alignment, we perform peak-calling to determine protein-binding locations in the ChIP-seq data. The peak-calling step of ChIP-Seq involves finding differentially binding sites between the two ChIP-Seq signals (input and immunoprecipitate

Use cases
Target audience This educational pipeline is designed for students without previous programming experience who are looking for an introduction to the acquisition, processing, analysis, and visualization of either RNA-seq or ChIP-seq data. Students of next-generation sequencing analysis may range the academic spectrum, from undergraduates to professors, all of whom share an interest in learning to analyze sequencing data. SeqAcademy also offers a useful introduction to the core steps of RNA/ChIP-seq analysis for use by bioinformatics educators who are teaching a class or mentoring students. Motivated individual learners, for instance a graduate student who is attempting RNA-seq analysis, may also benefit by working through SeqAcademy. The tutorial is completely self-contained, so users do not need to manage additional input files or tools beyond what is provided directly in the notebook document-every line of code to be run has already been written and tested. Thus, this flexible tutorial may be a suitable introduction to RNA-seq and ChIP-seq analysis for workshops, graduate school classes, or motivated individual learners. We also hope that fellow bioinformatics educators will build off of SeqAcademy to teach intermediate and advanced bioinformatics concepts and skills. The pipeline is simple and modular, so it can easily be adapted to analyze different datasets and customized to meet different user needs.

Learning objectives
The learning objectives of SeqAcademy are two-fold. The first and most immediate or practical objective is for a student to learn how to conduct the core steps of an RNA/ChIP-seq analysis, beginning with a search for publicly available sequencing data and ending with biologically meaningful results. The second objective is to foster a greater understanding of the concepts behind each step. This includes biological reasons behind why certain experiments such as ChIP-Seq and RNA-Seq are run, and the logic behind alignment, differential gene expression, and peak-calling. The tutorial pipeline is purposefully simple, as this will introduce an important component of next generation sequencing more gently, and will encourage students to build off of it to create more advanced pipelines that will meet the unique goals of the student.   Table 1 and Table 2 illustrate the sample input yeast data for RNA-Seq and ChIP-Seq, respectively. The RNA-Seq data examines aneuploidy while the ChIP-Seq data shows induction by 3-Amino-1,2,4-triazole (3-AT). Results of the principal component analysis, an unsupervised data reduction technique, of the RNA-Seq data are shown in Figure 2a. The slight clustering of the data into two different groups, euploid and aneuploid can be observed. A volcano plot is used to visualize significant differentially expressed genes between two groups, in this case euploid and aneuploid (Figure 2b). Figure 2c displays the enrichment of chromosome X for differentially expressed genes, consistent with the aneuploid sample having an extra X chromosome. Figure 3 shows an IGV screenshot of how peaks of protein-enrichment are distributed across the    yeast genome. The corresponding genes can be examined to determine proteins involved in 3-AT induction.

Conclusion and next steps
Limitations and future directions There are several limitations to take into account with this tutorial and future directions for further work. In this tutorial, we focused on using RNA-seq on "bulk" or homogenate tissue samples, as opposed to single-cell RNA-seq, which has distinct analytical considerations. Our pipeline is currently limited to only two of the various next generation sequencing analyses, and we would like to broaden the scope to also include DNA sequencing and other epigenetic sequencing protocols, such as whole-genome bisulfite sequencing. Our platform can also be developed further to incorporate more advanced features, such user interfaces for performing bioinformatics analyses from the web browser, login systems for users to keep track of their own progress, and forums and messaging systems for community feedback. We would also like to translate the pipeline into other languages to broaden its scope. In subsequent improvements, we plan to make the pipeline easily individualized for a user's own data sources by adjusting SraRunTables. Future hackathons may offer a useful setting to further improve this developing resource. Despite these limitations, SeqAcademy provides a solid starting foundation for beginners to learn the fundamentals.

Summary
We have presented a novel, standalone educational tool for two types of next generation sequencing data: RNA-Seq and ChIP-Seq data. This project offers a simple guidebook to an introductory analysis pipeline used in RNA-Seq and ChIP-Seq data. We introduced a cutting-edge bioinformatics tools frequently used for the acquisition, alignment, processing, analysis, and visualization of large-scale sequencing data and referenced further resources for continued learning. SeqAcademy meets the need for an educational analysis pipeline which can be used to teach undergraduate and graduate students with limited bioinformatics experience how to analyze publically available sequencing data.

Data availability
Use case data is available for the NCBI Sequence Read Archive Run Selector under accession numbers -SRP132584 and SRP106028 Center for Computational Biology, Flatiron Institute, Simons Foundation, New York, NY, USA I appreciate the efforts that the authors have made to make the RNA/ChIP-seq data processing easy by developing this pipeline. I can foresee the educational value of this paper. However, after reading this paper, I have several comments, especially about explaining 'why' each processing step is needed to users.

Software availability
It is not clear to me what's the meaning of protein interactions. Are you talking about protein-DNA interactions or protein-protein interactions? ChIP-seq is mainly used to identify protein-DNA interactions. If multiple ChIP-seq profiles are jointly analyzed, co-binding relationship between multiple proteins can be studied. However, this co-binding relationship is different from proteinprotein interaction. The author should say "protein-DNA interactions" instead of "protein interactions" as the latter usually refer to protein-protein interactions.
After reading the abstract and introduction, it is not clear to me why RNA-seq and ChIP-seq data are jointly discussed in this paper. Most of the reference papers are just talking about one data type. Please clarify why you choose these two. I know people can integrate both data types to identify regulatory networks (doi: 10.1093/bioinformatics/btx827) 1 . Is this your motivation? But this is not mentioned in this paper. As this manuscript is prepared to guide beginners lacking enough experience in this filed, it is necessary to tell them why these two data types are jointly discussed here. Although integrating these two data types is not the main focus of this paper, users can preprocess the raw data using this pipeline and do further downstream analysis by using other tools. Definitely, current description of each data type is necessary to let readers know each data type provides unique information.
Quality Control should be placed right after alignment as it is a general approach to check read quality and alignment rate. It is very misleading to readers that they need to get peaks or gene expression and then do QC check on BAM files.
The following part of RNA-seq data analysis should be improved "Afterwards, we demonstrate how to extract biological significance from these various analyses, by showing students how to visualize gene expression patterns and undertake exploratory data analysis with principal component analysis (PCA). Finally, we show how to undertake differential expression analysis In this paragraph, the authors say a lot of 'how', but to beginners, they need to know 'why' each step is needed.
The same issue to ChIP-seq data. Why do users need the intersection step using bedtools? If they get multiple replicates, the intersection step will return reliable peak loci. If data of multiple factors is obtained, they can use the interaction step to extract co-binding events. And gene annotation is usually needed after getting those peaks, which is missing in this pipeline.
The authors have made significant improvements to the tool in this revision, addressing a number of points raised in my initial review. However, sadly I was still unable to get the tutorial code to work myself. A number of the steps generated syntax and execution errors making progress impossible. The notebook has some statuses from the authors' previous execution steps with similar error messages, so I don't think that this is just my system (downloading the genome & building the index onwards).
There seems to be some active development work ongoing in the repository -both a `devel` workbook (better to use git branches instead of separate files), but also some duplicate code blocks in the main notebook. There also seems to be a discrepancy between the introduction text on the main website and on the repository.
In summary -I still think that the intention and design of the tool is excellent, and the manuscript text good. However, whilst the tutorial itself seems to be broken and not usable by others, I cannot approve this article. If these issues can be fixed so that it's possible for others to easily run through the tutorial, then I will be happy to approve.

Is the rationale for developing the new software tool clearly explained? Yes
Is the description of the software tool technically sound? Yes

Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool? Yes
Are the conclusions about the tool and its performance adequately supported by the findings presented in the article? Yes

Philip Ewels
Stockholm University, Department of Biochemistry and Biophysics, Stockholm, Sweden Bioinformatics training is of critical importance for the field of genomics, especially given the current shortage of experienced bioinformaticians and the utility of cross-training lab scientists in data analysis. In this manuscript, the authors describe a training tool developed to help teach people new to bioinformatics how to handle data from RNA and ChIP high throughput sequencing experiments.
Such a tool is valuable to the genomics community -it is a resource for interested parties to use in their own time and could also be very useful for those designing bioinformatics courses.

### Manuscript
The manuscript does a good job of describing the resource, as well as the reason it is helpful. In my opinion, it is longer than needed -the Introduction could be more concise and the Use cases section would feel more natural in the introduction. The figures could also easily be condensed into one and the table be a supplementary figure. Making the manuscript shorter would make it more accessible.

### Tutorial
The authors take an excellent approach to the design of their tutorial, using a Jupyter notebook and software installations from bioconda. These are appropriate, user friendly and provide an excellent starting point for well documented and reproducible analyses. Some dependencies are required, but these are well described on the tool's homepage and commonly used in bioinformatics. I attempted to run the tutorial on my laptop, but I ran into a few issues that required some debugging.
Broken steps I came across a couple of problems in the notebook: I couldn't get the MultiQC commands to work -the mix of Python variables in the bash commands failed in my Jupyter notebook. I changed the command to use the following command: !multiqc test/{"* test/".join(RNASeqoutrun)}* --quiet --outdir test/multiqc_rnaseq --force Next, the Bedtools steps didn't work -the intersect command refers to a sample SRR6703663, which does not seem to be included anywhere in the tutorial before this point. I gave up at this point, so I'm not sure about the final steps. If the authors think that I did something wrong, I'd be happy to have another go.
I recommend trying to get some colleagues to run this tutorial from scratch to get more feedback.

External scripts
My main concern with the tutorial is that much of the downstream processing is done within R scripts outside of the Jupyter notebook. For me, this sort of misses the point -it would be preferable to have this code mixed in with the tutorial without any external scripts. It should be possible to have both Python (bash) and R in the same notebook, and this would make those steps part of the exercise rather than an abstract "runDeseq.R" script that requires the user to go and find and read through.

Data subsampling
Although the authors have endeavoured to use a small dataset (Yeast, not Human), the data used is still large. The authors note in the tutorial that the alignment "will most likely take several hours.", which limits the utility of the tutorial in any teaching scenario.
It should be quite easy to use the -u/--qupto flag in HiSAT2 to limit the number of alignments performed. I imagine that most of the downstream analysis will still work with subsampled data, and will run much faster. I tried doing this with 1,000,000 reads and the alignment step completed in a couple of minutes.

### Tutorial -Minor points
Software requirements I had trouble getting the bioconda installation command to first install and then activate inside the the notebook. I would instead recommend creating the conda environment on the terminal first (including jupyter as a package to install), activating it and then running jupyter. This is what I ended up doing to get the tutorial to work.

Streaming Input Data
The approach of streaming data from the SRA for the alignment step is a little non-standard. The reasoning is described in the paper as helping to reduce the disk space footprint required for the analysis. Whilst its cool that this step works, I would argue that it lessens the value of the tutorial, as most bioinformatics projects can not use such a method. Starting with regular static FastQ files would be more typical and useful for future analyses. This would also allow the additional step of FastQC for raw data analysis, which is pretty typical for such pipelines. The downside of this is that it invalidates my above subsampling idea! But SRA data could still be streamed to FastQ files with a similar subsampling approach.