1.1 An Example Dataset

The first two parts of this workshop will analyze a small set of cells lying along a developmental trajectory. We will start from raw data objects obtained from a standard transcriptome alignment pipeline. Raw data and important reference data can be loaded directly from the workshop package.

library(bioc2016singlecell)

## Load Example Data
data(ws_input)
ls()

##  [1] "batch"          "bio"            "breaks"         "cl"            
##  [5] "cl3"            "cl7"            "cm"             "cont_levels"   
##  [9] "contrMat"       "counts"         "de"             "defaultMar"    
## [13] "fit"            "genesDendro"    "genesF"         "genesOneAll"   
## [17] "genesPairs"     "hk"             "manual"         "mast_res"      
## [21] "norm"           "norm_logcounts" "pca"            "plotCMar"      
## [25] "qc"             "res1"           "res2"           "res3"          
## [29] "rs"             "sca"            "se"             "wh_rm"

The counts data frame contains feature-level read counts from tophat alignments of 96 single-cell libraries to the mm10 reference genome (Trapnell, Pachter, and Salzberg 2009). qc contains library alignment metrics obtained from Picard.

colnames(qc)

##  [1] "NREADS"               "NALIGNED"             "RALIGN"              
##  [4] "TOTAL_DUP"            "PRIMER"               "PCT_RIBOSOMAL_BASES" 
##  [7] "PCT_CODING_BASES"     "PCT_UTR_BASES"        "PCT_INTRONIC_BASES"  
## [10] "PCT_INTERGENIC_BASES" "PCT_MRNA_BASES"       "MEDIAN_CV_COVERAGE"  
## [13] "MEDIAN_5PRIME_BIAS"   "MEDIAN_3PRIME_BIAS"

de and hk are positive and negative control gene sets derived from population studies. batch and bio are factors labeling batch and time point respectively. Consider the joint distribution of these factors:

## Joint distribution of batches and biological conditions (time after induction)
table(batch,bio)

##       bio
## batch  time0 time1 time2 time3 time4 time5
##   Y01     10     0     0     0     0     0
##   Y02      0     0     0     6     0     0
##   Y03      0     0     0     8     0     0
##   Y04      8     0     0     0     0     0
##   Y05A     0    10     0     0     0     0
##   Y05B     0     9     0     0     0     0
##   Y06      0     0     5     0     0     0
##   Y07A     0     0     9     0     0     0
##   Y07B     0     0    10     0     0     0
##   Y08      0     6     0     0     0     0
##   Y10      0     0     0     0     2     0
##   Y11      0     0     0     0     4     0
##   Y12A     0     0     0     0     0     2
##   Y12B     0     0     0     0     0     2
##   Y13      0     0     0     0     0     5

Notice that each time-point is composed of multiple batches. This nested design is common among scRNA-Seq studies due to limitations on the number of cells that can be harvested/sequenced concurrently.

1.2 Technical Variability and Batch Effects

We will first visualize cell level quality readouts, starting with the ratio of reads aligned to the genome:

# Color scheme
cc <- c(brewer.pal(9, "Set1"), brewer.pal(8, "Set2"),brewer.pal(9,"Set3"))

# Barplot of read proportion mapping to the genome
barplot(qc$RALIGN[order(batch)], col=cc[batch][order(batch)], border=cc[batch][order(batch)], main="Percentage of mapped reads, colored by batch")
legend("bottomleft", legend=levels(batch), fill=cc,cex=0.4)

We see that there aren’t important differences between batches, while there are a few cells of particularly low alignment efficiency. These cells can be removed via filtering procedures. We can alternatively consider the number of reads in each library:

# Barplot of total read number
barplot(qc$NREADS[order(batch)], col=cc[batch][order(batch)], border=cc[batch][order(batch)], main="Total number of reads, colored by batch")
legend("topright", legend=levels(batch), fill=cc, cex=0.4)

We see that coverage varies substantially between batches, and though some of this technical variability may be addressed with filtering, we will still need to carefully normalize the data to make batches comparable.

