\documentclass{article}
\usepackage{graphicx}
\usepackage[margin=2.5cm]{geometry}
\usepackage[labelfont=bf]{caption}
\usepackage{subcaption}
\usepackage{amsmath}
\usepackage{textcomp}

% 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}

% Pushing all floats to the bottom.
\renewcommand{\textfraction}{1.0}
\renewcommand{\floatpagefraction}{.9}
\renewcommand{\topfraction}{0.9}    % max fraction of floats at top
\renewcommand{\bottomfraction}{0.8} % max fraction of floats at bottom
\renewcommand{\textfraction}{0.07}  % allow minimal text w. figs

% Renaming float captions.
\makeatletter 
\renewcommand{\fnum@figure}{Supplemental \figurename~\thefigure}
\renewcommand{\fnum@table}{Supplemental \tablename~\thetable}
\makeatother

% URL linking.
\usepackage{url}
\urlstyle{same}

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

% Other convenience commands.
\newcommand{\figvarestimates}{2}
\newcommand\variance{\mbox{var}}

\begin{document}

\begin{titlepage}
\vspace*{3cm}
\begin{center}

{\LARGE
Assessing the reliability of spike-in normalization for analyses of single-cell RNA sequencing data
\par}

\vspace{0.75cm}

{\Large 
    \textsc{Supplemental Materials}
\par
}
\vspace{0.75cm}

\large
by


\vspace{0.75cm}
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}$

\vspace{1cm}
\begin{minipage}{0.9\textwidth}
\begin{flushleft} 
$^1$Cancer Research UK Cambridge Institute, University of Cambridge, Li Ka Shing Centre, Robinson Way, Cambridge CB2 0RE, United Kingdom \\[6pt]
$^2$Wellcome Trust and MRC Cambridge Stem Cell Institute, University of Cambridge, Wellcome Trust/MRC Building, Hills Road, Cambridge CB2 0XY, United Kingdom \\[6pt]
$^3$EMBL European Bioinformatics Institute, Wellcome Genome Campus, Hinxton, Cambridge CB10 1SD, United Kingdom \\[6pt]
$^4$Wellcome Trust Sanger Institute, Wellcome Genome Campus, Hinxton, Cambridge CB10 1SA, United Kingdom \\[6pt]
\end{flushleft}
\end{minipage}

\vspace{1.5cm}
{\large \today{}}

\vspace*{\fill}
\end{center}
\end{titlepage}

\section{Further interpretation of the mathematical terms}
$R_{is}$ was introduced as the average capture efficiency in well $i$ for all transcripts in set $s$.
It ranges from 0 to 1 and scales $r_{t_s}$ to determine the actual capture rate for each $t_s$.
The most obvious interpretation of $R_{is}$ (and $r_{t_s}$) is that of the efficiency of reverse transcription, but it also describes the efficiency of PCR amplification and tagmentation in Smart-seq2.
Thus, no additional variables are necessary for the latter steps.

The term $l_s L_i$ describes the rate at which reads are obtained from cDNA fragments during high-throughput sequencing.
The $l_s$ constant represents the average sequencing efficiency for transcripts in $s$, as well as factors such as mappability that affect the final counts.
The interpretation of $L_i$ depends on whether library quantification was performed to equalize the amount of cDNA from each well prior to sequencing.
If not, $L_i$ will be constant across wells, with its exact value depending on the sequencing depth.
However, if quantification was performed, $L_i$ will theoretically depend on the other variables that contribute to $T_{is}$.
Specifically,
\[
    L_i \approx D_i \left[ \sum_{s \in \{1, 2\}} \left( l_s V_{is} R_{is} \sum_{t_s} r_{t_s} c_{t_s} \right) + l_0 R_{i0} \sum_g r_g N_{ig} \right]^{-1}
\]
where $D_i$ is a random variable representing the total sequencing depth for well $i$ (in reads);
$N_{ig}$ is a random variable specifying the number of transcripts for each endogenous gene $g$ in $i$;
and $l_0$, $R_{i0}$ and $r_g$ are the equivalents to $l_s$, $R_{is}$ and $r_{t_s}$ for endogenous genes.
In practice, $L_i$ is effectively independent of $V_{is}$ and $R_{is}$ for any particular spike-in set $s$.
This is because the denominator of the above expression is dominated by the vastly larger number of cDNA fragments from the set of endogenous genes.
Any correlation between $L_i$ and the other terms would be negligible compared to the biological variance of expression, i.e., $\variance{(N_{ig})}$.

