\documentclass{article}
\usepackage[margin=3cm]{geometry}
\usepackage{graphicx}
\usepackage{float}
\usepackage{amsmath} 
\usepackage{textgreek}

% Hypperef settings.
\usepackage[hidelinks]{hyperref}
\urlstyle{same}

% Line space settings.
\usepackage{setspace}
\spacing{1.5}

% Biblatex customizations to hack into GR style.
\usepackage[terseinits=true,style=authoryear,firstinits=true,maxnames=10,maxcitenames=2,mincitenames=1,minbibnames=10,natbib,dashed=false]{biblatex}
\renewbibmacro{in:}{}
\DeclareFieldFormat[article, inbook, incollection, inproceedings, misc, thesis, unpublished]
{title}{#1}
\DeclareFieldFormat[article]
{pages}{#1}
\DeclareFieldFormat[article]
{volume}{\textbf{#1}}
\DeclareFieldFormat[article, inbook, incollection, inproceedings, misc, thesis, unpublished]
{number}{\mkbibparens{#1}} 
\renewbibmacro*{volume+number+eid}{
    \iffieldundef{volume}
    {}{\printfield{volume}\printfield{number}\printunit{\addcolon}}
}
\renewbibmacro*{date+extrayear}{
  \printunit{\addperiod\space}\printfield{year}
}
\renewcommand*{\revsdnamepunct}{} 
\renewcommand*{\finalnamedelim}{\addcomma\space}
\DeclareNameAlias{sortname}{last-first}

\bibliography{refnorm}

% References to supplementary files.
\newcommand{\suppfigschem}{1}
\newcommand{\suppfignorm}{2}
\newcommand{\suppfigtotals}{3}
\newcommand{\suppfigorder}{4}
\newcommand{\suppfigcell}{5}
\newcommand{\suppfigbiophys}{6}
\newcommand{\suppfignoise}{7}
\newcommand{\suppfignoisier}{8}
\newcommand{\suppfigsizevar}{9}
\newcommand{\suppfiglostsim}{10}
\newcommand{\suppfigspikeconc}{11}
\newcommand{\suppfigspikeave}{12}
\newcommand{\suppfigindex}{13}

\newcommand{\suppsecmath}{1}
\newcommand{\suppsecnoise}{2.1}
\newcommand{\suppsecnoisier}{2.2}
\newcommand{\suppsecindex}{3}
\newcommand{\suppsecdata}{4}
\newcommand{\suppsecsim}{5}

\newcommand{\supptabstats}{1}

% Revision commands. 
\usepackage{color}
\newcommand{\revised}[1]{\textcolor{red}{#1}}

% Convenience function for writing variance in math mode.
\newcommand\variance{\mbox{var}}

% Title information.
\title{Assessing the reliability of spike-in normalization for analyses of single-cell RNA sequencing data}
\author{Aaron T. L. Lun$^1$, Fernando J. Calero-Nieto$^2$, Liora Haim-Vilmovsky$^{3,4}$, \\ Berthold G\"ottgens$^2$, John C. Marioni$^{1,3,4,*}$}
\date{
    \begin{minipage}{0.9\textwidth}
        \begin{flushleft} 
            \begin{small}
                $^1$Cancer Research UK Cambridge Institute, University of Cambridge, Li Ka Shing Centre, Robinson Way, Cambridge CB2 0RE, United Kingdom \\
                $^2$Wellcome Trust and MRC Cambridge Stem Cell Institute, University of Cambridge, Wellcome Trust/MRC Building, Hills Road, Cambridge CB2 0XY, United Kingdom \\
                $^3$EMBL European Bioinformatics Institute, Wellcome Genome Campus, Hinxton, Cambridge CB10 1SD, United Kingdom \\
                $^4$Wellcome Trust Sanger Institute, Wellcome Genome Campus, Hinxton, Cambridge CB10 1SA, United Kingdom \\
                $^*$Corresponding author (email: marioni@ebi.ac.uk)
            \end{small}
        \end{flushleft}
    \end{minipage}\\[0.2in]
    \today{}
}

\begin{document}

\maketitle

\begin{quote}
\textbf{Abstract:}
By profiling the transcriptomes of individual cells, single-cell RNA sequencing provides unparalleled resolution to study cellular heterogeneity.
However, this comes at the cost of high technical noise, including cell-specific biases in capture efficiency and library generation.
One strategy for removing these biases is to add a constant amount of spike-in RNA to each cell, and to scale the observed expression values so that the coverage of spike-in \revised{transcripts} is constant across cells.
This approach has previously been criticized as its accuracy depends on the precise addition of spike-in RNA to each sample.
Here, we perform mixture experiments using two different sets of spike-in RNA to quantify the variance in the amount of spike-in RNA added to each well in a plate-based protocol.
We also obtain an upper bound on the variance due to differences in behaviour between the two spike-in sets.
We demonstrate that both factors are small contributors to the total technical variance and have only minor effects on downstream analyses such as detection of highly variable genes and clustering.
Our results suggest that \revised{scaling} normalization \revised{using spike-in transcripts} is reliable enough for routine use in single-cell RNA sequencing data analyses.
\end{quote}

\section*{Introduction}
Single-cell RNA sequencing (scRNA-seq) is a powerful technique for studying transcriptional activity in individual cells.
Briefly, RNA is isolated from single cells, reverse transcribed into cDNA and sequenced using massively parallel sequencing technologies \autocite{shapiro2013singlecell}.
This can be performed using microfluidics platforms like the Fluidigm C1 \autocite{pollen2014lowcoverage}; 
    with protocols such as Smart-seq2 \autocite{picelli2014full} that use microtiter plates;
    or with droplet-based technologies \autocite{klein2015droplet,macosko2015highly} that can profile thousands of cells.
Gene expression is quantified by mapping read sequences to a reference genome and counting the number of reads mapped to each annotated gene.
To avoid amplification biases, individual transcript molecules can also be tagged with unique molecular identifiers (UMIs) \autocite{islam2014quantitative}, such that sequencing to saturation and counting UMIs will yield the number of transcripts of each gene in a cell.
Regardless of whether reads or UMIs are used, not all transcript molecules will be captured and sequenced due to cell-specific inefficiencies in reverse transcription \autocite{stegle2015computational}.
The presence of these cell-specific biases compromises the direct use of the read/UMI count as a quantitative measure of gene expression.
Normalization is required to remove these biases before the gene counts can be meaningfully compared between cells in downstream analyses.

A common normalization strategy for RNA-seq data uses a set of genes that have constant expression across cells.
This set can consist of pre-defined ``house-keeping'' genes, or it can be empirically defined under the assumption that most genes are not differentially expressed (DE) between cells \autocite{lun2016pooling,anders2010differential,robinson2010tmm}.
Any systematic differences in expression between cells for this non-DE set of genes must, therefore, be technical in origin, e.g., due to differences in library size or composition bias \autocite{robinson2010tmm}.
Counts are scaled to eliminate these differences, yielding normalized expression values for downstream analyses.
This gene-based approach works well for bulk sequencing experiments where the population-wide gene expression profile is stable.
However, it may not be suitable for single-cell experiments where strong biological heterogeneity complicates the identification of a reliable non-DE set. 
For example, house-keeping genes may be turned on or off by transcriptional bursting, while processes like the cell cycle may trigger large-scale changes in the expression profile that preclude a non-DE majority.

An alternative normalization approach is to use spike-in RNA for which the identity and quantity of all transcripts is known \autocite{stegle2015computational,bacher2016design}.
The same amount of spike-in RNA is added to each cell's lysate, and the spike-in transcripts are processed in parallel with their endogenous counterparts to generate a sequencing library.
This yields a set of read (or UMI) counts for both endogenous and spike-in transcripts in each cell.
Normalization is performed by scaling the counts for each cell such that the counts for the spike-in genes are, on average, the same between cells \autocite{katayama2013samstrt}.
The central assumptions of this approach are that (i) the same amount of spike-in RNA is added to each cell, and (ii) the spike-in and endogenous transcripts are similarly affected by cell-to-cell fluctuations in capture efficiency.
Under these assumptions, any differences in the coverage of the spike-in transcripts between cells must be artifactual in origin and should be removed by scaling.
One particular advantage of this strategy is that it does not make any assumptions about the endogenous expression profile, unlike the non-DE approach described above.
This means that spike-in normalization can be applied in situations where large-scale changes in expression (e.g., related to changes in total RNA content, or involving highly heterogeneous populations containing many cell types) are expected and of interest \autocite{lun2016stepbystep,nestorowa2016single}.

There are two common criticisms of spike-in normalization that challenge the validity of its central assumptions.
The first is that the same quantity of spike-in RNA may not be consistently added to each sample \autocite{robinson2010tmm}, and the second is that synthetic spike-in transcripts may not behave in the same manner as endogenous transcripts \autocite{grun2015design} (i.e., \revised{the two sets of transcripts} have unequal capture efficiencies, caused by differences in their biophysical properties).
Any differences in spike-in quantity or behaviour across cells will compromise the accuracy of spike-in normalization \autocite{risso2014normalization}.
In some cases, it may also be difficult to gauge how much spike-in RNA should be added, especially if the quantity of endogenous RNA per cell is unknown, resulting in insufficient spike-in coverage for normalization.
These criticisms may contribute to the limited use of this normalization strategy in the scRNA-seq literature \autocite{bacher2016design}.
However, if one were to dismiss the use of spike-in normalization, there would be no general alternative for removing cell-specific biases in scRNA-seq data sets where a non-DE majority of genes cannot be assumed.
Thus, it is of particular interest whether or not the aforementioned criticisms of spike-in normalization are relevant to real scRNA-seq experiments.
To our knowledge, this has yet to be rigorously studied.

% The advantage of spike-ins is that it can capture cell size and doesn't require a non-DE majority.
% However, is the cell size interesting, or should it be normalized out (in which case accurate spike-in addition wouldn't be required at all)?
% It's most obviously useful for cell cycle considerations -- possibly also cancer and in some aspects of differentiation.
% The Islam non-UMI dataset is one example where there seems to be a 20-50-fold increase in counts in the MEFs vs mESCs.
% 
% I guess spike-in normalization is more technically correct, even if you end up getting every gene as being DE (technically true in terms of molecules, but possibly irrelevant).
% Rankings definitely change - spike-in normalization would favour genes that are upregulated on top of an increase in RNA content (as these would be hugely significant).
% Standard normalization would have no preference as you seek the balance between the two extremes.
% On the other hand, any non-DE genes with no change in the molecule number will be called as significant in standard normalization, which is definitively wrong.
% You could have your cake and eat it by separating the DE list into up/down changes for easier examination of the up/down changes.

\section*{Results}

\subsection*{Overview}
In this paper, we conduct a series of experiments to estimate the reliability of spike-in normalization in single-cell transcriptome studies employing plate-based protocols.
We use mixtures of two distinct spike-in RNA sets to quantify the variance of the added spike-in volume across cells, 
and show that it is quantitatively negligible in real experiments across a range of conditions. 
We also obtain an upper bound on the cell-to-cell variability in the differences in behaviour (i.e., the fold-changes in the capture efficiencies) between the two spike-in sets.
Simulations indicate that both factors have only minor effects on the results of downstream analyses such as detection of DE and highly variable genes.
These results suggest that spike-ins can be safely used for routine normalization of scRNA-seq data.

\revised{We emphasize that we are only interested in the performance of spike-in RNA for scaling normalization.
This involves the calculation of cell-specific scaling factors to remove relative biases between cells.
We are not investigating the performance of spike-in RNA for the absolute quantification of endogenous transcript molecules \autocite{svensson2017power},
which would require estimation of the absolute bias in each cell.
We are also not studying the use of spike-ins for batch correction \autocite{tung2016batch}, which would require modelling of gene-specific batch effects beyond simple cell-specific scaling.
Both of these tasks are separate to scaling normalization and will not be addressed here.
}

\subsection*{Description of the mixture experiments}
We aimed to assess the variability in the added spike-in quantity across cells.
To do so, we performed mixture experiments using two distinct spike-in sets (Figure~\ref{fig:expdesign}) -- the External RNA Controls Consortium (ERCC) set and the Spike-in RNA Variants (SIRV) set.
An equal volume of each spike-in set was added separately to all wells of a 96-well microtiter plate.
Each well contained a single lysed mouse cell -- a mouse 416B myeloid progenitor cell or trophoblast stem cell (TSC) -- thus mimicking real experimental conditions.
The resulting pool of endogenous/spike-in RNA in each well was used to generate a cDNA library, using a modified version of the Smart-seq2 protocol (see Methods).
This process was repeated for all wells and high-throughput sequencing was performed on all libraries.

\begin{figure}[tbp]
\begin{center}
\includegraphics[width=0.8\textwidth]{pics/plate_setup.pdf}
\end{center}
\caption{Schematic of the experimental design to assess the variability of spike-in addition in a plate-based scRNA-seq protocol.
(a) A cell is sorted into each well of a plate and lysed.
For one set of wells, an equal volume of each spike-in set is added separately, along with the reverse transcription (RT) reagents.
For another set of wells, an equal volume of a pooled mixture of the two spike-ins is added into each well (done twice to keep the protocol consistent).
Reverse transcription, PCR amplification, library generation and sequencing were then performed.
(b) The log$_2$-ratio between the total counts of the two spike-in sets was computed for each well.
The variance of the log-ratio was estimated from all wells with separate addition of spike-ins, and from wells with addition of the premixed pool.
The difference between these two estimates represents the variance attributable to volume addition.
}
\label{fig:expdesign}
\end{figure}

For each library, reads were mapped to the genome and assigned to genes to quantify expression.
The total count was computed across all transcripts of each spike-in set in each well.
The log$_2$-ratio of the totals between the two sets was computed for each well, and the variance of this log-ratio was computed \textit{across} wells.
Any variability in spike-in volume addition should manifest as an increase \revised{in} the variability of the log-ratio, given that the spike-in sets were added independently to each well.

We also repeated the experiment by adding volumes of ``premixed'' spike-in solution where the two spike-in sets had been pooled at a 1:1 ratio.
This ensures that there is no well-to-well variability in the relative quantities of RNA from the two spike-in sets.
The variance of the log-ratio across these premixed-addition wells provides a baseline level of variability in the protocol (e.g., due to sequencing noise).
The variance of volume addition was then estimated as the difference in the variance estimates from the premixed-addition wells and from the wells with separate addition of spike-ins.

We stress that the use of two different spike-in sets in each well is critical to this experiment.
Any well-specific biases should cancel out when the log-ratio is computed between sets in the same well.
This allows the contribution of the variance of volume addition to be quantified separately from other factors such as the variability of capture efficiency and sequencing depth across wells.

We performed both the premixed and separate-addition experiments on the same plate to avoid plate effects \autocite{hicks2015widespread,tung2016batch}.
For the separate-addition experiment, we also reversed the order of addition of the two spike-in sets to determine if this affected the variance estimate.
Finally, we generated data from replicate plates to ensure our results were reproducible.
This was done in a range of conditions, i.e., using different cell types, by different operators and with sequencing at different locations.

We used a protocol based on microtiter plates rather than microfluidics as it is easier to customise the spike-in addition step in the former.
Our experimental design requires two separate additions of spike-in RNA to each reaction (see Methods).
This is not straightforward to achieve on, say, the Fluidigm C1 chip where the added volume for each reagent depends on the design on the reaction chamber.
%Plate-based protocols also tend to be cheaper and are compatible with existing laboratory techniques such as indexed fluorescently-activated cell sorting (FACS).
Our focus on data from plate-based protocols reflects their widespread use in single-cell studies \autocite{segerstople2016single,islam2011characterization,wilson2015combined,scialdone2016resolving}.
Obviously, the procedure we describe here can be adapted to any protocol where the spike-in addition can be easily modified, e.g., plate-based CEL-seq \autocite{hashimshony2016celseq2} or STRT-seq \autocite{islam2011characterization}.

\subsection*{Mathematical framework for variance decomposition}
Denote the log$_2$-transformed total read count for well $i$ and spike-in set $s$ as
\[
T_{is} = \log_2 \left[ L_i l_s V_{is} R_{is} \sum_{t_s} r_{t_s} c_{t_s} \right] + \varepsilon_{is}
\]
where the sum is taken over all unique transcripts $t_s$ in $s$.
The other terms are defined as follows:
\begin{itemize}
    \item $c_{t_s}$, a constant specifying the concentration (in terms of transcripts per unit of volume) of $t_s$.
    \item $r_{t_s}$, a constant specifying the optimal transcript molecule-to-cDNA fragment capture rate for $t_s$.
    \item $R_{is}$, a random variable representing the average capture efficiency in $i$ for all transcripts in $s$.
    \item $V_{is}$, a random variable representing the volume of solution of $s$ added to $i$.
    \item $L_i$, a random variable representing the baseline cDNA fragment-to-read conversion rate for $i$.
    \item $l_s$, a constant that scales $L_i$ depending on the ``sequenceability'' of transcripts in $s$.
\end{itemize}
The product of all of these terms defines the expected number of reads for each $t_s$ in well $i$, and the sum of the products across all $t_s$ is the expected total count of set $s$ in $i$.
In addition, $\varepsilon_{is}$ represents the effect of sequencing noise on the log-total count, where $E(\varepsilon_{is})=0$ and $\variance(\varepsilon_{is})= \sigma^2_{lib(s)}$.

We assume that $R_{is}$, $V_{is}$ and $\varepsilon_{is}$ are mutually independent of each other, as they describe separate steps in the protocol.
We also assume that $V_{i1}$ and $V_{i2}$ are independent for sets $s=1$ and 2, as each spike-in set is added separately to each well.
Similarly, $\varepsilon_{i1}$ and $\varepsilon_{i2}$ are assumed to be independent as sequencing noise for each transcript should be unaffected by that of other transcripts.
(However, $R_{i1}$ and $R_{i2}$ are not independent due to well-specific factors affecting capture efficiency for all transcripts).
Further details on these variables are provided in Section~\suppsecmath{} of the Supplemental Materials.

Let $s=1$ represent the ERCC spike-in set and $s=2$ represent the SIRV spike-in set.
The log$_2$-total count across all spike-in transcripts in the ERCC and SIRV set is $T_{i1}$ and $T_{i2}$, respectively.
In the experiment where each spike-in set is added separately to each well, we denote the log$_2$-ratio of the total counts between the two sets as $\theta_i = T_{i1} - T_{i2}$ for well $i$.
This can also be written as
\[
    \theta_i = \log_2(V_{i1}) + \varepsilon_{i1} - \log_2(V_{i2}) - \varepsilon_{i2} + F_i + \log_2\left[ \frac{l_1 \sum_{t_1} r_{t_1} c_{t_1}}{l_2\sum_{t_2} r_{t_2} c_{t_2}} \right]
\]
where $F_i = \log_2(R_{i1}/R_{i2})$ and represents the log-fold change in the average capture efficiency between the two sets (i.e., the difference in behaviour of the transcripts).
Computing the variance of $\theta_i$ yields
\[
\variance(\theta_i) = 2 \sigma^2_{vol} + \sigma^2_{lib(1)} + \sigma^2_{lib(2)} + \variance(F_i)
\]
where $\sigma^2_{vol}$ is the variance of both $\log_2(V_{i1})$ and $\log_2(V_{i2})$.
The volume addition procedure is the same for each spike-in set, so $V_{i1}$ and $V_{i2}$ should have the same distribution.
We consider the variance of $F_i$ because $R_{i1}$ and $R_{i2}$ are not independent (due to well-specific factors, as previously mentioned).

In the experiment where the spike-in sets are premixed before addition, $V_{i1}=aV_{i2}$ for some constant $a$ representing the proportions in which the two sets are mixed.
(This should be close to unity.)
If the same premixed solution is added to each well, the relative volume of ERCC spike-ins to SIRV spike-ins must be constant for all wells.
This means that the log$_2$-ratio for the premixed experiment is 
\[
    \theta^*_i = \log_2(a) + \varepsilon_{i1} - \varepsilon_{i2} + F_i + \log_2\left[ \frac{l_1 \sum_{t_1} r_{t_1} c_{t_1}}{l_2\sum_{t_2} r_{t_2} c_{t_2}} \right] \;.
\]
As $a$ is constant for all $i$, the variance of $\theta^*_i$ becomes
\[
\variance(\theta^*_i) = \sigma^2_{lib(1)} + \sigma^2_{lib(2)} + \variance(F_i) \;.
\]
This represents the technical variance attributable to the rest of the scRNA-seq protocol.
To obtain an estimate of the variance of the volume addition step, simple arithmetic yields
\[
\sigma^2_{vol} = \frac{\variance(\theta_i) - \variance(\theta^*_i)}{2} \;.
\]
It should be stressed that this variance estimate is relevant to all experiments using the same protocol for spike-in addition, even if the identity or concentration of the spike-in set is different.

\revised{Generally, scaling normalization of RNA-seq data is performed by dividing all counts in each library by a library-specific constant, known as the ``size factor''.
For spike-in normalization, the size factor for cell $i$ is directly proportional to the sum of counts for the spike-in transcripts, i.e., $2^{T_{is}}$.
This reflects the fact that spike-in normalization aims to eliminate systematic differences in the coverage of spike-in set $s$ between cells, thus correcting for well/cell-specific technical biases.
(We assume each well contains a cell and will use ``cell'' and ``well'' interchangeably in the following text.)
Any variance due to volume addition ($\sigma^2_{vol}$) or technical noise ($\sigma^2_{lib(s)}$) will reduce the precision of $T_{is}$ and of the size factor estimates (Supplemental Figure~\suppfigschem{}), thus reducing the effectiveness of spike-in normalization.
}

\subsection*{Estimating the variance of volume addition}
Using our mathematical framework, we estimated the variance components based on the data from our mixture experiments.
We observed that the log-ratios $\theta_i$ and $\theta^*_i$ computed from each plate were roughly normally distributed (Supplemental Figure~\suppfignorm{}).
Thus, we fitted a linear model to each set of log-ratios and used the residual variance of the fit as our estimate of $\variance(\theta_i)$ or $\variance(\theta^*_i)$.
Linear models are particularly useful as they allow blocking on additional structure in the experimental design.
We used a one-way layout to account for shifts in the mean \revised{log-ratio} due to addition order or oncogene induction (see Methods).
The value of $T_{is}$ was also similar between wells with premixed or separate addition of spike-ins, which simplifies the calculation of $\sigma^2_{vol}$ (see Supplemental Figure~\suppfigtotals{}, Section~\suppsecmath{} of the Supplemental Materials for details).
Finally, the order of spike-in addition did not significantly affect the variance estimates for the separate-addition wells in most plates (Supplemental Figure~\suppfigorder{}).

Our results indicate that $\sigma^2_{vol}$ is consistently smaller than $\variance(\theta^*_i)$, i.e., the variance in the rest of the protocol (Figure~\ref{fig:varestimates}a).
Indeed, no significant difference was detected between the estimated $\variance(\theta_i)$ and $\variance(\theta^*_i)$ of each plate.
This indicates that variability of spike-in volume addition is a minor contributor to the technical variability of the spike-in counts.
To put these estimates into context, consider that the variance of the log-size factors $T_{is}$ across cells is at least one order of magnitude larger than $\sigma^2_{vol}$ (Figure~\ref{fig:varestimates}b).
\revised{This indicates that the error in the size factors due to variable volume addition is negligible relative to the amount of scaling that is performed to account for differences in sequencing depth and capture efficiency across wells, i.e., $\variance{(\log_2 L_{i})}$ and $\variance{(\log_2 R_{is})}$.}
We also computed the variance of the log$_2$-ratio of total counts for the mouse genes against one of the spike-in sets.
This represents the biological fluctuations in total RNA content across cells and was, again, at least an order of magnitude larger than $\sigma^2_{vol}$ (Supplemental Figure~\suppfigcell{}).
\revised{These results show that the variance of volume addition is small compared to other technical and biological sources of variability in a scRNA-seq experiment, and thus is unlikely to have a major effect on spike-in normalization.}

\begin{figure}[btp]
    \begin{center}
        \includegraphics[width=\textwidth,trim=0mm 15mm 0mm 0mm,clip]{../real/pics/variance_exp.pdf}
    \end{center}
    \caption{Variance estimates of (a) the log$_2$-ratio between the ERCC and SIRV total counts across wells, or (b) the log$_2$-size factors computed from those totals.
        For the separate/premixed experiments, each estimate is the residual variance of a linear model fitted to the log-ratios across the corresponding wells on each plate.
        The variance of volume addition is $\sigma^2_{vol}$ in our framework.
        For the log-size factors, each estimate is the residual variance of a linear model fitted to all cells on each plate.
        Results are shown for experiments with 416B cells or TSCs, with two replicate plates for each cell type.
        Error bars represent the standard errors of the estimates, assuming log-values are normally distributed.
        Numbers represent the residual degrees of freedom used for each estimate -- for (b), this was the same for each spike-in set.
        Differences between the separate-addition and premixed estimates for each batch were assessed using a one-sided F-test, yielding $p$-values of 0.28, 1.00, 1.00 and 0.06 from left to right.
    }
    \label{fig:varestimates}
\end{figure}

\subsection*{Estimating the variance of differential behaviour}
The variance of $F_i$ is also relevant as it determines the effect of differences in behaviour between distinct sets of transcripts.
Even when the average capture efficiency differs between sets, spike-in normalization is still appropriate \textit{provided that the fold change in efficiency is the same in all wells}.
Consider a situation where there is a consistent increase in efficiency in the spike-in set relative to endogenous transcripts \autocite{svensson2017power}.
This scales up the counts for the spike-in transcripts in all \revised{wells} by the same factor, which ultimately cancels out between \revised{wells}, i.e., the log-fold changes of endogenous or spike-in transcripts between \revised{wells} are unaffected.
However, if the fold change in efficiency varies across wells, the accuracy of spike-in normalization is compromised.
This is because specific changes in efficiency for the spike-in transcripts are confounded with general changes in efficiency for all transcripts in the well.
Differences in the coverage of spike-in transcripts may not represent technical biases affecting other transcripts, precluding their use for normalizing all counts.

The variance of $F_i$ quantifies the extent to which spike-in normalization is affected by well-to-well differences in efficiency between the spike-in sets.
In our mathematical framework, the variance of $\theta^*_i$ provides an upper bound for the variance of $F_i$.
Our estimate of $\variance(\theta^*_i)$ is an order of magnitude lower than the estimated variances of the log-size factors in each plate (Figure~\ref{fig:varestimates}) and of cellular RNA content (Supplemental Figure~\suppfigcell{}).
This indicates that the potential variance in the differences in spike-in behaviour, while greater than $\sigma^2_{vol}$, is still relatively small compared to other biases in the system, e.g., fluctuations in cellular RNA content, well-to-well variability in global capture efficiency.
To elaborate, consider that the maximum $\variance(F_i)$ corresponds to an error of $2^{\sqrt{0.015}} \approx 8$\% in the size factor estimates.
\revised{This error is small, especially} when we consider that spike-in normalization involves scaling the counts for each cell by at least $2^{\sqrt{0.15}} \approx 30$\% \revised{in our data sets.
These results suggest that variance in spike-in behaviour across cells is unlikely to have a strong effect on scaling normalization.}  

Here, $F_i$ is only computed between two spike-in sets.
In practice, the more relevant differences are those between synthetic spike-in and endogenous transcripts.
The variance of such differences is likely to be larger than $\variance(F_i)$, given the greater variability in sequence composition and length of endogenous transcripts.
Nonetheless, the SIRV and ERCC spike-ins do exhibit some variability in their biophysical properties (Supplemental Figure~\suppfigbiophys{}).
For example, the SIRV transcripts have more variable length and lower GC content compared to the ERCC transcripts.
This suggests that $F_i$ will include at least some of the differences in behaviour between synthetic and endogenous RNA, such that $\variance(F_i)$ can be used as a rough estimate of the magnitude of the associated variability.

\subsection*{Quantifying the effect of stochastic noise during sequencing}
We also performed simulations to gauge the contribution of $\sigma^2_{lib(s)}$ to $\variance(\theta^*_i)$ (see Section~\suppsecnoise{} of the Supplemental Materials).
Counts for spike-in transcripts were simulated such that any variability in the log-ratios was only caused by stochastic sampling noise, i.e., $\sigma^2_{lib(1)} + \sigma^2_{lib(2)}$.
Our results suggest that much of the estimated variance of $\theta^*_i$ in Figure~\ref{fig:varestimates} is driven by sampling noise (Supplemental Figure~\suppfignoise{}).
Specifically, we estimated the variance due to sampling noise to be 0.005-0.012 (using the original spike-in coverage for each plate), compared to estimates of 0.010-0.015 for $\variance(\theta^*_i)$ in Figure~\ref{fig:varestimates}.
Both $\sigma^2_{lib(1)} + \sigma^2_{lib(2)}$ and  $\variance(F_i)$ contribute to  $\variance(\theta^*_i)$, so these results suggest that the contribution of sampling noise is comparable to or greater than the impact of differences in spike-in behaviour.

We also observed that the variance due to sampling noise was robust to moderate decreases in the coverage of the spike-in transcripts in this simulation.
In ideal experiments, spike-in transcripts would take up 5-10\% of the library size for each cell (50,000-100,000 reads in our data).
Upon decreasing coverage \textit{in silico}, we observed an increase in $\sigma^2_{lib(1)} + \sigma^2_{lib(2)}$ due to the elevated effect of noise at low counts (Supplemental Figure~\suppfignoise{}).
However, even at 40-50\% coverage, the variance due to noise was still an order of magnitude lower than the variance due to cell-specific biases (Figure~\ref{fig:varestimates}b) or biological variability (Supplemental Figure~\suppfigcell{}).
These results suggest that spike-in normalization is still reliable when lower amounts of spike-in RNA are added.
This is especially relevant to data sets where the spike-in coverage is lower than recommended, due to difficulties in determining the appropriate concentration of spike-ins to add to each cell when the quantity of endogenous RNA is unknown.

\revised{Finally, we performed simulations to assess the effect of noise on the precision of the spike-in size factors themselves (see Section~\suppsecnoisier{} of the Supplemental Materials).
    This was performed using our 416B and TSC data sets as well as public data from existing studies \autocite{scialdone2015computational,zeisel2015brain,islam2014quantitative,hashimshony2016celseq2,buettner2015computational,grun2014validation,kolod2015single}.
In each data set, we observed that sampling noise resulted in $\approx 5$\% error in the estimates for the spike-in size factors (Supplemental Figure~\suppfignoisier{}).
By comparison, the size factors routinely varied by more than 30\% across cells (Supplemental Figure~\suppfigsizevar{}).
Thus, the loss of precision due to noise is small and can probably be ignored during spike-in normalization.
}

\subsection*{Assessing the downstream effect of variability with simulations}
We assessed whether the results of downstream analyses using spike-in normalization were sensitive to fluctuations in the total spike-in counts due to variability in spike-in addition, behaviour or sequencing noise.
First, we obtained data from \revised{scRNA-seq experiments that used spike-in RNA}.
This included \revised{a number of} public data sets \autocite{scialdone2015computational,islam2014quantitative,buettner2015computational,grun2014validation,kolod2015single,segerstople2016single} as well as our 416B and TSC data.
We then performed analyses such as detection of differentially expressed genes (DEGs) and highly variable genes (HVGs), as well as dimensionality reduction and clustering of cells.
This was done without any modification of the data to obtain a set of ``original results''.

Next, we designed simulations based on each of the real data sets (see Methods).
Briefly, the total spike-in count for each well was rescaled by a randomly sampled factor with variance equal to our experimental estimate of spike-in variance.
Counts for the individual spike-in transcipts were rescaled to reflect this new total, thus yielding a simulated data set.
\revised{Downstream} analyses were performed \revised{using the original counts for the endogenous genes and the simulated counts for the spike-in transcripts.
The new results were then compared to the original set of results from each analysis.}
Any differences indicate that the analysis is sensitive to spike-in variability in real experiments.
The advantage of this simulation design is that only the spike-in counts are modified.
\revised{No simulations or resampling were performed for the counts of the endogenous genes}, preserving the realistic nature of the data in each simulation \revised{and ensuring that only spike-in variability can cause differences in the analysis results.}

For DEG detection, we applied edgeR \autocite{robinson2010edgeR} and MAST \autocite{finak2015mast} to the original and simulated data after spike-in normalization.
edgeR represents methods designed for DE analyses of bulk RNA-seq data, while MAST represents bespoke single-cell methods.
In both cases, we observed only minor ($<$ 5\%) changes to the set of significant DEGs upon introducing spike-in variability in each data set (Figure~\ref{fig:setchange}a). 
Similar results were also observed in the top 200 DEGs with the smallest $p$-values, with fewer than 10\% of the genes in the set changing across iterations in all scenarios.
For HVG detection, we used methods based on the coefficient of variation \autocite{brennecke2013accounting} or the variance of log-expression values \autocite{lun2016stepbystep}.
Again, only minor changes were observed in most data sets (Figure~\ref{fig:setchange}b), for both the set of significant HVGs and for the top 200 HVGs with the smallest $p$-values.
These results suggest that the detection and ranking of DEGs and HVGs are largely robust to variability in spike-in volume or behaviour.
Indeed, genes that were not consistently detected across simulation iterations tended to have weak log-fold changes for DEGs, or small biological components for HVGs (Supplemental Figure~\suppfiglostsim{}).
This is expected as genes on the borderline of significance are more susceptible to random fluctuations in the size factors.

\begin{figure}[btp]
    \begin{center}
        \includegraphics[width=\textwidth,trim=0mm 5mm 0mm 5mm,clip]{../simulations/pics/setplot.pdf}
    \end{center}
    \caption{Effect of spike-in variability on DEG or HVG detection in simulated data.
        (a) The percentage change in the set of DEGs detected in each data set at a FDR of 5\% by edgeR or MAST.
        This was also calculated for the top set of 200 DEGs with the smallest $p$-values.
        Simulations were performed to detect DEGs \revised{in our 416B data set after inducing expression of a \textit{CBFB-MYH11} oncogene compared to a mCherry control (see Methods);
        between mouse embryonic stem cells (mESCs) in G1 and G2/M phases of the cell cycle \autocite{buettner2015computational};
        or between mESCs cultured in different conditions - serum, ground state (2i) or alternative ground state (a2i) \autocite{grun2014validation,kolod2015single}.}
        (b) The percentage change in the set of HVGs detected in each data set at a FDR of 5\%,
        using the Brennecke \textit{et al.} method based on the squared coefficient of variation (CV$^2$) or with a method based on the variance of log-expression.
        This was also calculated for the top set of 200 HVGs with the smallest $p$-values.
        Simulations were performed to detect HVGs in our 416B and TSC data sets, \revised{in liver cells \autocite{scialdone2015computational} and in mESCs \autocite{kolod2015single}.}
        All values represent the mean of 20 simulation iterations, and error bars represent standard errors.
    }
    \label{fig:setchange}
\end{figure}

For dimensionality reduction, we restricted ourselves to principal components analysis (PCA) on the normalized expression profiles of all cells. 
While $t$-distributed stochastic neighbour embedding \autocite{van2008visualizing} is commonly used, its robustness is difficult to evaluate due to its randomness.
\revised{We used a scRNA-seq study of the human pancreas \autocite{segerstople2016single} to generate} PCA plots of the first three principal components for both the original and simulated data.
At each simulation iteration, coordinates of all cells in the simulated plots were mapped onto the corresponding original plots to determine the sensitivity of the original locations to spike-in variability.
Figure~\ref{fig:dimclust}a indicates that changes in the location of each cell across simulation iterations were generally minor.
In particular, movement of cells across iterations did not compromise the separation of different cell types. 
Thus, spike-in variability does not appear to affect the visual interpretation of PCA plots.

\begin{figure}[btp]
    \begin{center}
        \includegraphics[width=\textwidth]{../simulations/clustering/clust_effect.pdf}
    \end{center}
    \caption{Effect of spike-in variability on dimensionality reduction and clustering in simulated data,
        based on real scRNA-seq data for cells extracted from a healthy human pancreas \autocite{segerstople2016single}.
        (a) PCA plots of the first three principal components, where each cell is coloured according to its annotated cell type from the original study. 
        The circle around each cell contains 95\% of remapped locations across the simulation iterations, and represents the deviation in location due to spike-in variability.
        (b) Clusters were identified from the original data by hierarchical clustering with Ward's criterion, followed by a tree cut with $k$ of 2, 5 or 10.
        This was repeated at each simulation iteration, and the maximum Jaccard index between each original cluster and any of the simulated clusters at the same $k$ was computed.
        Each value represents the mean of 20 simulation iterations, and the error bars represent standard errors.
        (c) The maximum Jaccard index for each original cluster generated with Ward's criterion compared to clusters generated from complete-linkage clustering of the original data.
    }
    \label{fig:dimclust}
\end{figure}

Finally, we performed hierarchical clustering and applied a tree cut to identify clusters of cells in the original data.
This was repeated at each simulation iteration to obtain a corresponding set of simulated clusters.
For each original cluster, we computed the Jaccard index with respect to each of the simulated clusters and recorded the maximum value across all simulated clusters.
A large maximum Jaccard index means that most of the cells in the original cluster are still grouped together in the simulation,
i.e., the original cluster is (mostly) successfully recovered in one of the simulated clusters.
We observed that the maximum Jaccard indices were moderate to large (Figure~\ref{fig:dimclust}b), with values above 0.6 for most of the original clusters.
To put this into context, we re-clustered the original data using a different algorithm.
This yielded smaller Jaccard indices for all clusters (Figure~\ref{fig:dimclust}c), indicating that spike-in variability has less effect on the results than the choice of clustering method.

\section*{Discussion}
In this article, we performed mixture experiments to quantify the variability of spike-in RNA addition across wells in a plate-based scRNA-seq protocol.
We also obtained a rough estimate of the well-to-well variability in the differences in behaviour between two different sets of spike-in transcripts.
Both values were at least an order of magnitude smaller than the variance of spike-in coverage across cells, suggesting that differences in spike-in volume or behaviour were not major sources of error in the context of spike-in normalization.
This was supported by simulations where the introduction of realistic levels of spike-in variance yielded only minor changes in the results of DEG and HVG analyses as well as PCA and clustering.
Our results indicate that spike-in normalization is reliable enough for routine use in scRNA-seq data analyses.
The common criticisms of using spike-in RNA for \revised{scaling} normalization are only weakly relevant, if at all, to single-cell transcriptome studies, and can generally be ignored.

Our conclusions differ from those of \citet{risso2014normalization}, where spike-in normalization is not considered reliable enough for analyses of bulk RNA-seq data.
We speculate that this difference may be due to the difficulty of adding an appropriate amount of spike-in RNA at the population level.
For example, should spike-in RNA be added at a constant ratio with respect to the concentration of endogenous RNA, or to the number of cells in the sample?
If the endogenous RNA concentration or the number of cells determines the amount of spike-in RNA to be added, these will need to be experimentally quantified for each sample.
In that case, how accurate is the quantification, and what effect do errors have on the downstream analysis?
These questions are not relevant to single-cell experiments where the obvious approach is to add the same amount of spike-in RNA to each individual cell.

It is worth mentioning another common criticism of the use of spike-in RNA -- namely, that the optimal concentration (to reach the suggested 5-10\% of total library size) depends on the amount of endogenous RNA in each cell.
This is not straightforward to gauge for biological systems that are not well studied.
If insufficient spike-ins are added, sampling noise will increase and the precision of spike-in normalization will deteriorate.
However, we do not consider this to be a matter of reliability.
It is \revised{not} surprising that suboptimal performance is obtained when inappropriate concentrations of reagents are used, and spike-in RNA is no exception.
Pilot experiments can be performed to identify the most suitable spike-in concentration to use for a given biological system, just like they would be used to determine the optimal dissociation conditions, lysis buffer, amplification cycles \revised{and so on}.
Once the optimal concentration is determined, it can be used for all experiments on that system and the choice of concentration ceases to be an issue.
In contrast, the variability of volume addition and of spike-in behaviour cannot be easily controlled even if all other parameters are optimized. 

We have used the Smart-seq2 protocol in our study to reflect its widespread use in the scRNA-seq literature.
However, our estimate of $\sigma^2_{vol}$ is agnostic to how reverse transcription, amplification and sequencing were performed, as these steps are represented by other mathematical terms.
Thus, we expect our conclusions to be broadly applicable to any scRNA-seq protocol where spike-in RNA is added in a similar manner (using repeater pipettes, see Methods).
Different results will be obtained using other methods for spike-in addition, e.g., with robotics systems or microfluidics, where volume handling may be even more precise.
Our experimental framework may also be useful for evaluating the precision of spike-in addition when developing new scRNA-seq protocols or setting up existing protocols in new laboratories, to ensure that spike-in RNA is added correctly to each cell.

The term $\variance(F_i)$ represents the variability in the difference in behaviour between the SIRV and ERCC spike-in sets across wells.
However, arguably a more relevant quantity is the variability in the difference $P_{is}$ between synthetic spike-in and endogenous RNA, as this affects the accuracy of normalization.
It may be possible to obtain a rough estimate of $\variance(P_{is})$ by using pooled cellular RNA from another organism as one of the spike-in sets \autocite{brennecke2013accounting}, so that $\variance(\theta^*_i)$ provides an upper bound on the variance in the differences in behaviour between synthetic and endogenous RNA.
We chose not to do so because of the difficulty in reproducibly using the same pool of cellular RNA across batches, and in calibrating the concentration of RNA to be added to each well.
Use of UMI counts may also provide a tighter bound on $\variance(F_i)$ or $\variance(P_{is})$ by reducing the contribution of amplification noise to $\sigma^2_{lib(s)}$ and $\variance(\theta^*_i)$.
Another consideration with endogenous RNA is the variability of lysis between cells, which we neglect to consider in our framework -- this is inherently difficult to assess with external spike-in RNA and may require other methods to quantify.

One interesting question is how to choose between spike-in normalization and approaches that assume a non-DE majority of genes.
This choice depends on whether total RNA content in each cell is of interest \autocite{lun2016stepbystep}.
Spike-in normalization will preserve changes in total RNA content between cells, whereas non-DE methods will treat such changes as bias (as a majority of genes are affected) and remove them.
This suggests that spike-in normalization is preferable in applications where changes in total RNA content can be easily associated with a biological process, e.g., T cell activation,   cell cycling.
In contrast, non-DE normalization may be more suitable for comparisons between distinct cell types, where the up- or downregulation of specific genes (conditional on the total RNA content of the cell) is more informative.
Obviously, this choice is subject to the experimental context.
If a non-DE majority cannot be assumed, spike-in normalization should be used; conversely, spike-in RNA cannot be easily added in droplet-based techniques, thus requiring non-DE methods.

We stress that our study only examines the reliability of spike-ins for ``relative'' normalization, i.e., to make counts comparable across cells.
We do not consider the reliability of spike-ins for absolute quantification, i.e., to determine the number of molecules of each transcript in each cell.
This is more difficult to evaluate as accuracy is affected by the magnitude of the differences in the behaviour of spike-in and endogenous transcripts.
In contrast, relative normalization is only affected by variability in the differences in behaviour across wells, as \revised{we have previously} discussed.

\section*{Methods}

\subsection*{Obtaining and culturing 416B cells and TSCs}

The murine multipotent myeloid progenitor cell line 416B \autocite{dexter1979isolation} was stably transduced with a TetOn construct of the \textit{CBFB-MYH11} (CM) oncogene (type A cDNA), using an in-frame F2A-mCherry protein as a reporter. 
As a control, cells were alternatively transduced with a version of the construct lacking the CM cDNA. 
Cells were maintained in RPMI medium, supplemented with 10\% fetal calf serum and antibiotics.
Expression of the CM oncogene or the mCherry control was induced by treatment with 1 \textmu{}g/ml of doxycycline, and induction was confirmed after 24 hours by measurement of mCherry levels by fluorescence activated cell sorting (BD Fortessa).

Murine TSCs were kindly provided by Dr.\ Jennifer Nichols (Wellcome Trust and MRC Cambridge Stem Cell Institute) and cultured by Liliana Antunes (Wellcome Trust Sanger Institute) on mouse embryonic fibroblast (MEF) feeders with TSC culturing medium (a combination of 70\% MEF conditioned media (R\&D systems) and 30\% RPMI 1640, supplemented with 20\% FBS, 2 mM L-glutamine, 1 mM sodium pyruvate, 100 \textmu{}M \textbeta{}-mercaptoethanol, 25 ng/mL human recombinant FGF4 (R\&D systems) and 1 \textmu{}g/mL heparin (Tocris Bioscience)). 
To prepare for single-cell sorting, cells were harvested with trypsin and MEF feeders were depleted by plating the cells onto a gelatinised plate followed by incubation for 1h at 37$^{\circ}$C on TSC culturing medium. 
The supernatant containing TSCs was used for sorting.

\subsection*{Spike-in mixture experiments with Smart-seq2}
% On a single plate, so no plate effects.
% Mention cutting up of plate, especially for the two conditions.
% No effect of different cell treatments.

Single-cell RNA sequencing was performed using an adaptation of the previously described Smart-seq2 protocol \autocite{picelli2014full}.
Single 416B cells or TSCs were sorted into individual wells of a 96-well microtiter plate.
Each well contained 2.3~\textmu{}l of lysis buffer with RNAse inhibitor (Ambion) in a 0.2\% (v/v) Triton X-100 solution. 
Reverse transcription (RT) was performed in a final volume of 13.2~\textmu{}l per well, containing 1~\textmu{}M of oligo-dT (Sigma-Aldrich), 1.04~mM of each dNTP (ThermoFisher), 100~U of SuperScript II retrotranscriptase (Invitrogen/ThermoFisher), 5~U of RNase inhibitor (Ambion), 5~mM of DTT, 1~M of Betaine (Sigma-Alrich), 6~mM of MgCl$_2$ (Ambion) and 1~\textmu{}M of TSO primer (Exiqon).
Preamplification was performed in a total volume of 27~\textmu{}l that contained 13.5~\textmu{}l of HiFi Hotstart ReadyMix (2$\times$; KAPA Biosystems) and 0.1~\textmu{}M of IS PCR primer (Sigma-Aldrich). 
After 23 cycles of amplification, samples were cleaned with 80\% (v/v) of Ampure beads (Beckman Coulter). 
Sequencing libraries were prepared using the Nextera XT DNA sample preparation kit (Illumina).
This was repeated to obtain several batches of sequencing data, with each batch consisting of one plate of cells of the same type.

To perform the mixture experiments, spike-in RNA was mixed into the RT reagent solution and added to each well.
This was done such that each well contained 0.1~\textmu{}l of a 1:3,000,000 dilution of the ERCC RNA Spike-In Mix (Invitrogen/ThermoFisher) and 0.12~\textmu{}l of a 1:3,000,000 dilution of the Spike-in RNA Variant (SIRV) Control Mix E0 (Lexogen). 
Two separate solutions of RT reagents were prepared for the different spike-in sets.
For one third of the wells, addition of the two spike-in sets was performed separately with the RT+ERCC solution first and the RT+SIRV solution second.
For another third of the wells, the order was reversed, i.e., with the RT+SIRV solution first and the RT+ERCC solution second.
For the remaining wells, the RT+SIRV and RT+ERCC solutions were premixed in a 1:1 ratio and the RT+SIRV+ERCC mixture was added twice to each well.
Each addition was performed independently for each well, using a repeater pipette dispensing 2~\textmu{}l at a time.

Sequencing of the 416B libraries was performed by the Genomics Core facility at the Cancer Research UK Cambridge Institute.
The first batch of libraries was sequenced on an Illumina HiSeq 2500 machine generating 125 bp single-end reads, while the second batch was sequenced on an Illumina HiSeq 4000 machine generating 50 bp single-end reads.
Sequencing of the TSC libraries was performed at the Wellcome Trust Sanger Institute after library preparation by the Single Cell Genomics Core facility.
Both batches were sequenced on an Illumina HiSeq 4000 machine generating 75 bp paired-end reads.

\subsection*{Data analysis for the mixture experiments}
Reads were mapped to the mm10 build of the mouse genome, including sequences of transcripts in the ERCC (\url{https://tools.thermofisher.com/content/sfs/manuals/ERCC92.zip}) and SIRV (\url{https://www.lexogen.com/wp-content/uploads/2015/11/SIRV_Sequences_151124.zip}) spike-in sets.
(The sequence of the \textit{CBFB-MYH11} oncogene was also included in the reference when aligning data from 416B cells.)
Mapping was performed using the subread aligner v1.5.1 \autocite{liao2013subread} in RNA-seq mode with unique alignment.
The 416B data were aligned in single-end mode while the TSC data were aligned in paired-end mode.
Reads with mapping qualities greater than or equal to 10 were assigned to exonic regions of genes using the featureCounts function in the Rsubread package v1.24.1 \autocite{liao2014featurecounts}.
Genes were defined using Ensembl v82 annotation for the GRCm38 mouse assembly and annotation for the ERCC and SIRV transcripts.
This yielded a count for each endogenous gene and spike-in transcript in each well.
Mapping and counting statistics for each batch of libraries are summarized in Supplemental Table~\supptabstats{}.

To evaluate spike-in quality, we verified that the total spike-in count (ERCC + SIRV) in each well comprised 5-10\% of the total library size (Supplemental Figure~\suppfigtotals{}).
This corresponds to the amount of spike-in RNA that we aimed to add to each well, and was consistent across wells within each plate.
The coverage of each ERCC transcript was directly proportional to its theoretical concentration in the spike-in mixture (Supplemental Figure~\suppfigspikeconc{}),
and the distribution of average read counts across spike-in transcripts or endogenous genes was consistent across plates (Supplemental Figure~\suppfigspikeave{}).
These diagnostics indicate that the spike-in transcripts were successfully captured, sequenced and processed into counts for most wells.
We removed any wells where the \revised{log-}total count for either spike-in set or for the endogenous genes was more than three median absolute deviations below the median value for each plate. 
It is likely that capture or sequencing failed for these wells, so they were not used for variance estimation.
In addition, we examined the effect of index switching \autocite{sinha2017index} in each data set generated on the HiSeq 4000, and found it to be negligible (see Section~\suppsecindex{} of the Supplemental Materials, Supplemental Figure~\suppfigindex{}).

Variance components were estimated from the libraries generated from a single plate.
In each well, the sum of counts across all transcripts in each spike-in set was computed, and the log$_2$-ratio between the ERCC and SIRV sums was calculated.
To estimate $\variance(\theta_i)$, a linear model with a one-way layout was fitted to the log-ratios for all wells where the two spike-in sets were added separately.
In each plate of the 416B data set, each combination of treatment (control or oncogene-induced) and spike-in addition order (ERCC or SIRV first) was treated as a group in the one-way layout.
In each plate of the TSC data, only the spike-in addition order was used to define the groups.
After fitting the model, the mean of the squared residual effects was used as an estimate of $\variance(\theta_i)$.
This was repeated for $\variance(\theta^*_i)$ using all wells where premixed spike-ins were added.
Here, addition order was irrelevant so the one-way layout contained only the two treatment groups in the 416B data set.
Similarly, only a single group was defined for the TSC data.
Linear modelling ensures that any changes in the mean log-ratio across groups do not inflate the variance estimate.
Note that we fit linear models to each plate separately, to check whether the estimates are consistent across replicate plates.

To detect differences in the variance estimates for premixed and separate addition, an F-test for the equality of variances was applied.
Under the null hypothesis of equal variances computed from independent data, the ratio of the variances $\sigma^2_1/\sigma^2_2$ should follow a F-distribution on $n_1$ and $n_2$ degrees of freedom, where $n_1$ and $n_2$ are the residual degrees of freedom used to estimate $\sigma^2_1$ and $\sigma^2_2$, respectively.
This can either be one-sided (i.e., $\sigma^2_1 \le \sigma^2_2$ under the null), in which case the upper tail probability at the observed ratio is taken as the $p$-value;
or it can be two-sided, in which case the $p$-value is defined as twice the smaller of the two tail probabilities.
Significant differences were defined by rejecting the null hypothesis at a type I error rate of 5\%.
We calculated $\sigma^2_{vol}$ from estimates of $\variance(\theta_i)$ and $\variance(\theta^*_i)$, using the expression described above.
However, if the difference between $\variance(\theta_i)$ and $\variance(\theta^*_i)$ was negative, $\sigma^2_{vol}$ was set to zero instead.
To assess the effect of the order of spike-in addition, a linear model was fitted to the subset of relevant wells on each plate to obtain an order-specific variance estimate.

\subsection*{Simulation design for resampling spike-in variability}
For each data set, we compute $T_{is}$ for each cell $i$ and spike-in set $s$.
To simplify the design of the simulations, we only consider the ERCC spike-in set here, i.e., $s=1$.
The variance of $T_{is}$ is
\[
    \variance(T_{is}) \approx \sigma^2_{lib(s)} + \sigma^2_{vol} + \variance(\log_2 R_{is}) + \variance(\log_2 L_i)
\]
where the approximation assumes that $L_i$ is independent of the other random variables that contribute to $T_{is}$.
(This is discussed in more detail in Section~\suppsecmath{} of the Supplemental Materials.)
Let $R_{is} = R_{i0}P_{is}$, where $R_{i0}$ is a random variable representing the well-specific average capture efficiency of endogenous transcripts and $P_{is}$ is the fold change in average efficiency of the transcripts in $s$ over their endogenous counterparts.
We assume that $R_{i0}$ and $P_{is}$ are independent for each well, and that $\variance(\log_2 P_{is})$ can be approximated with $\variance(F_i)$,
i.e., the well-to-well variability in relative capture efficiency between the two spike-in sets is similar to that between spike-ins and endogenous transcripts.
This means that
\[
    \variance(T_{is}) \approx \sigma^2_{lib(s)} + \sigma^2_{vol} + \variance(F_i) + \variance(\log_2 R_{i0}) + \variance(\log_2 L_i) \;,
\]
i.e., the variance of $T_{is}$ is a sum of the variances of its component terms.
The above approximation allows us to account for the measured $\sigma^2_{vol}$, \revised{$\sigma^2_{lib(s)}$} and $\variance(F_i)$ when simulating new values for $T_{is}$.

% If library quantification has been performed, then L_i will be more likely to be correlated with volume because everything competes against each other.
% At least it won't be correlated with global capture efficiency, as this should cancel out.
% If it hasn't been performed, then L_i should be constant so it shouldn't matter.
% The independence assumption means that we do not have to rescale the cellular counts to reflect variable undersampling.
% Doing so would be problematic, as direct scaling of the counts would distort the empirical mean-variance relationship of the count data.

Let us denote $x^2 = \sigma^2_{vol} + \variance(F_i) + \sigma^2_{lib(s)}$, representing the total variance in the log$_2$-total count of one spike-in set $s$ due to variable addition, capture efficiency and sequencing noise.
We use the estimated $\variance(\theta^*_i) \approx 0.015$ in Figure~\ref{fig:varestimates}a as our estimate $\hat{x}^2$ of the upper bound of $x^2$.
This is based on the fact that $\sigma^2_{vol}$ is near-zero in Figure~\ref{fig:varestimates}a, while $\variance(\theta^*_i) = \sigma^2_{lib(1)} + \sigma^2_{lib(2)} + \variance(F_i)$ and thus provides an upper bound on $\variance(F_i) + \sigma^2_{lib(s)}$ for any $s$.
We also denote $\hat\sigma^2_s$ as the estimate of $\variance(T_{is})$ across wells, and $\hat\mu_s$ as the estimate of $E(T_{is})$.
For each well $i$, we compute a simulated log$_2$-total $T^*_{is}$ as
\[
    T^*_{is} = (T_{is} - \hat\mu_s)\sqrt{1-\frac{ \hat{x}^2}{\hat\sigma^2_s}} + \hat\mu_s + X_i
\]
where $X_i \sim \mbox{Normal}(0, \hat{x}^2)$ and is independently sampled for each well.
This approach ensures that $\variance(T^*_{is}) = \hat\sigma^2_s$.
In contrast, if $X_i$ were directly added to $T_{is}$, the variance of $T^*_{is}$ would be inflated as $x^2$ is already present in $\variance(T_{is})$, i.e., the contribution of spike-in variance would be doubled.

Counts for the library generated from each well were rescaled to reflect the new, simulated log-total.
A quantile adjustment approach was used to preserve the empirical mean-variance relationship.
Briefly, a negative binomial generalized linear model (NB GLM) was fitted to the counts across all wells for each spike-in transcript, using the glmFit function in edgeR \autocite{mccarthy2012differential, robinson2010edgeR} with a design matrix containing all experimental factors in the current data set.
The value of $T_{is}/\log_2(e)$ was used as the offset for well $i$ during GLM fitting.
The NB dispersion was also estimated for each transcript using the estimateDisp function without empirical Bayes shrinkage.
For each transcript $t$, we assumed that the count $y_{ti}$ for well $i$ was sampled from a NB distribution with mean equal to the corresponding fitted value of the GLM and dispersion equal to the estimated transcript-specific value.
We scaled the NB mean by $2^{T^*_{is} - T_{is}}$ to obtain a modified NB distribution.
Using the q2qnbinom function \autocite{robinson2008small}, we calculated the lower tail probability of $y_{ti}$ in the original distribution and identified the corresponding quantile with the same tail probability in the modified distribution.
This new quantile was used as the simulated count for transcript $t$ in $i$.

% The estimation of the true parameters conditions on the observed spike-in totals, and doesn't make the assumption that spike-in totals are constant.
% For example, if it turned out that you added half the amount of spike-in to one well, it wouldn't distort the estimation of the (conditional) mean.
% Okay, maybe the trended dispersion would be a bit weird, because technical variability should depend on the amount of RNA, but some inaccuracy there is forgivable.

\subsection*{Evaluating the robustness of DEG detection}
\revised{We used a number of data sets to test the effect of spike-in variability on DEG detection.
This included our 416B data, where DEGs were detected between control and oncogene-induced cells;
and public data sets involving mESCs, where DEGs were detected between G1 and G2/M phases of the cell cycle \autocite{buettner2015computational}
or between different culture conditions \autocite{grun2014validation,kolod2015single}.
Access to each public data set is described in Section~\suppsecdata{} of the Supplemental Materials.}
In each study, DEGs were detected between conditions using edgeR and MAST.
Implementation details of each method are provided in Section~\suppsecsim{} of the Supplemental Materials. 
Briefly, normalization was performed by scaling the counts (explicitly or via offsets) such that the spike-in totals were the same between cells.
The set of DEGs in the original data was then identified at a FDR of 5\%.
This procedure was repeated for the simulated data, and the number of genes that were detected in the original results and not in the simulated results (or vice versa) was recorded as a proportion of the total number of original DEGs.
The proportion of the top 200 genes with the smallest $p$-values that were \revised{different} between the original and simulated results was also computed.
This was repeated for 20 simulation iterations and the average proportion across iterations was reported for each method.

\subsection*{Evaluating the robustness of HVG detection}
\revised{We used several data sets to test the effect of spike-in variability on HVG detection.
This included our 416B and TSC data sets, as well as public data sets involving mESCs \autocite{kolod2015single} or liver cells \autocite{scialdone2015computational}.}
In each data set, spike-in normalization was performed and HVGs were detected using two approaches based on spike-in counts.
The first approach is based on the method of \cite{brennecke2013accounting} where the squared coefficient of variation for each gene is tested for a significant increase above technical noise.
The second approach is based on the variance of the log-normalized expression values \autocite{lun2016stepbystep}, which provides some more robustness against outlier expression patterns.
Each method was applied on the original and simulated data, and a set of significant HVGs was detected at a FDR of 5\%.
The proportion of HVGs \revised{that were detected in the original results and not in the simulated results (or vice versa) was computed.}
\revised{Similarly, we computed} the proportion of the top 200 genes with the lowest $p$-values \revised{that differed between the original and simulated results.}
This was repeated for 20 simulation iterations and the average proportion across iterations was reported for each method.
\revised{See Section~\suppsecdata{} of the Supplemental Materials for details on public data access and Section~\suppsecsim{} for the implementation details of each HVG detection method.}

\subsection*{Evaluating dimensionality reduction and clustering}
\revised{We obtained count data from a study of pancreatic islet cells \autocite{segerstople2016single} (see Section~\suppsecdata{} of the Supplemental Materials).}
Spike-in normalization was performed and a set of HVGs was defined using the variance-of-log-expression method.
PCA plots of the first three components were constructed from the matrix of log-expression values for the HVGs.
This process -- including HVG detection -- was repeated with the simulated data after introducing spike-in variability.
To compare each simulated PCA plot to the original plot, the coordinates of each cell in the former were mapped onto the latter by rescaling and rotation.
Robustness was assessed based on the spread of remapped coordinates across all simulation iterations for each cell.
See Section~\suppsecsim{} in the Supplemental Materials for details.

To test the robustness of clustering, the matrix of Euclidean distances between cells was computed from the HVG log-expression values. 
Hierarchical clustering was performed using the Ward criterion and the resulting dendrogram was cut into 2, 5 or 10 clusters.
(This was done using the hclust and cutree commands, respectively, from the stats package.)
This process was repeated with the simulated data, and the Jaccard index between every pair of simulated and original clusters was computed.
For each original cluster, the maximum Jaccard index across all simulated clusters was recorded at each simulation iteration.
This value represents the extent to which the membership of the original cluster was preserved in the most similar simulated cluster.
We also compared the original clusters to those generated from complete-linkage clustering of the original HVG log-expression values.

\section*{Data access}
Data are available from the ArrayExpress database (\url{https://www.ebi.ac.uk/arrayexpress/}) using the accession E-MTAB-5522.
The R code used for the statistical analysis and simulations are available in the Supplemental Code, or at \url{https://github.com/MarioniLab/SpikeIns2016}.

\section*{Acknowledgements}
We thank Jennifer Nichols and Liliana Antunes for supplying the TSCs.
We also thank Victoria Moignard and Wajid Jawaid for helpful discussions about the experimental design.

\section*{Author contributions}
ATLL proposed the mixture experiments and performed the statistical analysis and simulations.
FJCN adapted the Smart-seq2 protocol for the spike-in mixtures and generated the 416B data.
LHV generated the TSC data.
BG and JCM provided direction and guidance on the project.
All authors wrote and approved the manuscript.

\section*{Disclosure declaration}
No conflicts of interest are declared.

\section*{Funding statement}
This work was supported by Cancer Research UK (core funding to JCM, award no.\ A17197), the University of Cambridge and Hutchison Whampoa Limited.
JCM was also supported by core funding from EMBL.
LHV was supported by an EMBL Interdisciplinary Postdoctoral fellowship, supported by H2020 Marie Sklodowska-Curie Actions.
Work in the G\"ottgens group was supported by Cancer Research UK, Bloodwise, the National Institute of Diabetes and Digestive and Kidney Diseases, the Leukemia and Lymphoma Society and core infrastructure grants from the Wellcome Trust and the Medical Research Council to the Cambridge Stem Cell Institute.

\section*{Supplemental materials}
The Supplemental Materials is a single PDF file that consists of Sections~1-\suppsecsim{} and contains Supplemental Figures~1-\suppfigindex{} and Supplemental Table~1.
The Supplemental Code is a compressed TAR file that contains the code and instructions required to reproduce the results in this paper.

\printbibliography

\end{document}
