Automated analysis of small RNA datasets with RAPID

Basic module

The first of four RAPID modules is rapidStats, which performs sequence (FastQ) alignment, with or without contaminant removal, using Bowtie2 (Langmead, Salzberg & Langmead, 2013). After alignment, RAPID obtains read statistics such as read length distribution, soft-clipped nucleotides, strandedness, and nucleotide content. RAPID can skip the alignment and directly use alignment files (BAM/SAM) as well. To efficiently process, capture and store the aforementioned statistics, RAPID uses SAMtools (Li et al., 2009), BEDtools (Quinlan & Hall, 2010), and custom Perl, Shell, and R scripts. The statistics captured by this module serve as input for other modules.

Normalization module for multi-sample comparison

RAPID aims to facilitate an unbiased comparison of genes or regions across multiple sRNA samples. Other than the sequencing depth itself, sRNA studies pose an additional challenge during normalization. For instance, to understand RNA interference (RNAi) mechanisms and how the siRNA homeostasis is maintained, often a gene or siRNA region is knocked down. One such knockdown strategy is to introduce large amounts of siRNAs, called primary siRNAs, against the knockdown gene or any siRNA region. Consequentially, secondary siRNA production is triggered by the primary siRNAs. These primary and secondary siRNAs, which are also sequenced, can add up to millions of reads in the total library size.

To our knowledge, there are no normalization methods specialized for knockdown based small RNA-seq studies. However, many methods have been proposed to normalize mRNA-seq data, which can be broadly categorized in two classes: (i) total count scaling (TCS) methods and (ii) methods which utilize quantities like median log-fold change, among all genes between mRNA-seq experiments. To be able to use the latter methods, sRNA loci annotation should be available, and should assume that most of the sRNA loci between samples are not differentially expressed. In model organisms like Paramecium tetraurelia, little is known about the localization, and expression variability of endogenous small RNA loci. Hence, the second class of methods may not be applicable. However, the disadvantage of TCS methods is that the used normalization factors were shown to be biased by highly expressed genes in the dataset (Dillies et al., 2013). In case of knockdown samples, TCS methods will be heavily biased because of the millions of primary, and secondary siRNAs associated with the knockdown gene or region.

In mRNA-seq data, a variant of TCS method (Sultan et al., 2008) was introduced, where normalization is achieved by scaling through a factor that estimates the difference in the number of reads mapped between samples. We previously proposed a variant of the TCS method in a knockdown based siRNA study (Götz et al., 2016). Here, we term this variant as KnockDown Corrected Scaling (KDCS) method, where we remove from the estimated total library size, all small RNA reads that map against the knockdown genes, this quantity is denoted as K below. Assume read count R for a region of interest that we want to compare between samples. T is the total number of reads mapping to the genome, and K is the number of small RNA reads mapping to the knockdown gene. We compute the normalized read count $\hat{R}$: (1)${\displaystyle \hat{R}=R\cdot \frac{M}{T-K},}$ where M is the maximum over all values (T₁ − K₁), …, (T_n − K_n) over all n samples. RAPID uses the KDCS method, by default. Hence, in the absence of knockdown genes, the normalization works as the normal TCS method. However, in order to provide flexibility with the choice of normalization for knockdown free analysis,we have also incorporated size factor-based normalization from DESeq2 (Love, Huber & Anders, 2014). If an user can safely assume that most of the genes or regions between samples are not differentially expressed, in a small RNA based study, then they can use the DESeq2 normalization.

Visualization module

