aaaffiliationtext: Merck & Co., Inc., Kenilworth, NJ, USAbbaffiliationtext:
MRC Biostatistics Unit, School of Clinical Medicine, University of Cambridge,
UK

# Online FDR Control for RNAseq Data

Lathan Liou  Milena Hornburg  David S. Robertson

###### Abstract

Motivation: While the analysis of a single RNA sequencing (RNAseq) dataset has
been well described in the literature, modern research workflows often have
additional complexity in that related RNAseq experiments are performed
sequentially over time. The simplest and most widely used analysis strategy
ignores the temporal aspects and analyses each dataset separately. However,
this can lead to substantial inflation of the overall false discovery rate
(FDR). We propose applying recently developed methodology for online
hypothesis testing to analyse sequential RNAseq experiments in a principled
way, guaranteeing FDR control at all times while never changing past
decisions.

Results: We show that standard offline approaches have variable control of FDR
of related RNAseq experiments over time and a naively composed approach may
improperly change historic decisions. We demonstrate that the online FDR
algorithms are a principled way to guarantee control of FDR. Furthermore, in
certain simulation scenarios, we observe empirically that online approaches
have comparable power to offline approaches.

Availability and Implementation: The onlineFDR package is freely available at
http://www.bioconductor.org/packages/onlineFDR. Additional code used for the
simulation studies can be found at
https://github.com/latlio/onlinefdr_rnaseq_simulation

Contact<EMAIL_ADDRESS>

Keywords: RNAseq, differential expression, false discovery rate, multiple
hypothesis testing

##  1 Introduction

RNA sequencing (RNAseq) is a powerful and widely used tool to profile the
expression of thousands of transcripts or genes in parallel, allowing
scientists to study a biological system at the transcriptome level. Most
RNAseq experiments aim to compare two biological states such as a treated and
untreated system to better understand the biological differences between the
two states. At the core of this analysis is the differential expression
analysis that compares the expression of genes between those two states. Since
hundreds to thousands of genes are simultaneously tested in high-throughput
experiments for their differential expression between two biological states,
multiple hypothesis testing correction is essential to control the rate of
type I errors (i.e. a false positive differential expression result). Earlier
methods aimed to control a metric known as the family-wise error rate (FWER),
but it is often highly conservative, controlling the probability of any false
positives at the cost of greatly reduced power to detect true positives
[Kor+19]. The false discovery rate (FDR), or expected proportion of
discoveries which are falsely rejected [BH95], was more recently proposed as
an alternative metric to the FWER in multiple testing control. FDR control has
become a largely accepted standard in genomics [CRS10]. Standard RNAseq
analysis approaches like DESeq2 and limma typically rely on the classic
Benjamini–Hochberg procedure to control the FDR [LHA14] [Rit+15].

While the analysis of a single RNAseq dataset has been well described in the
literature, the widespread application of sequencing in modern research often
results in a temporal accumulation of potentially related RNAseq datasets. If
datasets have the same comparison groups of interest and sufficient metadata
to adjust for any confounding differences in experimental or study design,
they may be combined in some way to pool information. The simplest and most
widely used analysis strategy is to just ignore the temporal aspect and
analyse each RNAseq dataset separately. However, this can lead to substantial
inflation of the overall FDR [Zrn+21]. Intuitively, we want to control FDR
over all the RNAseq datasets that are related in some way (e.g. transcriptome
of the same biological population). Open-source repositories of high-
throughput gene expression data such as Gene Expression Omnibus are examples
of growing datasets where similar RNAseq datasets are collected [EDL02].

An alternative analysis strategy is to pool together the differential
expression results of the RNAseq datasets as they become available, and then
apply classical methods of FDR control on the _p_ -values. However, the
critical disadvantage of this approach is that the decisions (i.e. the
hypotheses rejected) for a particular dataset may change over time as new
datasets are added, which is clearly undesirable in most application areas
(see Section 3).