It can be very helpful to visualize distributions of single metrics, but it’s important to note that these metrics are often correlated. Sometimes it may be more useful to consider principal components of the quality matrix, identifying principal factors of library variation:

qpc = prcomp(qc,center = TRUE,scale. = TRUE)
barplot(cumsum(qpc$sdev^2)/sum(qpc$sdev^2), border="gray", xlab="PC", ylab="proportion of variance", main="Quality PCA")

Even though 14 quality metrics have been quantified, PCA shows us that only a small number of PCs are needed to described a majority of the variance (e.g. 3 to explain 74%). We will now visualize the distribution of the first PC in the context of batch:

barplot(qpc$x[,1][order(batch)], col=cc[batch][order(batch)], border=cc[batch][order(batch)], main="Quality PC1, colored by batch")
legend("bottomleft", legend=levels(batch), fill=cc, cex=0.8)

This PC appears to correlate with batch (particularly Y04), but it also highlights 5 outlier libraries that are significantly different from the rest. While we won’t discuss the application of these PCs to filtering in this workshop, we will show how they can be applied to normalization.

1.3 Drop-out Characteristics

Next we consider a uniquely single-cell problem: drop-outs. hk contains a list of genes that are believed to be ubiquitously and uniformly expressed across the target tissue. Because we assume these genes are expressed in all cells, we can label all zeroes as drop-out events. Below we model detection failures as a logistic function of mean expression, in line with the standard logistic model for drop-outs employed by the field.

# Mean log10(x+1) expression
mu_obs = rowMeans(log10(counts[hk,]+1))

# Drop-outs
drop_outs = counts[hk,] == 0

# Logistic Regression Model of Failure
ref.glms = list()
for (si in 1:dim(drop_outs)[2]){
  fit = glm(cbind(drop_outs[,si],1 - drop_outs[,si]) ~ mu_obs,family=binomial(logit))
  ref.glms[[si]] = fit$coefficients
}

The list ref.glm contains the intercept and slope of each fit. We can now visualize the fit curves and the corresponding Area Under the Curves (AUCs):

par(mfrow=c(1,2))

# Plot Failure Curves and Calculate AUC
plot(NULL, main = "False Negative Rate Curves",ylim = c(0,1),xlim = c(0,6), ylab = "Failure Probability", xlab = "Mean log10 Expression")
x = (0:60)/10
AUC = NULL
for(si in 1:ncol(counts)){
  y = 1/(exp(-ref.glms[[si]][1] - ref.glms[[si]][2] * x) + 1)
  AUC[si] = sum(y)/10
  lines(x, 1/(exp(-ref.glms[[si]][1] - ref.glms[[si]][2] * x) + 1), type = 'l', lwd = 2, col = cc[batch][si])
}

# Barplot of FNR AUC
barplot(AUC[order(batch)], col=cc[batch][order(batch)], border=cc[batch][order(batch)], main="FNR AUC, colored by batch")
legend("topleft", legend=levels(batch), fill=cc, cex=0.4)

You may have noticed that the 5 outlier cells seen here are the same as the 5 outliers highlighted in the PCA of quality metrics above. Drop-out summaries such as the AUC can very useful for assessing single-cell library quality.

1.4 The scone Workflow

So far we have only described potential problems with our data set. Now we will take steps to address them! The basic qc and normalization pipeline we will use today will allow us to:

• Filter out poor libraries using the metric_sample_filter function.
• Run and score many different normalization workflows (different combinations of normalization modules) using the main scone function.
• Browse top-ranked methods and visualize trade-offs with the biplot_colored function.

The SCONE’s normalization workflow template is composed of 3 modules:

1. Data imputation: replacing zero-abundance values with expected values under a drop-out model. We will only consider identity imputation (no imputation) in this example.
2. Scaling or quantile normalization: either i) normalization that scales each sample’s transcriptome abundances by a single factor or ii) more complex offsets that match quantiles across samples. Examples: TMM or DESeq scaling factors, upper quartile normalization, or full-quantile normalization.
3. Regression-based approaches for removing unwanted correlated variation from the data. Examples: RUVg (Risso et al. 2014) or regression on Quality Principal Components described above.