As visualization enables better understanding of data, the rapidVis module of RAPID automatically generates insightful plots from the output of previous modules. RAPID makes use of Rmarkdown (http://rmarkdown.rstudio.com) to create easily navigable HTML reports. This module contains two modes: statistics and comparison mode. The statistics mode accepts input from the rapidStats output file, and provides various single category plots detailing on the distribution of read length, strandedness, soft-clipped nucleotides, and coverage plots for each gene/region analyzed. In addition, this report also provides combinations of the aforementioned properties. For instance, how does strandedness differ across different read lengths. Comparison mode accepts the rapidNorm analysis output file, to equip the user with qualitative reports (Heatmaps, PCA, MDS) of samples. Further, sample and gene/region wise comparison plots of the properties inherent in the data. All plots are shown both in normal and log scale such that the user can directly incorporate them into publications.

Differential analysis module

Differential Expression (DE) analysis is one of the common downstream analysis in comparative studies. RAPID equips the user with this functionality by incorporating the DESeq2 package. Upon invoking the rapidDiff module, raw counts are utilized from the output of the rapidStats module to perform DE analysis, with default parameters of DESeq2. Results of the DE analysis include intuitive plots (such as MA Plot, Heatmap, PCA) and the list of DE genes/regions.

Usage and availability

We strongly recommend using RAPID from https://anaconda.org/bioconda/rapid as a conda recipe. However, it can also be freely accessed from https://github.com/SchulzLab/RAPID. A detailed use case based documentation is provided at http://rapid-doc.readthedocs.io/en/latest/.

Definitions

Here we describe the formula used in the case studies.

Coefficient of variation (CoV): Coefficient of variation is the ratio of standard deviation to mean of the data set. For a gene or region of interest, i, with n samples (2)${\displaystyle Co{V}_i=\frac{{\sigma }_i}{{\mu }_i},}$ where σ_i, and μ_i are the standard deviation, and mean of the gene or region of interest i in the n samples respectively.

AntiSense Ratio (ASR):

Antisense ratio is the ratio of the number of antisense reads to the total number of reads in a gene or region. If R, and AS are the total, and antisense read counts of a region of interest, i, respectively, antisense ratio is calculated as follows: (3)${\displaystyle AS{R}_i=\frac{A{S}_i}{R_i}.}$

DatasSets

We show the application of RAPID, using two different datasets which are briefly described below.

Comparison of Paramecium tetraurelia serotypes

P. tetraurelia is a unicellular, free-living ciliate commonly found in fresh-water lakes. They show nuclear dimorphism and have a wide range of phenotypes. One such phenotype that depends on an epigenetically controlled, mutually exclusive expression of members of a multigene family is called serotypes (Cheaib et al., 2015). With a dimorphic nucleus and more than 11 serotypes, P. tetraurelia sRNAs are involved in regulatory mechanisms at the post-transcriptional, and epigenetic level (Götz et al., 2016; Cheaib et al., 2015; Carradec et al., 2015).

Our first example is an analysis on four sRNA-seq datasets (ENA: PRJEB25903) from wildtype serotypes (51A, 51B, 51D, and 51H) of P. tetraurelia. We were interested in sRNAs produced in the rDNA cluster producing 17S, 5.8S, 25S ribosomal RNAs. A simple genomic visualization of the different components of rDNA cluster regions we quantified can be seen in Fig. 2. We can observe from the strand-specific read distribution plots in Fig. 2 that the regions, namely External Transcribed Spacer (ETS), Internal Transcribed Spacer 1(ITS1), and Internal Transcribed Spacer 2 (ITS2), which get excised in the processing of polycistronic pre-rRNA, accumulate 23nt antisense small RNAs. It is known from yeast, that rRNA maturation involves co-transcriptional endonucleolytic cleavage and highly concerted trimming events to subsequently process the final rRNAs (Henras et al., 2015). Our data here suggests that in P. tetraurelia these elimination processes are associated with antisense siRNAs, possibly produced from RNA-dependent RNA Polymerase activity.

Analysis of our case studies, and their respective figures can be reproduced with the help of ( Data S1).

Normalization case study in Paramecium tetraurelia knockdowns

One of the unique features of RAPID is the KDCS normalization that can correct for the excess of sRNAs introduced in knockdown experiments in experimental approaches utilized in many diverse organisms. To demonstrate the effect of KDCS normalization, we utilized the ICL knockdown data sets from the study by Götz et al. (2016). This study investigates the molecular mechanisms of different sets of trans-acting RNAi components in P. tetraurelia. ICL is a gene in P. tetraurelia, which is not involved in the RNAi machinery. In the original study, as a control, the ICL gene is knocked down by introducing primary siRNAs against ICL (see ‘Methods’). In our work, we quantified the sRNA read counts of four example sRNA regions (which are in the original study different constructs of the ND169 gene) from the mentioned datasets. In this setup we expect that all datasets behave the same, as these are biological replicates of the same system.

As described earlier, very little is known about the localization, and expression variability of endogenous small RNA loci in P. tetraurelia. Therefore, in these sRNA knockdown samples, normalization methods, such as DESeq2, may be inappropriate, due to the assumption that the majority of regions are unchanged. We compared the effect of the TCS, DESeq2 and KDCS normalization approaches to using no normalization. We used the coefficient of variation (CoV) to measure the performance of the normalization method (see Methods; Eq. (2)) as normalization should reduce the variance in read count per region. A smaller CoV suggests a better performance of the normalization method.

Figure 3 shows the CoV values of the raw, and normalized sRNA read counts, for four example regions that had been studied by Götz et al. (2016). We can observe from Fig. 3, that the KDCS method performs better in all the regions, compared to the generic TCS method. It also performs as good or slightly better than the normalization of DESeq2 for this example. All normalization approaches are better than using no normalization, which strongly argues for their use. This experiment suggests that our KDCS method is a better alternative to the TCS method and is applicable when few regions are known.

Analysis of Schizosaccharomyces pombe knockdowns

The fission yeast, or Schizosaccharomyces pombe, is another widely studied unicellular eukaryotic model organism, where RNAi pathways are prevalent. We explored the time point datasets of WT, and three different knockdowns of S.pombe from the study by (Joh et al., 2016). In this study, the authors investigated quiescence associated changes in small RNA transcriptomes and epigenetic modifications to identify the key players involved in quiescence. They also examined the role of RNAi proteins, by knocking out three of them, namely, Ago1, Clr4, and Dcr1. One of the key findings was that during quiescence in S.pombe, a set of sRNA-enriched genes were identified as crucial elements for the survival of the organism (Joh et al., 2016). We explored these sRNA-enriched genes using RAPID, to demonstrate the discovering potential of RAPID.

With a simple RAPID analysis comparing the different knockout samples, it was easy to screen for interesting properties in the data. We discovered that a subset of these sRNA-enriched genes have relatively higher antisense ratio (ASR, see ‘Methods’ Eq. (3)) (Fig. 4). This increased ASR was observed in different sets of genes in various knockouts. The subgroup of genes with higher ASR, could play a cis-regulatory role to silence the genes. While further investigation is necessary, our RAPID analysis (Fig. 4) suggests an involvement of different sRNA mechanisms in ensuring the survival of S.pombe in quiescence.