The error term $\varepsilon_{is}$ denotes the variability due to sequencing noise in the counts for spike-in set $s$.
We defined its variance as $\sigma^2_{lib(s)}$, which implicitly assumes that there is no relationship between the variance and the mean of $T_{is}$.
This is a strong assumption given that mean-variance relationships are often observed in RNA-seq data \citep{mccarthy2012differential,law2014voom}.
However, we note that the distribution of $T_{is}$ across wells with separate addition of spike-ins is similar to that across wells with premixed addition on the same plate (Supplementary Figure~\ref{fig:totals}).
If the mean of $T_{is}$ does not change between separate/premixed experiments, neither will the value of $\sigma^2_{lib(s)}$, regardless of the nature of the mean-variance relationship. 
This suggests that $\sigma^2_{lib(s)}$ will not change between $\variance(\theta_i)$ and $\variance(\theta^*_i)$, allowing calculation of $\sigma^2_{vol}$ from their difference.

\section{Computing the contribution of stochastic noise}

\subsection{Calculating the effect of noise on the log-ratios}
We performed simulations to obtain a rough estimate of the contribution of $\sigma^2_{lib(s)}$ to the variance of $\theta^*_i$.
For each plate in each data set (416B or TSC), we extracted the spike-in count data for the premixed wells.
For each spike-in set (SIRV or ERCC), the log-transformed sum of counts across all transcripts in that set was used as the GLM offset. 
The NB dispersion was estimated for each spike-in transcript using the estimateDisp function from edgeR.
We fitted a GLM to the counts for each transcript using a design matrix with a one-way layout.
This was set up using two groups of induced or control cells for each 416B plate, and with one group containing all cells for each TSC plate.
To simulate counts, we recalculated the fitted values of the GLM after setting the offsets of all cells to be equal to the mean offset.
The cell- and transcript-specific fitted values and the transcript-specific dispersion were used as the parameters for a NB distribution from which counts were resampled.
In this manner, simulated counts were obtained for each spike-in transcript in each cell.

Our aim is to compute the variance of the log-ratios of the total counts between spike-in sets from the simulated data.
By forcing the offsets of all wells to be equal, we removed differences in the log-ratios between wells due to cell-specific capture efficiency and sequencing depth.
This means that any variance in the simulated log-ratios across wells must be due to stochastic sampling noise.
Estimates of this variance represent the approximate contribution of $\sigma^2_{lib(1)} + \sigma^2_{lib(2)}$ to $\variance(\theta^*_i)$ (Supplemental Figure~\ref{fig:sampledepth}, at 100\% coverage).

We also downscaled the fitted values to determine the effect of reducing spike-in coverage on the sampling noise.
Specifically, we repeated the simulations after scaling the fitted values down to 1\% of their original values, to mimic the addition of less spike-in RNA to each well.
To preserve the mean-variance relationship, we fitted a loess curve of degree 1 to the original log-transformed NB dispersion against the log-mean of the fitted values across all cells.
We used this trend to determine the new NB dispersion for each transcript after downscaling.
Count sampling was then performed using the scaled fitted values and new dispersions.

We did not examine the effect of reducing spike-in coverage on $\sigma^2_{vol}$ or $\variance(F_i)$.
Both of these terms refer to experimental processes that should be independent of the concentration of spike-in RNA.
For example, the process of adding a volume of a solution is unrelated to the solution's contents.
Similarly, the process of capturing a transcript molecule is unrelated to the number of other transcript molecules, excepting the presence of limiting reagents.
Thus, the true values of $\sigma^2_{vol}$ or $\variance(F_i)$ should generally be unaffected by the depth of coverage.
In contrast, sampling noise during sequencing depends on the coverage and will contribute to the true value of $\sigma^2_{lib(s)}$.
This motivated our simulations to test the effect of reducing coverage on $\sigma^2_{lib(s)}$.