mfilt_report <- metric_sample_filter(expr = counts,
                                     nreads = qc$NREADS,ralign = qc$RALIGN,
                                     suff_nreads = 10^5,
                                     suff_ralign = 90,
                                     pos_controls = hk,
                                     zcut = 3,mixture = FALSE, plot = TRUE)

2.1 On Threshold Selection

Let’s take a closer look at the computation behind the ralign filter. In selecting the threshold value 90, metric_sample_filter is taking 4 values into account:

1. hard_ralign, the default “hard” threshold at 15 - rather forgiving…
2. 3 (zcut) times the standard deviation below the mean ralign value.
3. 3 (zcut) times the MAD below the median ralign value.
4. suff_ralign, the sufficient threshold set to 90.

hist(qc$RALIGN, breaks = 0:100)
# Hard threshold
abline(v = 15, col = "yellow", lwd = 2) 
# 3 (zcut) standard deviations below the mean ralign value
abline(v = mean(qc$RALIGN) - 3*sd(qc$RALIGN), col = "green", lwd = 2) 
# 3 (zcut) MADs below the median ralign value
abline(v = median(qc$RALIGN) - 3*mad(qc$RALIGN), col = "red", lwd = 2)
# Sufficient threshold
abline(v = 90, col = "grey", lwd = 2)
# Final threshold is the minimum of 1) the sufficient threshold and 2) the max of all others
thresh = min(90,max(c(15,mean(qc$RALIGN) - 3*sd(qc$RALIGN),median(qc$RALIGN) - 3*mad(qc$RALIGN))))
abline(v = thresh, col = "blue", lwd = 2, lty = 2)

legend("topleft",legend = c("Hard","SD","MAD","Sufficient","Final"),lwd = 2, col = c("yellow","green","red","grey","blue"),lty = c(1,1,1,1,2), cex = .5)

We see here that the 3rd threshold exceeds the sufficient threshold, and so metric_sample_filter settles for the sufficient threshold. Note that if mixture=TRUE an additional criterion is included: bi-modal metric distributions are fit with to a two component mixture model, and a threshold is defined with respect to the mean and standard deviation of the “better” component. As metric_sample_filter relies on a maximum of candidate thresholds, we recommend users treat this function as a liberal filter.

2.2 Applying the sample filter

With metric_sample_filter output in hand, filtering out the poor samples is fairly straightforward:

# Which thresholds are missing? (breadth)
is_na_filt = unlist(lapply(is.na(mfilt_report),any)) 

# Identify samples failing any threshold
m_sampfilter = !apply(simplify2array(mfilt_report[!is_na_filt]),1,any)

# Filter Samples
fcounts = counts[,m_sampfilter]
fqc = qc[m_sampfilter,]
fbatch = batch[m_sampfilter]
fbio = bio[m_sampfilter]

# Simple gene filter
filterCount <- function(counts, nRead=5, nCell=5){
  filter <- apply(counts, 1, function(x) length(x[x>=nRead])>=nCell)
  return(filter)
}
genefilter <- filterCount(fcounts)

# Filter genes
fcounts = fcounts[genefilter,]
fhk = hk[hk %in% rownames(fcounts)]
fde = de[de %in% rownames(fcounts)]

3.1 Selecting Workflow Parameters with run=FALSE

The arguments to scone allow for a lot of flexibility, but sometimes you will need to run very specific combinations of modules. For this purpose, scone can be run in run=FALSE mode, returning only a data frame of workflows to be performed.

params <- scone(expr = as.matrix(fcounts),scaling = c(none = identity,deseq = DESEQ_FN, tmm = TMM_FN, uqp = UQ_FN_POS, fq = FQT_FN),
                ruv_negcon = fhk, k_ruv = 3,
                qc = as.matrix(fqc), k_qc = 3,
                bio = fbio,adjust_bio = "yes",
                batch = fbatch,adjust_batch = "yes",
                run = FALSE)
head(params)

