This is the simulation function of the package PROPER; it carries out simulations for a specified number of runs and for a given sample size (two groups of equal size). At each run, a dataset which has the characteristics of the data being investigated is generated and analyzed, i.e., transcripts differentially expressed between the two groups are detected, using one of the differential expression (DE) analysis methods currently available in PROPER (edgeR, DSS and DEseq). After each run the names of the transcripts differentially expressed are returned along with their p-value, log fold change and FDR.

The extensions we are proposing consist of changing the argument/parameter that specifies the DE method to add EBSeq to the list of methods and altering the core of the function to enable DE analysis using EBSeq. Since for isoform analysis, unlike gene-level analysis, isoform names are required in addition to the gene name, we inserted a check to ensure this requirement is met to prevent the function from crashing if input tables similar to those for edgeR, DEseq, and DSS are provided while EBSeq is specified as DE analysis method.

In the implementation of EBSeq we enabled within runSims, transcripts significantly differentially expressed are those with PPDE ≥ (1–alpha) where alpha (type I error rate) is set to 0.05; a value that can be changed in the code we provided. In this setting significance is based on alpha so that all transcripts with FDR ≤ alpha are considered significantly differentially expressed and FDR is similar to the idea of bonferroni correction. However it is important to note that EBSeq does not require extra multiple test adjustment since it uses a single large model to test all transcripts simultaneously rather than a model testing one gene at a time (Leng et al., 2013). Alpha is also called target FDR which in the setting we described is based on a hard threshold. EBSeq offers the possibility to compute a soft threshold via the function crit_fun described in its R manual.

This is the function we wrote to carry out EBSeq isoform analysis within runSims. It is analogous to the other three functions that can be called to run DE analysis using egdeR, DEseq, and DSS. The function analyzes the simulated isoform expression data and outputs a matrix with all the transcript, p-values, log fold changes, FDR values, and other statistics. This information is then subsequently passed on respectively to comparePower, powerPlot and plotFDR to evaluate power and generate graphs. comparePower is central to PROPER as it is the function that computes power by evaluating the success rate over the multiple simulation runs. Success rate is determined by comparing the p-values to the specified type I error rate (alpha) to reject to null hypothesis. More details about this function are available from PROPER tutorial document.

Our changes to the function powerPlot consisted mainly in amendments to allow for the plotting of just one curve instead of multiple curves (one per sample size considered) and minor additions related to the visualization of the curves. The function plotFDR is similar to powerPlot save it plots FDR curves. There is a function of the same name in PROPER and hopefully, the code we generated can be integrated into that function to enable the additional features we introduced.

The aim of this short guide is to demonstrate the use of the new extensions to compute statistical power for isoform expression data using the EBSeq approach. An excellent PROPER tutorial is available online at: https://www.bioconductor.org/packages/devel/bioc/vignettes/PROPER/inst/doc/PROPER.pdf

This documentation along with the R manual of the PROPER package provides a detailed explanation of the parameters we set in the code on this document. Therefore, the arguments required by PROPER functions used in this tutorial will not be explained with great details. Similarly we refer the user to EBSeq documentation (https://www.rdocumentation.org/packages/EBSeq/versions/1.12.0) for a deeper understanding of that method.

The data used in this guide (see Supplementary Material Datasheet 1) and the scenario investigated are just meant for demonstration and are by no means biologically meaningful. They are only meant to help understand the formats of the input data and the steps to evaluate statistical power and visualize the results.

It is important to note that a real world analysis typically includes more than 10,000 genes and a PROPER analysis of such large data can be computationally slow. For example an analysis of 15,000 genes undertaken with 50 simulation runs on a machine with 8GB memory took 2 hours to complete. This is because, in addition to the iterative process of PROPER, each model fit in EBSeq takes 2 or more iterations to reach convergence. The high computational cost represents a limitation; but we believe this weakness is worth the more realistic results achieved through this approach.

The data simulated and analyzed is a matrix where the first 2 columns of the table hold isoforms and transcripts IDs. The process is stopped and a message printed on screen if the information in these 2 columns is missing while the user has specified EBSeq as DE analysis method. For the extensions presented in the manuscript to run, the names of the first two columns must be “transcript_id” and “gene_id” and the IDs in the first column must be unique; here we used UCSC IDs but other names such as RefSeq IDs are also valid. The rest of the columns holds sample IDs. As shown below, rows store isoform expression values for the samples in the columns.

The toy data used in this demonstration (Table 1) consist of the expression of 208 isoforms across 8 samples (2 groups: 4 cases and 4 controls). The first 4 samples (ID1 to ID4) represent to control group whilst ID5 to ID10 are the case group. PROPER requires an equal size for the 2 groups (e.g., cases and controls).

Along with this table we provide 3 variables (AveLogCPM, logDispersion, libCount) which represent respectively the log average Count Per Million [average expression of each isoform, CPM = count/sum (counts) x 1 million] log taqwise (genewise) dispersion (dispersion value for each isoform) and the depth of the data (the total count for each sample). This information gives the characteristics of “real” data (actually some hypothetical data) we are analyzing. As already mentioned these parameters are explained with ample details in the PROPER tutorial and Manual. For this tutorial we obtained the characteristics of our data from a pilot analysis but as explained by Wu et al. one may also set these parameters based on literature or other sources.

In PROPER, power analysis is carried out in 3 steps: First the simulation scenario is set up then the simulated data are generated with the characteristics of the real data and finally the power is evaluated.

Let us begin by loading the R libraries and functions required. We load the package edgeR because we use the Trimmed mean of M-values normalization (TMM) method embedded in that library to normalize the data in the EBSeq analysis step.

Now we load the characteristics of our “real” data and the estimated log fold change of the differentially expressed isoforms which represent a proportion (propDE) of 6% of all the isoforms.

In a real world power analysis it is probably better to obtain these characteristics and a rough estimate of the number of differentially transcripts and their fold change difference from the real data by running a preliminary DE analysis. In such case the analyst can:

Set the simulation i.e. generate an object that stores the parameters of the simulations.

We now can run the simulation, generating n datasets, where n is the number of simulation runs, and analyzing each dataset. We set the number of runs to 20 but in real world example one should select at least 50 runs and the larger the number of runs the more reliable the results are. The argument “Nreps” represents the number of case-control pairs (4 in the toy data). The function returns p-values and FDR values for each simulation run.

Finally power is evaluated across all the runs and the results graphically displayed. The arguments and output of the below function are described in details in the R manual for PROPER; briefly: The function returns overall and marginal values for power, type I error rate and FDR. In this example we used only the arguments relevant for the isoform analysis process ran by calling run.EBSeq from within this function; DE significance is based on FDR and the cut-off for significance is 0.05. We assess power for each expression strata (stratify.by = “expr”) and all the expression strata we specified (strata) are considered (strata.filtered = 0).

Finally we visualize the results by generating a graph with 2 plots (Figure 1), one for power and one for FDR. If the parameter “error bar” is set to “TRUE” uncertainty bars are displayed around the estimates in the plot. For this demonstration we used a very small dataset which results in large uncertainty that might cause R to throw a warning message if the argument is set to “TRUE.”

The code in this tutorial was run using R version 3.2.1, PROPER version 1.0.0, EBSeq version 1.10.0 and edgeR version 3.8.6.

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.