• Role of read depth in experiment design
  • "Do you have enough reads?": Introducing subSeq software
  • "Do people usually have enough reads?": subSeq results across many RNA-Seq experiments
  • "Exactly how many reads do we need?": Predicting the curve

• Multiplexed barcode sequencing allows multiple libraries to be run on the same flowcell
• Investigators can choose how to allocate their reads
  • Additional replicates
  • Additional conditions/time points
  • Multiple experiments
• Designing the experiment requires knowing effect of read depth

• When designing a sequencing experiment, choose a read depth and number of replicates to satisfy any of the following requirements:
  • Number of significant genes
  • Statistical power for a given gene or set of genes
  • Statistical power for a given effect size
• Given an existing experiment, determine if we have a high enough depth to trust the biological conclusions
• Test how analysis methods are affected by read depth

• Given an \(m \times n\) matrix \(X\) of counts
• Construct a matrix \(Y\), where every read has a probability \(p\) of being included:

\[ Y_{i,j} \sim \mbox{Binom}(X_{i,j}, p) \]

• Mathematically identical to (e.g.) DownsampleSam software

• Unnormalized matrix of read counts (e.g. from htseq-count)
  • One row per gene
  • One column per sample
• Choice of subsampling proportions in (0, 1]
  • More proportions = more running time, more detail
• Choice of one or more analysis methods
  • Anything you can do on the full count matrix in R

library(subSeq)
data(hammer)
library(Biobase)
countmatrix <- exprs(hammer)[, 1:4]
countmatrix <- countmatrix[rowSums(countmatrix) >= 5, ]
head(countmatrix)

##                    SRX020102 SRX020103 SRX020104 SRX020105
## ENSRNOG00000000001         2         4        18        24
## ENSRNOG00000000007         4         1         3         1
## ENSRNOG00000000008         0         1         4         2
## ENSRNOG00000000010        19        10        19        13
## ENSRNOG00000000012         7         5         1         0
## ENSRNOG00000000017         4         1        12        18

We compare a control group with rats treated with L5 SNL.

treatment <- phenoData(hammer)$protocol[1:4]

treatment

## [1] control control L5 SNL  L5 SNL 
## Levels: control L5 SNL

We recommend choosing proportions on a logarithmic scale.

proportions <- 10^seq(-2, 0, .25)
proportions

## [1] 0.01000 0.01778 0.03162 0.05623 0.10000 0.17783 0.31623 0.56234 1.00000

Choose one or more RNA-Seq differential expression methods for each step.

subsamples <- subsample(countmatrix, proportions, method=c("edgeR", "DESeq2"),
                       treatment=treatment)

subsamples.summ <- summary(subsamples)
head(subsamples.summ)

##      depth proportion method replication significant pearson spearman
## 1:  198567    0.01000  edgeR           1         201  0.3296   0.3892
## 2:  352812    0.01778  edgeR           1         425  0.3920   0.4603
## 3:  627759    0.03162  edgeR           1         827  0.4606   0.5373
## 4: 1117230    0.05623  edgeR           1        1534  0.5503   0.6222
## 5: 1985611    0.10000  edgeR           1        2468  0.6417   0.7127
## 6: 3530798    0.17783  edgeR           1        3502  0.7456   0.8001
##       MSE   estFDP     rFDP percent
## 1: 2.3728 0.001191 0.000000 0.02403
## 2: 2.0724 0.004430 0.007059 0.05046
## 3: 1.7405 0.007126 0.008464 0.09805
## 4: 1.3607 0.005091 0.003911 0.18271
## 5: 1.0577 0.009132 0.009319 0.29236
## 6: 0.7355 0.013354 0.013992 0.41289

plot(subsamples.summ)

• Perform multiple replications of subsampling at each depth
• Use your own custom analysis methods (anything that can be done on a count matrix in R)
  • Alternative splicing/differential exon usage
• Recreate the exact count matrix used at any depth

• Repository of processed gene expression experiments
• We used 25 contrasts across 12 experiments
  • edgeR and DESeq2
  • Performed subsampling at 41 proportions, evenly spaced on a log scale from .01 to 1
  • 3 replications at each proportion

• Distribution of effect sizes
• Distribution of gene expression level
• Distribution of gene length
• Experiment design (# of replicates)
• Amount of variation added at each stage

\[ G(p) = k\times \Phi(log(p - b), \mu, \sigma)\]

where

• \(G(p)\) is the number of significant genes at proportion \(p\)
• \(\Phi\) is the CDF of the normal distribution.
• \(k,\mu,\sigma,b\) are the four parameters fit for each curve, found through nonlinear least squares

• subSeq can predict the effect of increasing the read depth
  • Can predict effect of doubling quite accurately
  • Tripling with moderate accuracy
• Different experiments need different numbers of reads: The extent of saturation can't easily be predicted from the depth
• subSeq replications might be necessary when the number of significant genes is < 500, rarely necessary when > 1000.