##                                   imputation_method scaling_method
## none,none,no_uv,no_bio,no_batch                none           none
## none,deseq,no_uv,no_bio,no_batch               none          deseq
## none,tmm,no_uv,no_bio,no_batch                 none            tmm
## none,uqp,no_uv,no_bio,no_batch                 none            uqp
## none,fq,no_uv,no_bio,no_batch                  none             fq
## none,none,ruv_k=1,no_bio,no_batch              none           none
##                                   uv_factors adjust_biology adjust_batch
## none,none,no_uv,no_bio,no_batch        no_uv         no_bio     no_batch
## none,deseq,no_uv,no_bio,no_batch       no_uv         no_bio     no_batch
## none,tmm,no_uv,no_bio,no_batch         no_uv         no_bio     no_batch
## none,uqp,no_uv,no_bio,no_batch         no_uv         no_bio     no_batch
## none,fq,no_uv,no_bio,no_batch          no_uv         no_bio     no_batch
## none,none,ruv_k=1,no_bio,no_batch    ruv_k=1         no_bio     no_batch

In the call above, we have set the following parameters:

• scaling = list(…). This argument contains a list of scaling normalization functions that will be applied, including the identity (no-op), DESeq scaling, TMM normalization, scaling by the upper quartile of positive counts, and full-quantile normalization.
• ruv_negcon = fhk. A list of genes to be used as negative controls for RUVg normalization.
• k_ruv = 3. The maximum number of RUVg factors to consider.
• k_qc = 3. The maximum number of quality PCs (QPCs) to be included in a linear model, analogous to RUVg normalization. The qc argument must be provided.
• adjust_bio = “yes.” Time points will be included in RUVg or QPC normalization and restored to the residual of the regression. The bio argument must be provided.
• adjust_batch = “yes.” Batch will be included in RUVg or QPC normalization. The batch argument must be provided.

These arguments translate to the following set of options:

apply(params,2,unique)

## $imputation_method
## [1] "none"
## 
## $scaling_method
## [1] "none"  "deseq" "tmm"   "uqp"   "fq"   
## 
## $uv_factors
## [1] "no_uv"   "ruv_k=1" "ruv_k=2" "ruv_k=3" "qc_k=1"  "qc_k=2"  "qc_k=3" 
## 
## $adjust_biology
## [1] "no_bio" "bio"   
## 
## $adjust_batch
## [1] "no_batch" "batch"

While adjusting for biology prior to downstream analysis is generally problematic, it can be particularly dangerous in the single-cell context: there is significantly more sample-intrinsic biological variability in scRNA-Seq data, variability that cannot be captured by exposure. In our nested design scenario we would only consider adjustments for biology when batch effects are included. We can produce an updated params data frame reflecting this choice:

is_screened = (params$adjust_biology == "bio") & (params$adjust_batch != "batch")

params = params[!is_screened,]

3.2 Calling scone with run=TRUE

Now that we have selected our workflows, we can run scone in run=TRUE mode. This mode offers a few additional arguments, including an optional params argument to pass any results from the run=FALSE mode. In order to understand these arguments, we must first understand the 8 metrics used to evaluate each normalization. The first 6 metrics rely on a reduction of the normalized data down to 3 dimensions via PCA (default). Each metric is taken to have a positive (higher is better) or negative (lower is better) signature.

• BIO_SIL. The average silhouette width of clusters defined by bio, defined with respect to a Euclidean distance metric over the first 3 expression PCs. Positive signature.
• BATCH_SIL. The average silhouette width of clusters defined by batch, defined with respect to a Euclidean distance metric over the first 3 expression PCs. Negative signature.
• PAM_SIL. The maximum average silhouette width of clusters defined by PAM clustering, defined with respect to a Euclidean distance metric over the first 3 expression PCs. Positive signature.
• EXP_QC_COR. Maximum squared Spearman correlation between first 3 expression PCs and first k_qc QPCs. Negative signature.
• EXP_UV_COR. Maximum squared Spearman correlation between first 3 expression PCs and first 3 PCs of the negative control (specified by eval_negcon or ruv_negcon by default) sub-matrix of the original (raw) data. Negative signature.
• EXP_WV_COR. Maximum squared Spearman correlation between first 3 expression PCs and first 3 PCs of the positive control (specified by eval_poscon) sub-matrix of the original (raw) data. Positive signature.
• RLE_MED. The mean squared median Relative Log Expression (RLE). Negative signature.
• RLE_IQR. The mean inter-quartile range (IQR) of the RLE. Negative signature.