\subsection{Calculating the effect of noise on size factor precision}
The above simulations describe the contribution of stochastic noise to the log-ratios in our mixture experiments.
However, we can also perform simulations to quantify the effect of noise on the precision of the spike-in size factors themselves.
For a given data set, we first performed quality control by removing cells with small outlier values (i.e., more than 3 median absolute deviations below the median) for the log-transformed total endogenous counts or for the log-total spike-in counts.
We estimated the NB dispersion and fitted a GLM to the spike-in counts for each transcript, using the log-sum of spike-in counts as the GLM offset for each cell.
For each transcript in each cell, we used the corresponding fitted value and dispersion estimate to sample a new count from a negative binomial distribution, yielding a new set of simulated spike-in counts.

To estimate a size factor for each cell, we calculated the sum of simulated counts for all spike-in transcripts in each cell.
The count sums were scaled so that the average value across all cells in the data set was equal to unity.
The scaled count sum for each cell was then used as its spike-in size factor.
We repeated this procedure with a new set of simulated counts for multiple iterations, and we defined the estimation error as the standard deviation of the estimated size factors for the same cell across iterations. 
This was expressed as a percentage of the original value of the size factor (i.e., computed from the original counts) for each cell.
We performed these simulations on the data sets in Section~\ref{sec:datasets} to obtain Supplementary Figure~\ref{fig:sampling}.

We stress that these simulations only examine the variability in spike-in counts due to random sampling noise. 
It is not possible to make any precise conclusions regarding variability in the \textit{amount} of spike-in RNA added to each cell.
To do so requires the use of spike-in mixtures, which are not available in public data.

\section{Examining index switching on the HiSeq 4000}
Cells from three of our plates were multiplexed and sequenced on the Illumina HiSeq 4000 machine using the ExAmp chemistry.
A recent study showed that the index primers used for barcoding each sample could switch between samples \citep{sinha2017index}, thereby contaminating one sample with DNA from other samples.
We examined the effect of index switching on our 416B data, exploiting the fact that the first plate was generated on the HiSeq 2500 while the second was generated on the HiSeq 4000.
The variance estimates were very similar between the two plates (Figure~\figvarestimates{}), suggesting that index switching has minimal effect on our conclusions.
We also made use of the presence of the \textit{CBFB-MYH11} oncogene, the expression of which was upregulated in half the cells on each plate.
We observed similar log-fold changes between induced and control cells for both plates (Supplementary Figure~\ref{fig:switching}a, b) whereas we would have expected that the log-fold change for the second plate would be closer to zero if samples were homogenized due to index switching.

For the TSC data set, each plate contained a negative control well where no cell was added.
Any counts for endogenous mouse genes in those control wells could be attributed to index switching (or to other sources of contamination, e.g., cell-free RNA).
In each plate, the total count across all mouse genes for each negative control well was around 2-3 orders of magnitude lower than the total mouse counts in the majority of wells (Supplementary Figure~\ref{fig:switching}c, d).
This suggests that, at worst, index switching would result in 1\% of the library size for each well being attributable to contamination from other wells.
To put this result into context, the coefficients of variation for the total ERCC count in the first and second plate are 0.29 and 0.26, respectively.
This means the coverage routinely varies by 26-29\% across wells due to cell-specific biases.
By comparison, the 1\% (maximum) rate of contamination due to index switching is likely to have negligible effect.

% x <- loadin1$counts[loadin1$genes$spike1,]
% x <- colSums(loadin1$counts[loadin1$genes$spike1,])
% sqrt(var(x)/mean(x)^2)
% x <- loadin2$counts[loadin2$genes$spike1,]
% x <- colSums(loadin2$counts[loadin2$genes$spike1,])
% sqrt(var(x)/mean(x)^2)

