## ----knitr-options, echo=FALSE, message=FALSE, warning=FALSE---------------
## To render an HTML version that works nicely with github and web pages, do:
## rmarkdown::render("vignettes/vignette.Rmd", "all")
library(knitr)
opts_chunk$set(fig.align = 'center', fig.width = 6, fig.height = 5, 
    error=FALSE, warning=FALSE, message=FALSE, dev = 'png')
library(ggplot2)
theme_set(theme_bw(12))

## --------------------------------------------------------------------------
library(scater)
data("sc_example_counts")
data("sc_example_cell_info")
example_sce <- SingleCellExperiment(
    assays = list(counts = sc_example_counts), 
    colData = sc_example_cell_info
)
example_sce

## ----calc-qc-metrics-------------------------------------------------------
example_sce <- calculateQCMetrics(example_sce)
colnames(colData(example_sce))
colnames(rowData(example_sce))

## ----calc-qc-metrics-multi-feature-controls--------------------------------
example_sce <- calculateQCMetrics(example_sce, 
    feature_controls = list(ERCC = 1:20, mito = 500:1000),
    cell_controls = list(empty = 1:5, damaged = 31:40))

all_col_qc <- colnames(colData(example_sce))
all_col_qc <- all_col_qc[grep("ERCC", all_col_qc)]

## ----plot-highest, fig.asp=2, fig.wide=TRUE--------------------------------
plotQC(example_sce, type = "highest-expression", exprs_values = "counts")

## ----plot-qc-exprs-freq-vs-mean-default------------------------------------
plotQC(example_sce, type = "exprs-freq-vs-mean")

## ----plot-pdata-pct-exprs-controls-----------------------------------------
plotColData(example_sce, x = "total_features_by_counts",
    y = "pct_counts_feature_control", colour = "Mutation_Status") +
    theme(legend.position = "top") +
    stat_smooth(method = "lm", se = FALSE, size = 2, fullrange = TRUE)

## ----plot-sceset-blocking--------------------------------------------------
plotScater(example_sce, block1 = "Mutation_Status", block2 = "Treatment",
     colour_by = "Cell_Cycle", nfeatures = 300, exprs_values = "counts")

## --------------------------------------------------------------------------
example_sce2 <- example_sce
example_sce2$plate_position <- paste0(
     rep(LETTERS[1:5], each = 8), 
     rep(formatC(1:8, width = 2, flag = "0"), 5)
)
plotPlatePosition(example_sce2, colour_by = "Gene_0001",
    by_exprs_values = "counts") 

## ----plot-fdata------------------------------------------------------------
plotRowData(example_sce, x = "n_cells_by_counts", y = "mean_counts")

## ----plot-pdata-2, fig.wide=TRUE, fig.asp=0.3------------------------------
p1 <- plotColData(example_sce, x = "total_counts", 
    y = "total_features_by_counts")
p2 <- plotColData(example_sce, x = "pct_counts_feature_control",
    y = "total_features_by_counts")
p3 <- plotColData(example_sce, x = "pct_counts_feature_control",
    y = "pct_counts_top_50_features")
multiplot(p1, p2, p3, cols = 3)

## ----highest-2, fig.wide=TRUE----------------------------------------------
p1 <- plotQC(example_sce[, !example_sce$is_cell_control],
    type = "highest-expression")
p2 <- plotQC(example_sce[, example_sce$is_cell_control],
    type = "highest-expression")
multiplot(p1, p2, cols = 2)

## --------------------------------------------------------------------------
example_sce <- example_sce[,1:40]

## --------------------------------------------------------------------------
filter(example_sce, Treatment == "treat1")

## --------------------------------------------------------------------------
keep.total <- example_sce$total_counts > 1e5
keep.n <- example_sce$total_features_by_counts > 500
filtered <- example_sce[,keep.total & keep.n]
dim(filtered)

## --------------------------------------------------------------------------
keep.total <- isOutlier(example_sce$total_counts, nmads=3, 
    type="lower", log=TRUE)
filtered <- example_sce[,keep.total]

## ----plot-pca-outlier------------------------------------------------------
example_sce <- runPCA(example_sce, use_coldata = TRUE,
    detect_outliers = TRUE)
plotReducedDim(example_sce, use_dimred="PCA_coldata")

## --------------------------------------------------------------------------
summary(example_sce$outlier)

## --------------------------------------------------------------------------
keep_feature <- nexprs(example_sce, byrow=TRUE) >= 4
example_sce <- example_sce[keep_feature,]
dim(example_sce)

## ----plot-qc-expl-factors-all----------------------------------------------
example_sce <- normalize(example_sce)
plotQC(example_sce, type = "expl")

## ----plot-qc-expl-variables-select-variables-------------------------------
plotQC(example_sce, type = "expl",
    variables = c("total_features_by_counts", "total_counts",
        "Mutation_Status", "Treatment", "Cell_Cycle"))

## ----plot-qc-pairs-pc------------------------------------------------------
plotQC(example_sce, type = "expl", method = "pairs", theme_size = 6)

## --------------------------------------------------------------------------
sizeFactors(example_sce) <- librarySizeFactors(example_sce)
summary(sizeFactors(example_sce))

## --------------------------------------------------------------------------
example_sce <- normalize(example_sce)

## --------------------------------------------------------------------------
library(limma)
batch <- rep(1:2, each=20)
corrected <- removeBatchEffect(logcounts(example_sce), block=batch)
assay(example_sce, "corrected_logcounts") <- corrected

## --------------------------------------------------------------------------
sessionInfo()