res <- scone(expr = as.matrix(fcounts),scaling = c(none = identity, deseq = DESEQ_FN, tmm = TMM_FN, uqp = UQ_FN_POS, fq = FQT_FN),
             ruv_negcon = fhk, k_ruv = 3,
             qc = as.matrix(fqc), k_qc = 3,
             bio = fbio,adjust_bio = "yes",
             batch = fbatch,adjust_batch = "yes",
             run = TRUE,params = params,
             eval_poscon = fde, eval_kclust = 2:3, conditional_pam = TRUE)

In the call above, we have set the following parameters:

• eval_poscon = fde. A list of genes to be used as positive controls for evaluation.
• eval_kclust = 2:3. For PAM_SIL, range of k (# of clusters) to use when computing maximum average silhouette width of PAM clusterings.
• conditional_pam = TRUE. For PAM_SIL, apply separate PAM clusterings to each time point rather than across all time points. Average is weighted by time point group size.

Running the command will take a couple minutes. The output will contain a list of four elements:

names(res)

## [1] "normalized_data" "metrics"         "scores"          "params"

normalized_data contains a list of normalized expression data (log-scale); each expression matrix is named according to the same convention as seen in the params argument.

metrics contains the 8 raw metrics for each normalization. scores contains metrics multiplied by their signature - or “scores” - as well as a 9th column that contains the mean score for that normalization. Normalization workflows in normalized_data,metrics, and scores are sorted in decreasing order by mean score.

head(res$scores)

##                                       BIO_SIL BATCH_SIL   PAM_SIL
## none,fq,ruv_k=1,no_bio,no_batch    0.08653280 0.1508225 0.3741884
## none,fq,ruv_k=2,no_bio,no_batch    0.08126123 0.1558709 0.3716837
## none,fq,ruv_k=3,no_bio,no_batch    0.09865276 0.1415136 0.3629964
## none,deseq,ruv_k=2,no_bio,no_batch 0.09721160 0.1722529 0.3983736
## none,deseq,ruv_k=1,no_bio,no_batch 0.08828919 0.1640219 0.4215769
## none,fq,qc_k=1,no_bio,no_batch     0.08318459 0.1581687 0.4326760
##                                    EXP_QC_COR  EXP_UV_COR EXP_WV_COR
## none,fq,ruv_k=1,no_bio,no_batch    -0.2132825 -0.09864541  0.5405376
## none,fq,ruv_k=2,no_bio,no_batch    -0.2012192 -0.08994268  0.5313734
## none,fq,ruv_k=3,no_bio,no_batch    -0.2300302 -0.06228916  0.5779535
## none,deseq,ruv_k=2,no_bio,no_batch -0.2984754 -0.16153256  0.5042606
## none,deseq,ruv_k=1,no_bio,no_batch -0.3437892 -0.16409015  0.4916737
## none,fq,qc_k=1,no_bio,no_batch     -0.2336496 -0.28368051  0.5126642
##                                         RLE_MED    RLE_IQR  mean_score
## none,fq,ruv_k=1,no_bio,no_batch    -0.002169450 -0.9977360 -0.01996901
## none,fq,ruv_k=2,no_bio,no_batch    -0.002856720 -1.0625194 -0.02704359
## none,fq,ruv_k=3,no_bio,no_batch    -0.002686065 -1.1112346 -0.02814046
## none,deseq,ruv_k=2,no_bio,no_batch -0.007600576 -0.9843497 -0.03498245
## none,deseq,ruv_k=1,no_bio,no_batch -0.007537978 -0.9375725 -0.03592851
## none,fq,qc_k=1,no_bio,no_batch     -0.003546027 -0.9857528 -0.03999193