In light of these limitations, we propose applying recently developed
methodology for online hypothesis testing to analyse sequential RNAseq
experiments in a principled way, guaranteeing FDR control at all times while
never changing past decisions. In the framework of online hypothesis testing,
the investigator must decide whether to reject the current null hypothesis
without knowing the future $p$-values (or even the total number of hypotheses
to be tested), but only knowing the historic decisions to date. Recent
proposals have extended this framework to allow batches of hypotheses to be
tested in this way, which we exploit in our context of analysing RNAseq
datasets. Batches are simply sets of _p_ -values; for instance, a differential
expression analysis of one RNAseq dataset generates a batch of _p_ -values,
and the analysis of another RNAseq dataset generates another batch of _p_
-values. Robertson et al. describe scenarios in which online hypothesis
testing can be applied to growing public biological databases such as the
International Mouse Phenotype Consortium [Rob+19].

Through a simulation study of RNAseq datasets, we compare and contrast the
aforementioned three analysis strategies, showing that the online hypothesis
testing algorithms indeed control the FDR. Furthermore, in certain simulation
scenarios, we observe empirically that online approaches can even be
competitive with the other approaches in terms of power. We also illustrate
the proposed methodology on three real-world RNAseq datasets arriving over
time.

##  2 FDR Control Paradigms

In this paper, we compare three different FDR control methods and distinguish
between “offline” and “online” methods (Figure 1). We also want to delineate
the vocabulary we use in this paper: “dataset” or “batch” refers to a set of
differential expression _p_ -values (note that this is different from control
or experimental “samples” in a standard two-group differential expression
comparison experiment); “family” refers to a group of datasets.

Figure 1: Diagram of offline, pooled offline, and online FDR control
approaches

We use the term “offline” to refer to methods that take as input a single set
of data and output a set of rejection decisions for all hypotheses at once to
control the FDR. The offline approach assumes no knowledge of previous or
future datasets and only considers the dataset obtained at any given time
(Figure 1). This approach is what is widely performed in academia and
industry; we note that “offline” is not a common designation, but we describe
this conventional approach as such to distinguish it from “online” method,
which we describe later.

In our paper, for the offline approaches, we use the Benjamini-Hochberg (BH)
procedure [BH95] and Storey-BH (StBH) procedure [Sto02]. Given p-values
$p_{1},p_{2},\ldots,p_{N}$ and a significance level $\alpha$, the BH procedure
is as follows:

  1. 1.

Order the _p-_ values from smallest to largest, giving ordered _p-_ values
$p_{(1)}\leq\cdots\leq p_{(N)}$ and define $p_{(0)}=0$;

  2. 2.

Let $i^{*}$ be the maximal index such that $p_{(i^{*})}\leq i\alpha/N$;

  3. 3.

Reject null hypothesis $H_{j}$ for every test with $p_{j}\leq p_{(i^{*})}$.

StBH improves upon the original BH procedure by letting users choose a
parameter $\lambda\in(0,1)$, which is used to estimate the proportion of nulls
in the _p-_ value set, defined as:

$\frac{1+\Sigma^{N}_{i=1}\bm{1}\\{p_{i}>\lambda\\}}{N(1-\lambda)}$

The motivation for using StBH is that BH might be overly conservative when
there are many expected rejections with a strong signal, so StBH is more
adaptive to the set of _p-_ values at hand, which will provide important
intuition to the online analogues of these algorithms. It is important to note
that BH relies on the assumption that the _p-values_ are positively dependent
[YD01] and StBH requires independent _p_ -values [Sto02].

We use the term “pooled offline” to describe an extension of the classical
offline approach in which multiple datasets are appended into a single pooled
dataset and classical FDR control methods (e.g. BH are then applied to the
pooled dataset; Figure 1). The pooled offline approach may be expected, at
first glance, to properly control FDR across all datasets since it takes into
consideration previous hypotheses from prior datasets; however, we demonstrate
later that this is not the case.

To illustrate the pooled offline approach, let’s consider a toy example where
we have 5 datasets, which have arrived over time, each containing $N$
hypotheses. Each dataset $d$ contains $N$ _p-_ values
$p_{d1},p_{d2},\ldots,p_{dN}$. If we applied a pooled offline approach at the
time the fifth dataset arrived, we combine the 5 sets of _p-_ values
$p_{11},p_{12},\ldots,p_{21},\ldots,p_{51},\ldots p_{5N}$ and then apply
either the BH or Storey-BH procedure as described above. Note that the pooled
offline approach applied to a single dataset is exactly equivalent to the
offline approach. The important assumption with offline methods is that they
require knowledge of all the _p-_ values, which does not address the situation
where we do not know the total number of hypotheses or future _p-_ values.

