## ---- echo=FALSE, results="hide"---------------------------------------------- knitr::opts_chunk$set(error=FALSE, warning=FALSE, message=FALSE) library(BiocStyle) set.seed(10918) ## ----quickstart-load-data, message=FALSE, warning=FALSE----------------------- library(scRNAseq) example_sce <- ZeiselBrainData() example_sce ## ----quickstart-add-exprs, results='hide'------------------------------------- str(counts(example_sce)) ## ----------------------------------------------------------------------------- example_sce$whee <- sample(LETTERS, ncol(example_sce), replace=TRUE) colData(example_sce) rowData(example_sce)$stuff <- runif(nrow(example_sce)) rowData(example_sce) ## ----------------------------------------------------------------------------- library(scater) per.cell <- perCellQCMetrics(example_sce, subsets=list(Mito=grep("mt-", rownames(example_sce)))) summary(per.cell$sum) summary(per.cell$detected) summary(per.cell$subsets_Mito_percent) ## ----------------------------------------------------------------------------- colData(example_sce) <- cbind(colData(example_sce), per.cell) ## ----------------------------------------------------------------------------- plotColData(example_sce, x = "sum", y="detected", colour_by="tissue") ## ----plot-pdata-pct-exprs-controls-------------------------------------------- plotColData(example_sce, x = "sum", y="subsets_Mito_percent", other_fields="tissue") + facet_wrap(~tissue) ## ----------------------------------------------------------------------------- keep.total <- example_sce$sum > 1e5 keep.n <- example_sce$detected > 500 filtered <- example_sce[,keep.total & keep.n] dim(filtered) ## ----------------------------------------------------------------------------- keep.total <- isOutlier(per.cell$sum, type="lower", log=TRUE) filtered <- example_sce[,keep.total] ## ----------------------------------------------------------------------------- qc.stats <- quickPerCellQC(per.cell, percent_subsets="subsets_Mito_percent") colSums(as.matrix(qc.stats)) filtered <- example_sce[,!qc.stats$discard] ## ----------------------------------------------------------------------------- # Pretending that the first 10 cells are empty wells, for demonstration. per.feat <- perFeatureQCMetrics(example_sce, subsets=list(Empty=1:10)) summary(per.feat$mean) summary(per.feat$detected) summary(per.feat$subsets_Empty_ratio) ## ----------------------------------------------------------------------------- ave <- calculateAverage(example_sce) summary(ave) ## ----------------------------------------------------------------------------- summary(nexprs(example_sce, byrow=TRUE)) ## ----plot-highest, fig.asp=1, fig.wide=TRUE, eval=.Platform$OS.type!="windows"---- plotHighestExprs(example_sce, exprs_values = "counts") ## ----filter-no-exprs---------------------------------------------------------- keep_feature <- nexprs(example_sce, byrow=TRUE) > 0 example_sce <- example_sce[keep_feature,] dim(example_sce) ## ----------------------------------------------------------------------------- example_sce <- logNormCounts(example_sce) # see below. vars <- getVarianceExplained(example_sce, variables=c("tissue", "total mRNA mol", "sex", "age")) head(vars) ## ----------------------------------------------------------------------------- plotExplanatoryVariables(vars) ## ----------------------------------------------------------------------------- example_sce <- logNormCounts(example_sce) assayNames(example_sce) ## ----------------------------------------------------------------------------- summary(librarySizeFactors(example_sce)) ## ----------------------------------------------------------------------------- cpm(example_sce) <- calculateCPM(example_sce) ## ----------------------------------------------------------------------------- assay(example_sce, "normed") <- normalizeCounts(example_sce, size_factors=runif(ncol(example_sce)), pseudo_count=1.5) ## ----------------------------------------------------------------------------- agg_sce <- aggregateAcrossCells(example_sce, ids=example_sce$level1class) head(assay(agg_sce)) colData(agg_sce)[,c("ids", "ncells")] ## ----------------------------------------------------------------------------- agg_sce <- aggregateAcrossCells(example_sce, ids=colData(example_sce)[,c("level1class", "tissue")]) head(assay(agg_sce)) colData(agg_sce)[,c("level1class", "tissue", "ncells")] ## ----------------------------------------------------------------------------- agg_feat <- sumCountsAcrossFeatures(example_sce, ids=list(GeneSet1=1:10, GeneSet2=11:50, GeneSet3=1:100), average=TRUE, exprs_values="logcounts") agg_feat[,1:10] ## ----plot-expression, fig.wide=TRUE------------------------------------------- plotExpression(example_sce, rownames(example_sce)[1:6], x = "level1class") ## ----plot-expression-scatter-------------------------------------------------- plotExpression(example_sce, rownames(example_sce)[1:6], x = rownames(example_sce)[10]) ## ----plot-expression-col------------------------------------------------------ plotExpression(example_sce, rownames(example_sce)[1:6], x = "level1class", colour_by="tissue") ## ----plot-expression-many----------------------------------------------------- plotExpression(example_sce, rownames(example_sce)[1:6]) ## ----------------------------------------------------------------------------- example_sce <- runPCA(example_sce) str(reducedDim(example_sce, "PCA")) ## ----------------------------------------------------------------------------- example_sce <- runPCA(example_sce, name="PCA2", subset_row=rownames(example_sce)[1:1000], ncomponents=25) str(reducedDim(example_sce, "PCA2")) ## ----plot-tsne-1comp-colby-sizeby-exprs--------------------------------------- # Perplexity of 10 just chosen here arbitrarily. set.seed(1000) example_sce <- runTSNE(example_sce, perplexity=10) head(reducedDim(example_sce, "TSNE")) ## ----plot-tsne-from-pca------------------------------------------------------- set.seed(1000) example_sce <- runTSNE(example_sce, perplexity=50, dimred="PCA", n_dimred=10) head(reducedDim(example_sce, "TSNE")) ## ----------------------------------------------------------------------------- example_sce <- runUMAP(example_sce) head(reducedDim(example_sce, "UMAP")) ## ----plot-reduceddim-4comp-colby-shapeby-------------------------------------- plotReducedDim(example_sce, dimred = "PCA", colour_by = "level1class") ## ----plot-pca-4comp-colby-sizeby-exprs---------------------------------------- plotTSNE(example_sce, colour_by = "Snap25") ## ----plot-pca-default--------------------------------------------------------- plotPCA(example_sce, colour_by="Mog") ## ----plot-pca-4comp-colby-shapeby--------------------------------------------- example_sce <- runPCA(example_sce, ncomponents=20) plotPCA(example_sce, ncomponents = 4, colour_by = "level1class") ## ---- fig.wide=TRUE----------------------------------------------------------- ggcells(example_sce, mapping=aes(x=level1class, y=Snap25)) + geom_boxplot() + facet_wrap(~tissue) ## ----------------------------------------------------------------------------- ggcells(example_sce, mapping=aes(x=TSNE.1, y=TSNE.2, colour=Snap25)) + geom_point() + stat_density_2d() + facet_wrap(~tissue) + scale_colour_distiller(direction=1) ## ----------------------------------------------------------------------------- ggcells(example_sce, mapping=aes(x=sizeFactor, y=Actb)) + geom_point() + geom_smooth() ## ----------------------------------------------------------------------------- colnames(example_sce) <- make.names(colnames(example_sce)) ggfeatures(example_sce, mapping=aes(x=featureType, y=X1772062111_E06)) + geom_violin() ## ----------------------------------------------------------------------------- sessionInfo()