MOIRAI: a compact workflow system for CAGE analysis

Installation

MOIRAI can be installed on UNIX based operating systems (tested on MacOSX, Ubuntu and Fedora) running Apache and PHP services. To make installation easy, we prepared an automatic script to set up MOIRAI on user’s personal computer. Using it, users will be able to run a default "Hello World" test workflow included in the download package after minutes of installation. Other CAGE workflows require the installation of analysis software, reference genomes and configuration of input and output directories. We provide the tools and scripts to automate these steps.

Workflow overview

Each box in a MOIRAI workflow represents one process and an arrow between units describes a flow of the data. When a unit has multiple inputs, it will not be executed until all previous computations are completed. Furthermore, in the case of paired-end sequencing data, MOIRAI will verify that input files belong to the same sample by matching the file names. All units are colored to reflect three basic types of entities. Input files and reference databases are colored green, intermediate steps producing temporary files are gray and finally all outputs are shown in blue. For quality control purposes, we provide graphs embedded within the workflow to summarize the results at key steps. Similarly, thumbnails of windows looking similar to unix terminals allow users to check summary statistics and look at the most frequent sequences.

Importing tools

To add flexibility, MOIRAI contains a standardized web-based mechanism to incorporate new command line tools and methods into the system. In brief, users can make an entry into an online table. Each entry requires the user to set the name of the program, the in- and output file formats, options and default parameters. It is possible to include the same program several times but with different parameter options to customize the resulting unit for particular types of data. Once incorporated these units can be used to create new workflows or modify existing ones.

Available tools

CAGE tags from different libraries are tagged with barcode sequences and pooled into one sequencing lane to improve sequencing capacity and to enable accurate comparison between samples [12]. We wrote tools to demultiplexing these pooled tags into separate outputs according to barcode sequences for both single-end and paired-end CAGE. The MOIRAI workflow system will automatically recognize that demultiplexing units produce multiple files and will apply all downstream steps to all the files.

After demultiplexing we use two tools to remove unwanted sequences. TagDust [13] removes reads similar to primer dimers and other artifacts while rRNAdust removes reads matching ribosomal RNA sequences.

The tome database system is a simple implementation of a compressed sparse row data structure to store mapped CAGE reads from multiple experiments. The database can be queried to extract the number of CAGE reads within genomic boundaries given in bed6 format. The latter can be used to calculate the fraction of reads in each library mapping to promoters and to cluster biological samples.

SAMstat [14] displays mismatch, insertion/deletion error profiles, mapping rates and many other useful statistics of FASTQ/SAM/BAM sequence reads in HTML format. The report helps users to spot biases and problems in sequencing runs and protocols.

In addition we included several generic plotting tools to display the base distribution in intermediate sequence files and to list the most frequent sequences. These tools are extremely useful to spot and correct problems in complicated pipelines.

Summary of results

Once a workflow is completed, MOIRAI summarizes the results and additional information on the run itself in a web page. The latter includes paragraphs describing each step including the version and parameters of the software used. By default all standard error and output messages are recorded per process and are stored in a separate html page. The run time and resulting file sizes of each step can be viewed by moving the mouse over the corresponding box in the workflow. Finally, the main output files are accessible by clicking on blue output boxes.

In addition to the results, we include a copy of the workflow template used for the computation itself. The workflow can be imported into another user’s account or MOIRAI instance to precisely reproduce all the results.

Edit workflow

We included a workflow editor that can be accessed via the web-interface or as a stand alone application. The advantages of using the stand alone application include dragging and dropping data files directly from the desktop into the editor, offline editing of workflows as well as the ability to selectively execute only the recently added steps. The latter is desirable when debugging large workflows with time consuming initial processing steps. The mechanism for incorporating new tools and setting parameters is simple.

CAGE mapping workflow

The mapping workflow is organized into three basic tasks: (a) raw reads are demultiplexed and trimmed, (b) artifacts and reads corresponding to ribosomal RNA are filtered out and (c) the actual mapping of the remaining reads. We placed quality control units on the left and right hand side of the workflow (Figure 1).

A multiplex sequence and a "CAGCAG" linker sequence at the beginning can be readily identified in the top left plot showing the base distribution of the raw reads. After demultiplexing and read trimming the base distribution is much more even (middle left plot). By default, we remove reads containing ambiguous bases (N). We can see from the second plot on the left that these low quality bases (shown in black) accumulate at the 3’ end of the reads which is typical for the Illumina sequencer [15].

Finally, a bar graph at the end of a workflow summarizes how many reads passed each stage. In the example shown we conclude that the CAGE reads are of high quality since only 2.4% of sequences were removed by filtering steps, and more than 90% of the filtered tags were successfully mapped.

Based on our experience, we consider a CAGE library to be of acceptable quality if fewer than 10% of the reads correspond to ribosomal RNA and if the mapping rate higher than 70%.

Annotation workflow

The library annotation workflow takes mapped CAGE reads from multiple libraries and annotates them according to user specified rules (Figure 2). Here, we obtained gene models from UCSC [16] and used intersectBed to determine the fraction of CAGE reads at 100bp upstream, 5’ UTR exon, coding exon, 3’ UTR exon, intron, and 100bp downstream categories. We intersect the CAGE data with the annotation files hierarchically to make sure each CAGE read is only annotated to a single category. All remaining tags are categorized as intergenic. The results are summarized in a tab delimited table and shown as a graph in the workflow. While being designed with CAGE data in mind the same mechanism can be used to create any annotation pipeline based on any genome annotation.

In the case of the ENCODE CAGE libraries, it is evident that libraries obtained from nuclear sub-compartments have fewer reads mapped to promoters (100bp upstream and 5’UTR category) compared to the other libraries. As expected, the polyA+ libraries have high percentage of promoter regions compared to nuclear compartments as most measured RNA molecules are messenger RNAs.

Expression analysis workflow

Finally, the expression analysis workflow groups mapped reads into peaks using a parametric clustering algorithm implemented in the program Paraclu [17]. In brief, this method reports genomic intervals containing many more CAGE reads than surrounding regions. These regions can be contained within each other giving rise to hierarchies of clusters. From this we select all clusters of length ≤ 200bp and with a stability greater than two. Normalized expression values for these peaks in multiple samples is extracted using the tome program.

The tome database itself is saved and can be used independently of the computed workflow. Based on the expression within the CAGE peaks the workflow calculated the Pearson’s correlation between samples and plots the results as a dendrogram and heat map (Figure 3).

The dendrogram graph is drawn by pvclust where probability values (p-values) are computed for each cluster with bootstrap resampling techniques [18]. If approximately unbiased (AU) p-value and bootstrap probability (BP) value are not necessary, bootstrap can be switched to computationally quicker hierarchical clustering (hclust). Biological replicas cluster together and libraries obtained from the whole cell and cytoplasmic fraction cluster closer together than nuclear libraries. One cytosolic polyA+ replica is clustered far away from other two replicate libraries.

The heatmap is another way of visualizing the overall similarity among samples. The graph shows two clusters of nuclear sub-compartments and cytoplasmic fractions with nuclear polyA+ replicates in the middle. It is interesting to note that nucleus polyA+ libraries appear to be similar to both nuclear and cytoplasmic libraries. Light green patterns appearing between all polyA+ samples indicates large number of mRNAs with polyA+ in nucleus are transported out to cytosol.