In contrast, the online approach allows for informed hypothesis-driven
decision-making using prior datasets and offers guaranteed FDR control for
potential future datasets. With offline methods, there are distinct points in
time when an FDR control method is applied, but with online methods, the
application is continuous throughout the timeline of datasets: past, present,
and future (Figure 1). We describe the online paradigm for testing batches of
$N$ hypotheses arriving sequentially over time, where each batch represents a
single RNAseq experiment. At each time $t\in\mathbb{N}$, a batch of $N_{t}$
hypotheses arrives. In the online paradigm, no information about future
batches needs to be available, such as how many hypotheses each batch contains
or the total number of batches. We denote the set of hypotheses in the _t-_ th
batch by $\bm{H_{t}}:=\\{H_{t,1},\ldots,H_{t,N_{t}}\\}$. Each hypothesis has a
$p$-value associated with it. Let $\bm{P_{t}}$ denote the set of _p-_ values
corresponding to the _t-_ th batch of hypotheses, given by
$\bm{P_{t}}:=\\{P_{t,1},\ldots,P_{t,N_{t}}\\}$, where $P_{t,j}$ is the _j-_ th
_p-_ value in batch $t$. The FDR _up to the time t_ is then the expected false
discovery proportion (FDP) up to time $t$, i.e. the expected proportion of
false discoveries up to time $t$ among the rejected hypotheses up to time $t$,
for all batches $s$ up until time $t$. Letting $V(t)$ denote the number of
incorrectly rejected hypotheses in batch $t$ and $R(t)$ denote the number of
rejected hypotheses in batch $t$, the FDR is defined as follows:

$\displaystyle\text{FDR}(t)\equiv\mathbb{E}[\text{FDP}(t)]\equiv\mathbb{E}\left(\frac{\sum^{t}_{s=1}V(s)}{\max(\sum^{t}_{s=1}R(s),1)}\right).$

The batch procedures for online control FDR were originally described by Zrnic
et al. (2021) [Zrn+21] and implemented in the R package _onlineFDR_ [Rob+19]
which is freely available through Bioconductor at
http://www.bioconductor.org/packages/onlineFDR.

In brief, batch procedures aim to interpolate between the offline and online
settings and achieve high power, while guaranteeing FDR control at all times
$t\in\mathbb{N}$. Batch algorithms adaptively determine a _test level_
$\alpha_{t}$ based on information about past batches of hypothesis tests. We
use three Batch procedures: BatchBH, BatchStBH, and BatchPRDS. BatchBH and
BatchStBH are the online analogues of the offline BH procedure and Storey-BH
procedure respectively. Both BatchBH and BatchStBH assume that _p_ -values are
independent. BatchStBH aims to be more adaptive to the data than BatchBH,
given its requirement for a user-chosen constant $\lambda\in(0,1)$. BatchPRDS
is a special case of BatchBH that controls the FDR when the _p-_ values in a
single batch are positively dependent yet independent across batches. Further
details about these algorithms can be found in the original paper [Zrn+21].

##  3 Real World Data Application

We used three real-world datasets procured from internal studies of Merck &
Co., Inc., Kenilworth, NJ, USA. Datasets are anonymized as Batch 1, 2, and 3.
All three datasets were RNAseq studies of CT26 syngeneic mouse models treated
with murinized rat anti-mouse PD-1 antibody (muDX400) or vehicle control with
the same experimental design such as treatment duration and dose. The order of
datasets was determined by the chronological order that they were performed
over time at Merck & Co., Inc., Kenilworth, NJ, USA. We selected vehicle
control samples as our control arm and muDX400 samples as our treatment arm in
a two-sample differential expression analysis. Batch 1 contained 10 control
and 10 treatment samples. Batch 2 contained 10 control and 5 treatment
samples. Batch 3 contained 19 control and 18 treatment samples. After
filtering out ambiguously mapped genes and genes with low counts across
samples using the _edgeR_ package [CLS16], raw count data was voom-transformed
and differential expression analysis was performed using the _limma-voom_
workflow [Rit+15][Law+14]. 14708, 12230, and 14948 genes went into the
differential expression analysis for Batch 1, 2, and 3 respectively.