\section{Overview of data sets used in simulations}
\label{sec:datasets}
We obtained the following public data sets, containing counts for ERCC spike-in transcripts:
\begin{itemize}
    \item Liver cells from \citet{scialdone2015computational}, obtained from the ArrayExpress repository with the accession E-MTAB-3707.
        No experimental factors were considered.
    \item Brain cells from \citet{zeisel2015brain}, obtained from \url{http://linnarssonlab.org/cortex}.
        Only the subset of cells from the cortex was used.
        No experimental factors were considered.
    \item Mouse embryonic stem cells from \citet{buettner2015computational}, obtained from ArrayExpress with the accession E-MTAB-2805.
        Cell cycle phase (G1, G2M or S) was treated as the experimental factor.
    \item Mouse embryonic stem cells from \citet{kolod2015single}, obtained from \url{http://www.ebi.ac.uk/teichmann-srv/espresso}.
        The culture condition (2i, a2i or lif) was treated as the experimental factor.
        Only the subset of cells from batch 3 was used as these contained spike-ins for all conditions.
    \item Mouse embryonic stem cells from \citet{islam2014quantitative}, obtained from the Gene Expression Omnibus (GEO) with the accession GSE46980.
        No experimental factors were considered.
    \item Mouse embryonic stem cells from \citet{grun2014validation}, obtained from GEO with the accession GSE54695.
        The culture condition (2i or serum) was treated as a experimental factor.
    \item Fibroblasts from \citet{hashimshony2016celseq2}, obtained from GEO with the accession GSE78779.
        No experimental factors were considered.
        Only the subset of cells from the C1 experiment was used.
    \item Pancreatic islet cells from \citet{segerstople2016single}, obtained from ArrayExpress with the accession E-MTAB-5061.
        Cells annotated as low quality were removed.
        (For downstream analyses, we considered all remaining cells as high-quality, so no further cell-level quality control was applied.)
        We only used the cells extracted from a single individual (HP1502401, healthy male) for simplicity.
\end{itemize}
In each of our TSC and 416B data sets, data from both plates were analyzed together.
We treated the plate of origin for each cell as an experimental factor, along with the oncogene induction status for 416B cells.
Only the ERCC counts were used in simulations, and all SIRV transcripts were discarded for simplicity.

\section{Implementation details for the downstream analyses}

\subsection{Quality control on cells and genes}
Cell-level quality control was performed on each data set by removing cells with outlier values for various metrics \citep{lun2016stepbystep}.
These metrics included the log-transformed total read count across all genes and the log-transformed total number of expressed genes, where small outliers were removed; 
and the proportion of reads mapped to spike-in transcripts, where large outliers were removed.
Outliers were defined as values that were more than three median absolute deviations away from the median in the specified direction.
Genes were also removed if the average count across all cells was below 0.1.
This reduces computational work by filtering out low-abundance genes that do not contain enough information for reliable inference.

\subsection{Methods for detecting differentially expressed genes}
For DEG detection with edgeR v3.18.1, a NB GLM was fitted to the counts for each gene \citep{mccarthy2012differential} using a suitable design matrix.
The design matrix for each data set was constructed using an additive parameterization for all experimental factors, including the factor containing the conditions to be compared for differential expression.
The log-transformed total count for the spike-in transcripts was used as the offset for each library.
This is equivalent to spike-in normalization as each count is effectively downscaled by the spike-in total.
An abundance-dependent trend was fitted to the NB dispersions of all genes using the estimateDisp function.
For each gene, empirical Bayes shrinkage was performed towards the trend to obtain a shrunken NB dispersion, and the likelihood ratio test was applied to test for significant differences in expression between conditions.
The Benjamini-Hochberg correction was applied to control the FDR.

For MAST v1.2.1, an effective library size was defined by multiplying the total spike-in count by a constant value $C$ for each library.
We set $C$ to the ratio of the average total count for the endogenous genes to the average total count for the spike-in transcripts, where each average was computed across all libraries.
This procedure recapitulates the effect of spike-in normalization by ensuring that the fold-differences in the effective library sizes are equal to the fold-differences in the spike-in totals between cells.
Counts were converted to count-per-million (CPM) values using the effective library sizes, which were further log-transformed after adding a pseudo-count of 1.
For each gene, a hurdle model was fitted to the log-CPMs across all cells using the zlm function with default parameters.
This was parameterized as an additive model for all experimental factors, with the proportion of genes with non-zero counts in each library included as a covariate \citep{finak2015mast}.
Putative DE genes between the relevant conditions were identified using the lrTest function.

\subsection{Methods for detecting highly variable genes}
The first HVG detection method was based on the approach described by \citet{brennecke2013accounting}, with some modifications to the size factors to perform spike-in normalization.
Specifically, the spike-in totals were scaled so that their mean across all libraries was equal to unity.
The size factor for each library was defined as its scaled spike-in total.
We did \textit{not} compute separate size factors for the endogenous genes and spike-in transcripts, as this would require the use of non-DE normalization methods.
The rest of the method was implemented as originally described, using the technicalCV2 function in the scran package (v1.4.4) with min.bio.disp set to zero.
If any experimental factors were present, they were regressed out by log-transforming the counts; applying the removeBatchEffect function from the limma package \citep{ritchie2015limma} to the log-counts; and converting the corrected log-values back to the count scale, prior to using technicalCV2.

The second approach to detect HVGs was based on computing the variance of log-expression values \citep{lun2016stepbystep}.
For each count in each library, a normalized log-expression value was defined as the log-ratio of the count with the spike-in size factor for that library.
(A pseudo-count of 1 was added to avoid undefined values upon log-transformation.)
If any experimental factors were present, they were used to construct a design matrix with an additive parameterization -- otherwise, an intercept-only design matrix was used.
The design matrix was used to fit a linear model to the log-expression values to obtain the residual variance for each spike-in transcript.
A loess curve was fitted to the log-variance against the mean for all spike-in transcripts using the trendVar function in scran, representing the mean-variance relationship due to technical noise.
Variance estimation was repeated on the log-expression values for each endogenous gene, and the biological component of the variance and a $p$-value was computed for each gene using the decomposeVar function.

\subsection{Methods for dimensionality reduction and clustering}
Size factors for all libraries were defined from the spike-in totals as previously described. 
HVG detection was performed using the variance-of-log-expression method, where HVGs were defined as genes detected at a FDR of 5\% and with biological components above 0.5.
PCA was performed on the normalized log-expression values of the HVGs, using the prcomp function from the stats package with scaling.
The first two PCs were used as the coordinates for each cell in one PCA plot, and the first and third PCs were used as coordinates in another plot.
Each point was coloured according to its annotated cell type \citep{segerstople2016single}.

The procedure above was repeated at each simulation iteration with new spike-in counts.
Coordinates of all cells in each simulated PCA plot were mapped onto the original plot, after scaling and rotating the coordaintes to eliminate differences between plots that were not relevant to interpretation.
Specifically:
\begin{itemize}
    \item For each cell, the simulated coordinates were right-multiplied by a 2-by-2 transformation matrix
        \[
            \left[\begin{array}{cc}
                    b_x \cos(\psi) & - b_y \sin(\psi) \\
                    b_x \sin(\psi) & b_y \cos(\psi) 
                \end{array}
            \right]
        \]
        where $\psi$ is the angle around the origin and $b_x$ and $b_y$ scale the $x$- and $y$-coordinates, respectively.
    \item The squared Euclidean distance from the (scaled and rotated) simulated coordinates to the original coordinates was computed for each cell, and summed across all cells.
    \item The scaling and rotation parameters of the matrix were identified that minimized the sum of squared distances across all cells, using the optim function from the stats package.
\end{itemize}
The more obvious approach to remapping is to directly project the simulated log-expression data onto the space of the original plot.
However, we do not do this as it does not capture the variability in the identification of the PCs across iterations.
Upon completion of the simulation, each cell will have one original location and one remapped location per iteration. 
For each cell, the smallest circle centered at its original location was drawn that contained 95\% of the remapped locations.
This avoids inflated circles due to outliers.

\printbibliography

\newpage

\begin{figure}[btp]
    \begin{center}
        \includegraphics[width=\textwidth]{pics/explanation.pdf}
    \end{center}
    \caption{Schematic of the effect of variability in spike-in addition ($\sigma^2_{vol}$) or library preparation ($\sigma^2_{lib(s)}$) on the performance of spike-in normalization.
        The spike-in size factor for each cell is calculated from the sum of counts for spike-in transcripts.
        The size factors faithfully represent the true cell-specific biases when equal amounts of spike-in RNA are added to each well and technical noise in library preparation and sequencing is low.
        This is not true when spike-in addition is variable or technical noise is high, manifesting as errors relative to the true biases (i.e., large residual variance).
    Note that the true biases are unknown, so $\sigma^2_{vol}$ and $\sigma^2_{lib(s)}$ are instead estimated from spike-in mixtures as described in our experimental design.
    }
\end{figure}

\begin{figure}[btp]
    \begin{center}
        \begin{subfigure}{0.49\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 5mm 0mm 5mm,clip]{../real/pics/qq_separate.pdf}
            \caption{}
        \end{subfigure}
        \begin{subfigure}{0.49\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 5mm 0mm 5mm,clip]{../real/pics/qq_premixed.pdf}
            \caption{}
        \end{subfigure} \\[0.1in]
        \begin{subfigure}{0.49\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 5mm 0mm 5mm,clip]{../real/pics/dist_separate.pdf}
            \caption{}
        \end{subfigure}
        \begin{subfigure}{0.49\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 5mm 0mm 5mm,clip]{../real/pics/dist_premixed.pdf}
            \caption{}
        \end{subfigure}
    \end{center}
    \caption{
        Distribution of the log-ratios after separate or premixed addition of spike-in sets.
        For each plate, a linear model was fitted to the log-ratios of the corresponding wells to account for factors in the experimental design.
        (a, b) Quantile-quantile plots of the residuals of the fitted model for each plate.
        Residuals were standardized and plotted against the theoretical quantiles of a standard normal distribution.
        The dotted line represents equality between the sample and theoretical quantiles.
        (c, d) Density estimates of the residuals from each plate, computed using a Gaussian kernel with twice the default bandwidth.
    }
\end{figure}

\begin{figure}[btp]
    \begin{center}
        \begin{subfigure}{0.49\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 5mm 0mm 5mm,clip]{../real/pics/total_ercc.pdf}
            \caption{}
        \end{subfigure}
        \begin{subfigure}{0.49\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 5mm 0mm 5mm,clip]{../real/pics/total_sirv.pdf}
            \caption{}
        \end{subfigure}
        \begin{subfigure}{0.49\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 5mm 0mm 5mm,clip]{../real/pics/prop_ercc.pdf}
            \caption{}
        \end{subfigure}
        \begin{subfigure}{0.49\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 5mm 0mm 5mm,clip]{../real/pics/prop_sirv.pdf}
            \caption{}
        \end{subfigure}
    \end{center}
    \caption{Distribution of the coverage of the (a, c) ERCC or (b, d) SIRV spike-in set across wells, 
        shown as the total number of reads assigned to the spike-in transcripts (a, b) or as a proportion of the total number of reads in the corresponding library (c, d).
        For each plate, separate boxplots are shown for wells in which spike-ins were added separately or premixed before addition.
        Dots represent wells that are more than 1.5 interquartile ranges from the first or third quartile of the corresponding distribution.
    }
    \label{fig:totals}
\end{figure}

\begin{figure}[btp]
    \begin{center}
        \includegraphics[width=0.7\textwidth,trim=0mm 10mm 0mm 5mm,clip]{../real/pics/variance_order.pdf}
    \end{center}
    \caption{Estimated variance of the log-ratio of total counts between spike-in sets, computed across wells in which the ERCC spike-in set was added before the SIRV set or vice versa.
        This exploits the presence of a number of wells in each plate (indicated above each bar) for which the order of spike-in addition was reversed.
        Error bars represent standard errors of the variance estimates for normally distributed log-ratios.
        Differences between the ERCC- and SIRV-first estimates of each batch were assessed using a two-sided F-test, yielding $p$-values of 0.96, 0.66, 0.02 and 0.33 for the respective batches from left to right.
    }
\end{figure}

\begin{figure}[btp]
    \begin{center}
        \includegraphics[width=0.7\textwidth,trim=0mm 10mm 0mm 5mm,clip]{../real/pics/variance_cell.pdf}
    \end{center}
    \caption{Estimated variance of the log-ratio of total counts across all endogenous mouse genes to that across all ERCC transcripts, computed across all wells on each plate (numbers shown above each bar).
        Error bars represent standard errors of the variance estimates for normally distributed log-ratios.
    }
\end{figure}

\begin{figure}[btp]
    \begin{center}
        \includegraphics[width=0.49\textwidth,trim=0mm 10mm 0mm 10mm,clip,page=1]{../sequences/biophysical/comparison.pdf}
        \includegraphics[width=0.49\textwidth,trim=0mm 10mm 0mm 10mm,clip,page=2]{../sequences/biophysical/comparison.pdf}
    \end{center}
    \caption{Biophysical properties of transcripts in each of the two spike-in sets and for 2000 randomly selected transcripts from the mouse mm10 genome.
    Boxplots are shown for the distribution of lengths and GC contents of transcripts (not including the poly-A tail) in each set.
}
\end{figure}

\begin{figure}[btp]
    \begin{center}
        \begin{subfigure}{0.49\textwidth}
            \includegraphics[width=\textwidth,page=1]{../real/depth/depth_effect.pdf}
            \caption{}
        \end{subfigure}
        \begin{subfigure}{0.49\textwidth}
            \includegraphics[width=\textwidth,page=2]{../real/depth/depth_effect.pdf}
            \caption{}
        \end{subfigure}
    \end{center}
    \caption{Effect of spike-in coverage on the variance of the log-ratios of the ERCC and SIRV total counts, using simulations based on the (a) 416B or (b) TSC data sets.
        Simulations were performed where the only source of variation between cells was the stochastic sampling of counts from a NB distribution with parameters estimated from real data.
        This was repeated after scaling the mean count for each spike-in transcript to 1-100\% of the original coverage. 
        The last point of each line represents the variance at 100\% coverage on the corresponding plate.
        The total count of both spike-in sets was computed as an average across all cells on the plate, and is used here to quantify overall coverage of the spike-ins after scaling. 
        (See Supplemental Figure~\ref{fig:totals} for the original total count of each set.)
        Each point represents the mean variance from 20 simulation iterations, with error bars representing standard errors of the variance estimates.
    }
    \label{fig:sampledepth}
\end{figure}

\begin{figure}
    \begin{center}
        \includegraphics[width=\textwidth]{../simulations/sampling/collated_sampling.pdf}
    \end{center}
    \caption{Relative estimation error for the spike-in size factor of each cell, shown as a percentage of the value of the size factor.
        Boxplots represent the distribution of errors across cells in each data set.
        Errors were calculated after resampling spike-in counts from a negative binomial distribution with parameters estimated from the original counts.
        Calculations were performed using data for mouse embryonic stem cells (mESCs) from \citet{islam2014quantitative},
        \citet{buettner2015computational}, \citet{grun2014validation} and \citet{kolod2015single};  
        fibroblasts from \citet{hashimshony2016celseq2};
        trophoblast stem cells (TSCs) and 416B cells from our experiments;
        liver cells from \citet{scialdone2015computational};
        and brain cells from \citet{zeisel2015brain}. 
    }
    \label{fig:sampling}
\end{figure}

\begin{figure}
    \begin{center}
        \includegraphics[width=\textwidth]{../simulations/sampling/collated_sizefac.pdf}
    \end{center}
    \caption{Absolute deviation of the spike-in size factor of each cell, shown as a percentage of the mean size factor across cells.
        Deviations were calculated by fitting a linear model to the $\log_2$-size factors, using a design matrix containing all experimental factors in the data set.
    The deviation for each cell was defined as $|2^r-1| \times 100$ where $r$ is the corresponding residual for that cell in the fitted model.
    Each boxplot represents the distribution of deviations across cells in a particular data set, see Supplemental Figure~\ref{fig:sampling}.
    }
\end{figure}


\begin{figure}[btp]
    \begin{center}
        \begin{subfigure}{0.24\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 10mm 0mm 0mm,clip,page=1]{../simulations/pics/comparison_edgeR_DEG.pdf}
            \caption{}
        \end{subfigure}
        \begin{subfigure}{0.24\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 10mm 0mm 0mm,clip,page=2]{../simulations/pics/comparison_edgeR_DEG.pdf}
            \caption{}
        \end{subfigure}
        \begin{subfigure}{0.24\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 10mm 0mm 0mm,clip,page=1]{../simulations/pics/comparison_log_HVG.pdf}
            \caption{}
        \end{subfigure}
        \begin{subfigure}{0.24\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 10mm 0mm 0mm,clip,page=2]{../simulations/pics/comparison_log_HVG.pdf}
            \caption{}
        \end{subfigure}
    \end{center}
    \caption{Properties of genes that were detected as DEGs or HVGs in both the original and simulated data with spike-in variability (common), compared to genes that were no longer detected or additionally detected after simulation (lost/gained).
    (a, b) Absolute magnitudes of the log-fold changes of DEGs from edgeR in the 416B or mESC cell cycle data sets.
    (c, d) Biological components of variation of HVGs from the variance of log-expression method in the 416B and TSC data sets.
    The number of genes in each category is also shown for each data set.
    Similar results were obtained for the CV$^2$ method and MAST.
    }
\end{figure}

\begin{table}[btp]
    \caption{Alignment and counting statistics for each batch of scRNA-seq data, including the total number of fragments (reads for single-end data, read pairs for paired-end data), percentage of reads mapped to the reference genome and percentage of fragments assigned to genic regions.
    For each statistic, the median value across all wells in the batch is shown with first and third quartiles in brackets.}
    \begin{center}
        \begin{tabular}{l r r r}
            \hline
            \textbf{Data set} & \textbf{Total ($\times 10^6$)} & \textbf{Mapped (\%)} & \textbf{Counted (\%)} \\
            \hline
            416B (I)  & 2.80 (2.39-3.10) & 59.2 (56.6-61.6) & 46.3 (45.1-49.3) \\
            416B (II) & 2.82 (2.40-3.26) & 50.3 (47.3-53.1) & 39.0 (36.5-42.3) \\
            Tropho (I) & 2.02 (17.9-2.20) & 88.8 (88.1-89.3) & 74.9 (73.0-76.5) \\
            Tropho (II) & 2.33 (2.08-2.57) & 89.1 (87.6-89.7) & 62.8 (61.5-65.5) \\
            \hline
        \end{tabular}
    \end{center}
\end{table}

% blah <- read.table("my_qual.tsv", header=TRUE, comment="")
% totals <- rowSums(blah)
% mapped <- 1 - blah$Unassigned_Unmapped/totals
% counted <- blah$Assigned/totals
% summary(totals, digits=100)
% summary(mapped)
% summary(counted)

\begin{figure}[btp]
    \begin{center}
       \includegraphics[scale=0.5,page=1,trim=0mm 25mm 10mm 5mm,clip]{../real/pics/against_theoretical.pdf}
       \includegraphics[scale=0.5,page=2,trim=10mm 25mm 10mm 5mm,clip]{../real/pics/against_theoretical.pdf}
       \includegraphics[scale=0.5,page=3,trim=0mm 5mm 10mm 5mm,clip]{../real/pics/against_theoretical.pdf}
       \includegraphics[scale=0.5,page=4,trim=10mm 5mm 10mm 5mm,clip]{../real/pics/against_theoretical.pdf}
       \caption{Observed average read count for each transcript in the ERCC spike-in set, plotted against the theoretical concentrations (obtained as Mix A in ``ERCC Controls Analysis: ERCC RNA Spike-In Control Mixes (English)'' at https://www.thermofisher.com/order/catalog/product/4456740).
        Each point represents a spike-in transcript, with counts averaged across all wells in a plate.
        The red line represents a line of best fit with gradient 1.
        Equivalent plots for the SIRVs are not shown as the relationship between the observed count and the theoretical molar concentration is not straightforward in the presence of isoforms.
       }
    \end{center}
\end{figure}

\begin{figure}[btp]
    \begin{center}
        \includegraphics[width=0.5\textwidth]{../real/pics/feature_abundances.pdf}
    \end{center}
    \caption{Distribution of average counts across the ERCC and SIRV spike-in transcripts and endogenous mouse genes.
    For each feature, counts were averaged across all wells on each plate.}
\end{figure}

\begin{figure}[btp]
    \begin{center}
        \begin{subfigure}{0.49\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 10mm 0mm 0mm,clip,page=1]{../real/index_swapping/Calero_check.pdf}
            \caption{}
        \end{subfigure}
        \begin{subfigure}{0.49\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 10mm 0mm 0mm,clip,page=2]{../real/index_swapping/Calero_check.pdf}
            \caption{}
        \end{subfigure}
        \begin{subfigure}{0.49\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 5mm 0mm 0mm,clip,page=1]{../real/index_swapping/Liora_check.pdf}
            \caption{}
        \end{subfigure}
        \begin{subfigure}{0.49\textwidth}
            \includegraphics[width=\textwidth,trim=0mm 5mm 0mm 0mm,clip,page=2]{../real/index_swapping/Liora_check.pdf}
            \caption{}
        \end{subfigure}
    \end{center}
    \caption{Effect of index switching in each data set for the mixture experiments.
        (a, b) Log-fold change of \textit{CBFB-MYH11} oncogene expression in control and induced 416B cells on each plate, generated using the HiSeq 2500 (plate I) or 4000 (II).
        (c, d) Distribution of log-total counts for mouse genes for all wells on each TSC plate.
        The negative control well on each plate is marked by the red triangle.
    }
    \label{fig:switching}
\end{figure}

\end{document}