We note that in our two-sample differential expression analysis design for
this study, we do not aim to identify biologically meaningful insights, but
rather use a vehicle control versus muDX400 treatment as a way to showcase the
behavior of online FDR control methods compared to their offline counterparts.

We apply offline, pooled offline, and three Batch online approaches to three
real-world datasets (annotated as Batches 1, 2 and 3). As a benchmark, we see
that almost all of the FDR control methods are more conservative than the
uncorrected approach at an $\alpha$ of 0.05 at all time points (Table 1). The
offline StBH and pooled offline StBH approaches make more rejections compared
to their respective BH counterparts. This is expected because StBH is more
adaptive to the data at hand and has been shown empirically to yield higher
power [Sto02]. Zrnic et al. have shown theoretical guarantees that FDR will be
controlled across all batches and even future batches, assuming _p_ -values
are independent ([Zrn+21]. Although BatchBH and BatchStBH make fewer
rejections than their offline and pooled offline analogues, it is likely that
they are acting more conservatively to maintain overall FDR control. It is
worth noting that in this real-world example, BatchStBH rejecte more genes as
differentially expressed (DE) than offline and pooled BH at each time point
(Table 1).

Table 1: Cumulative number of rejections (proportion of total genes tested) in
real-world RNAseq data across batches. N is the total number of genes tested.
Blue color scales with the number of null hypothesis rejections with dark blue
for high number of rejections, light blue for low number of rejections.

The BatchBH procedure made the same rejections for DE genes 56% of the time as
its offline and pooled analogues across three accumulated datasets (Figure 2).
BatchStBH achieved a higher proportion of overlapping rejections of 66% across
the three datasets. Online approaches are able to detect a moderately high
percentage of the same DE genes as the offline approaches, but the online
approaches in this real-world data example appear more conservative.

Figure 2: Overlap in differentially expressed Genes Between Offline, Pooled
Offline, and Batch Approaches, $\alpha=0.05$

We next want to illustrate that naively composing offline FDR control
approaches is not appropriate, using Pooled BH between Batches 1 and 2 as an
example. In Batch 1, a set of 3469 genes were declared DE by Pooled BH (Table
1). In the pooled approach, we first performed the differential expression
analysis independently for Batch 1 and Batch 2, and then appended the
resulting p-values together and performed an FDR correction. Cumulatively
between Batches 1 and 2, 4600 genes were declared DE; however, only 3069 genes
from dataset 1 were declared DE. 400 genes “lost” their significance because
BH reorders all the _p-_ values. To further illustrate the change in historic
decisions in a pooled BH approach, we randomly sampled 50 genes declared DE
and 50 genes declared non-DE from Batch 1 100 times. Of the genes that were
declared DE in Batch 1, the number of genes that had decisions changed with a
pooled application ranged from 1 to 10 (average of 6). One of the random
samples is shown in Figure 3. When comparing the differential expression
decisions (e.g. null hypothesis rejections) for the sampled genes in Batch 1
to the decisions for the same sampled genes from Batch 1 but from the pooled
FDR correction across Batches 1 and 2 and across Batches 1-3, we observe how
the differential expression rejection decisions for the sampled genes from
Batch 1 have changed. It is likely that _p-_ values that are close to the
critical value as defined in the BH algorithm will have their decisions
changed. Potentially changing whether a gene is rejected as being DE each time
an additional batch of genes is tested is clearly infeasible from the
practical perspective of wanting to act in a timely manner on the hypothesis
test results (e.g. carrying out subsequent functional validation experiments).
It is important to emphasize that online methods do not change the historic
rejection decisions made when new datasets of _p_ -values arrive and are
pooled with existing datasets.

Figure 3: Schematic visualizing the pooled offline approach and how gene
differential expression results from Batch 1 are improperly changing
(highlighted in black boxes). Online methods would produce a visualization
akin to triplicating the t1 column.

##  4 Simulation Study

To explore the application of the online FDR control methods on a broader
scale, we performed a RNAseq simulation study using the compcodeR package
[Son14], which simulates read counts from a negative binomial distribution
(see Supplementary Appendix for more details). We used one of the standard
practice workflows for RNAseq differential expression analysis and applied the
voom method [Law+14] from the limma package (v3.50.0) [Rit+15].

In all our simulation experiments, 50 RNAseq datasets (referred to as
“batches”) each with 10,000 genes were simulated to represent a group of
batches arriving over time (which we henceforth refer to as a “family”). For
each family, we simulated a range of proportions of DE genes:

$\displaystyle\\{0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.2,0.3,0.4,0.5\\}$

We focus on two-group comparisons only. We specified 5 replicates for both the
control and experimental groups, a log-fold change simulation parameter of
1.5, a sequencing depth of $10^{7}$, minimum and maximum factors that are
multiplied with sequencing depth of $0.9$ and $1.1$ respectively, and a
proportion of upregulated genes set to 0.5, which were all mostly default
values in the _compcodeR_ package. We then ran _limma-voom_ as our
differential expression analysis method. This simulation setup was replicated
100 times independently for computational reasons.

For each family, we applied offline BH, offline StBH, pooled offline BH,
pooled offline StBH, and the online procedures BatchBH, BatchStBH, and
BatchPRDS. The FDP was calculated as the number of false positives divided by
the total number of genes declared DE. The average FDP across 100 trial runs
provides empirical approximation of the true FDR. Power was calculated as the
number of truly DE genes divided by the total number of genes declared DE and
averaged over 100 trial runs. For offline BH and StBH, this entailed applying
BH and StBH to each of the 50 datasets within a family, and then calculating
FDR and power. For pooled BH and StBH this entailed applying BH and StBH a
single time to all 50 datasets combined within a family. BatchBH, BatchStBH,
and BatchPRDS were applied as described in
https://dsrobertson.github.io/onlineFDR/articles/onlineFDR.html [Zrn+21].

Figure 4: FDR (A) and average power (B) versus proportion of DE genes ($\pi$)
comparing offline, pooled offline, and online Batch algorithms at batch size
of 10000. Number of batches set to 50 at a constant $\pi$. $\alpha$ set to
0.05. Simulated log fold-change set to 1.5. Shaded ribbons represent empirical
95% confidence bounds.

Comparing the FDR across a range of proportions of DE genes the offline and
pooled offline methods performed very similarly (Figure 4A). We observe that
when either the offline BH (as is commonly done) or pooled offline BH
approaches are performed on individual datasets, and the individual datasets
are naively composed as a family of datasets arriving over time, the standard
deviation of FDP is highly variable at low $\pi$ (e.g. $\pi\leq 0.1$),
exceeding the pre-specified significance threshold. At higher $\pi$, FDR
control is achieved by both offline BH and pooled offline BH. Offline StBH and
pooled offline StBH on average slightly inflate the FDR above the threshold
with a similarly high extent of variability at low $\pi$ as their BH
counterparts. BatchBH and BatchPRDS control FDR well below the significance
threshold throughout, with BatchPRDS exhibiting more conservatism. BatchPRDS
is also highly variable at very low $\pi$; however, it still controls FDR at
the nominal level as expected from its theoretical guarantees. In our
simulation, BatchStBH experiences a slight inflation of FDR at higher $\pi\geq
0.4$.

BatchBH and BatchStBH had close power to that of offline BH and pooled offline
approaches starting from $\pi=0.05$, which became closer as $\pi$ increased.
BatchPRDS was very conservative across the range of $\pi$ tested. In this
simulated scenario, all methods had low power at $\pi<0.1$. Figure 4 shows
that offline and pooled offline methods do not control the FDR in the online
framework where online methods do, while still maintaining practically
comparable levels of power.

##  5 The effect of ordering

Ordering of hypotheses can be crucial for the outcome of online testing, since
it can lead to either a loss or gain of power [JM18]. To increase power,
hypotheses can be ordered, using side information, such that those most likely
to be rejected are tested first. Hence, we sought to compare the performance
of offline and online methods in a simulated ordered setting.

We have simulated a ‘best-case’ scenario in which a scientist decides to order
RNAseq experiments so that batches with a higher proportion of DE genes are
tested first, where we assume that there is _a priori_ justification. We
represent this scenario by simulating the first 10 datasets (termed “Super
Batches”) as having a higher proportion $\pi$ of DE genes, followed by 40
datasets (termed “Regular Batches”) that have a lower proportion $\pi$ of DE
genes, such that the overall proportion $\pi$ of DE genes is one of
$\\{0.1,0.2,0.3,0.4,0.5\\}$:

Table 2: Ordered Setup of Number and Proportion of Super Batches and Regular Batches Family | Super Batch | Super Batch $\pi$ | Regular Batch | Regular Batch $\pi$ | Overall $\pi$
---|---|---|---|---|---
50 | 10 | 0.3 | 40 | 0.050 | 0.1
50 | 10 | 0.3 | 40 | 0.175 | 0.2
50 | 10 | 0.5 | 40 | 0.250 | 0.3
50 | 10 | 0.5 | 40 | 0.375 | 0.4
50 | 10 | 0.5 | 40 | 0.500 | 0.5

In this setting, offline StBH and pooled offline StBH do not control FDR at
the specified level, across different $\pi$ (Figure 5A). The ordering of
hypothesis tests has benefited the BH-based algorithms where both offline BH
and pooled offline BH now control the FDR. BatchBH and BatchPRDS still control
the FDR at various $\pi$; however, BatchStBH inflates FDR to comparable levels
as $\pi>0.2$. We return to this issue in the Discussion. The ordering has also
benefited all approaches in terms of power. The power achieved by the online
algorithms is very close to that achieved by their offline and pooled offline
counterparts across all $\pi$.

Figure 5: FDR (A) and average power (B) versus probability of non-null
hypotheses $\pi$ comparing offline, pooled offline, and online Batch
algorithms at batch size of 10000 in an ordered setting such that datasets
with a higher proportion of differentially expressed genes are tested first.
Number of batches set to 50. $\alpha$ set to 0.05. Simulated log fold-change
set to 1.5. Shaded ribbons represent empirical 95% confidence bounds

##  6 Discussion

In this paper, we demonstrate that the online FDR control algorithms provide a
principled way to guarantee control of FDR at a nominally specified level in
an online paradigm, e.g. RNAseq datasets arriving sequentially over time, and
at every time step, a batch of decisions is made via the BH or StBH algorithm.
We show that offline approaches have variable control of FDR of related RNAseq
experiments over time in certain settings simulations, and that a pooled
application may improperly change previous differential expression rejection
decisions, which can be construed as one form of p-hacking [Hea+15].
Furthermore, in some simulation scenarios, we observe empirically that online
approaches have comparable power to offline approaches.

Using three real RNAseq datasets that we considered as a “family” of datasets
arriving over time, we applied the online FDR control methodologies. Although
our datasets were limited in that we could not assess “ground truth”, we see
that the online methodologies are broadly more conservative than their offline
and pooled offline counterparts, which is likely driven by the theoretical
long-term FDR control conferred. Notably, as the number of datasets increase,
the online methods are able to recapitulate a higher proportion of genes as
detected by more conventional offline methods, which suggests that while the
online approaches may have lower statistical power in the short-term, in the
long-term, they may have practically comparable power to offline methods. We
selected the datasets that were similar enough in terms of identical time
points, treatment, and mouse model to reduce measurable confounding to the
best of our ability. But we also remark that the primary motivation behind
this analysis was not biologically driven, in that questions such as whether
similar gene biology was detected across the datasets were considered to be
beyond the scope of this paper. Generally, similar RNAseq datasets added to a
growing data repository such as the Gene Expression Omnibus (GEO) are
potentially indicated use cases for the application of online FDR algorithms.

In our unordered simulation setting, we observed that the online algorithms
successfully controlled FDR across a wide range of $\pi$. Is is notable that
particularly at low proportions of differentially expressed genes, the offline
approaches did not maintain control of FDR. The range of $\pi$ that we test
may be relevant for the range of number of DE genes that may be detected in
practice [Cor+17]. BatchStBH had slight inflation of FDR, but that was at a
very high $\pi$ of 0.5 in which the distribution of DE and non-DE genes in a
RNAseq dataset is symmetrical, which is likely very rare in practice. The
slight inflation of FDR could have also resulted from the subtle dependency
introduced in our simulation methods. The _compcodeR_ package draws $\mu$ and
$\phi$ estimates for read count generation from the same underlying real-world
data distribution, so it is possible that, although we have simulated the
datasets independently to the best of our ability, there still exists some
degree of dependency amongst the _p_ -values. BatchPRDS was more conservative
than BatchBH and BatchStBH; it guarantees FDR control when _p-_ values are
dependent within each batch but at the cost of a substantially lower power
than any of the other algorithms. In real RNAseq datasets, it is likely that
DE gene _p-_ values are dependent [Sub+05]. Future work should investigate the
performance of Batch online FDR control methods within varying dependency
structures. Subramanian et al. also proposed evaluating sequencing data at the
level of gene sets since single-gene analysis may miss important effects on
pathways. Thus, one other potential future research direction is to apply
online control methods to gene set analysis _p-_ values.

We also simulated an ordered setting, since in practice, scientists may order
batches of hypotheses in such a way that seeks to optimizes their power and
reflects their _a priori_ knowledge. We show in our simulations that ordering
hypotheses increases the power of online methods, but appears to inflate the
FDR of BatchStBH at higher $\pi$. We note a limitation of our simulation where
in practical scenarios, it may be a more nuanced decision between ordering
batches of RNAseq datasets or the hypotheses (i.e. gene _p_ -values)
themselves. In our application of online methods to our real-world data, we
did not explore _a priori_ ordering of datasets as biological rationale was
not the emphasis. However, it is possible that if we were to test Batch 2
after both Batches 1 and 3 (as it appeared inherently less powerful) for some
valid biological justification (e.g. this was a technical pilot of the RNAseq
experimental setup itself where high signal was not expected), we may achieve
better overall power across the three datasets together in an online
framework. Future studies could investigate ordering retrospectively (e.g. in
a database such as GEO) or prospectively (e.g. making decisions how to order
datasets arriving over time). We reiterate that the purpose of our simulation
was to represent how much more powerful online methodologies could be if we
happened to have inherently more powerful datasets ordered first. In addition,
in reality, we may not expect such a high proportion of DE genes, nevertheless
we caution the use of BatchStBH in an ordered setting as it might “greedily”
increase statistical power at the expense of FDR.

Additional future directions of exploring FDR control with RNAseq data can
explore and benchmark how varying batch size, the number of samples, and using
different gene filtering methods or differential expression methods would
affect performance of online algorithms. Using other FDR control methods
(besides BH and Storey-BH) interpolated in the online manner as we have
described is also possible.

## Acknowledgements

We would like to thank Eric Muise, Wendy Blumenschein, Doug Linn for their
permission to use the real-world datasets as well as their comments and
feedback. We would also like to thank Ola Olow and Uwe Mueller for their
detailed comments and suggestions to improve this manuscript.

## Conflicts of Interest

We have no conflicts of interest to declare.

## Funding

This work was supported by the UK Medical Research Council [MC_UU_00002/14 to
D.S.R.], the Biometrika Trust [to D.S.R.] and the NIHR Cambridge Biomedical
Research Centre [BRC1215-20014 to D.S.R.]. The views expressed in this
publication are those of the authors and not necessarily those of the NHS, the
National Institute for Health Research or the Department of Health and Social
Care (DHSC). For the purpose of open access, the author has applied a Creative
Commons Attribution (CC BY) licence to any Author Accepted Manuscript version
arising.

## Data availability statement

Code used for the simulation studies can be found at
https://github.com/latlio/onlinefdr_rnaseq_simulation. The real-world RNAseq
data underlying this article will be shared on reasonable request by the
corresponding author.

## References

  * [BH95] Yoav Benjamini and Yosef Hochberg  “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing” Publisher: [Royal Statistical Society, Wiley]  In _Journal of the Royal Statistical Society. Series B (Methodological)_ 57.1, 1995, pp. 289–300  URL: http://www.jstor.org/stable/2346101
  * [YD01] Yoav Benjamini and Daniel Yekutieli  “The control of the false discovery rate in multiple testing under dependency”  In _The Annals of Statistics_ 29.4, 2001, pp. 1165–1188  DOI: 10.1214/aos/1013699998
  * [EDL02] Ron Edgar, Michael Domrachev and Alex E. Lash  “Gene Expression Omnibus: NCBI gene expression and hybridization array data repository.”  In _Nucleic acids research_ 30.1, 2002, pp. 207–210  DOI: 10.1093/nar/30.1.207
  * [Sto02] John D. Storey  “A Direct Approach to False Discovery Rates” Publisher: [Royal Statistical Society, Wiley]  In _Journal of the Royal Statistical Society. Series B (Statistical Methodology)_ 64.3, 2002, pp. 479–498  URL: http://www.jstor.org/stable/3088784
  * [Sub+05] Aravind Subramanian et al.  “Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles.”  In _Proceedings of the National Academy of Sciences of the United States of America_ 102.43, 2005, pp. 15545–15550  DOI: 10.1073/pnas.0506580102
  * [CRS10] James J. Chen, Paula K. Roberson and Michael J. Schell  “The false discovery rate: a key concept in large-scale genetic studies.” Place: United States  In _Cancer control : journal of the Moffitt Cancer Center_ 17.1, 2010, pp. 58–62  DOI: 10.1177/107327481001700108
  * [Law+14] Charity W Law, Yunshun Chen, Wei Shi and Gordon K Smyth  “voom: Precision weights unlock linear model analysis tools for RNA-seq read counts” Publisher: BioMed Central  In _Genome biology_ 15.2, 2014, pp. R29–R29  DOI: 10.1186/gb-2014-15-2-r29
  * [LHA14] Michael I. Love, Wolfgang Huber and Simon Anders  “Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2”  In _Genome Biology_ 15.12, 2014, pp. 550  DOI: 10.1186/s13059-014-0550-8
  * [Son14] Charlotte Soneson  “compcodeR–an R package for benchmarking differential expression methods for RNA-seq data.” Place: England  In _Bioinformatics (Oxford, England)_ 30.17, 2014, pp. 2517–2518  DOI: 10.1093/bioinformatics/btu324
  * [Hea+15] Megan L. Head et al.  “The extent and consequences of p-hacking in science.”  In _PLoS biology_ 13.3, 2015, pp. e1002106  DOI: 10.1371/journal.pbio.1002106
  * [Rit+15] Matthew E Ritchie et al.  “limma powers differential expression analyses for RNA-sequencing and microarray studies” Edition: 2015/01/20 Publisher: Oxford University Press  In _Nucleic acids research_ 43.7, 2015, pp. e47–e47  DOI: 10.1093/nar/gkv007
  * [CLS16] Yunshun Chen, Aaron T L Lun and Gordon K Smyth  “From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline” Publisher: F1000Research  In _F1000Research_ 5, 2016, pp. 1438–1438  DOI: 10.12688/f1000research.8987.2
  * [Cor+17] Susan M. Corley et al.  “Differentially expressed genes from RNA-Seq and functional enrichment results are affected by the choice of single-end versus paired-end reads and stranded versus non-stranded protocols”  In _BMC Genomics_ 18.1, 2017, pp. 399  DOI: 10.1186/s12864-017-3797-0
  * [JM18] Adel Javanmard and Andrea Montanari  “Online rules for control of false discovery rate and false discovery exceedance” ISBN: 0090-5364 Publisher: JSTOR  In _The Annals of statistics_ 46.2, 2018, pp. 526–554
  * [Kor+19] Keegan Korthauer et al.  “A practical guide to methods controlling false discoveries in computational biology”  In _Genome Biology_ 20.1, 2019, pp. 118  DOI: 10.1186/s13059-019-1716-1
  * [Rob+19] David S Robertson, Jan Wildenhain, Adel Javanmard and Natasha A Karp  “onlineFDR: an R package to control the false discovery rate for growing data repositories” Publisher: Oxford University Press  In _Bioinformatics (Oxford, England)_ 35.20, 2019, pp. 4196–4199  DOI: 10.1093/bioinformatics/btz191
  * [Zrn+21] Tijana Zrnic, Daniel L Jiang, Aaditya Ramdas and Michael I Jordan  “The Power of Batching in Multiple Hypothesis Testing”  In _arXiv preprint arXiv:1910.04968_ , 2021