## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set( echo = TRUE, message = FALSE, warning = FALSE ) ## ----eval=FALSE--------------------------------------------------------------- # if (!require("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # # BiocManager::install("speckle") ## ----------------------------------------------------------------------------- library(speckle) library(SingleCellExperiment) library(CellBench) library(limma) library(ggplot2) library(scater) library(patchwork) library(edgeR) library(statmod) ## ----------------------------------------------------------------------------- sc_data <- load_sc_data() ## ----------------------------------------------------------------------------- commongenes1 <- rownames(sc_data$sc_dropseq)[rownames(sc_data$sc_dropseq) %in% rownames(sc_data$sc_celseq)] commongenes2 <- commongenes1[commongenes1 %in% rownames(sc_data$sc_10x)] sce_10x <- sc_data$sc_10x[commongenes2,] sce_celseq <- sc_data$sc_celseq[commongenes2,] sce_dropseq <- sc_data$sc_dropseq[commongenes2,] dim(sce_10x) dim(sce_celseq) dim(sce_dropseq) table(rownames(sce_10x) == rownames(sce_celseq)) table(rownames(sce_10x) == rownames(sce_dropseq)) ## ----------------------------------------------------------------------------- i.10x <- seq_len(ncol(sce_10x)) i.celseq <- seq_len(ncol(sce_celseq)) i.dropseq <- seq_len(ncol(sce_dropseq)) set.seed(10) boot.10x <- sample(i.10x, replace=TRUE) boot.celseq <- sample(i.celseq, replace=TRUE) boot.dropseq <- sample(i.dropseq, replace=TRUE) sce_10x_rep2 <- sce_10x[,boot.10x] sce_celseq_rep2 <- sce_celseq[,boot.celseq] sce_dropseq_rep2 <- sce_dropseq[,boot.dropseq] ## ----------------------------------------------------------------------------- sample <- rep(c("S1","S2","S3","S4","S5","S6"), c(ncol(sce_10x),ncol(sce_10x_rep2),ncol(sce_celseq), ncol(sce_celseq_rep2), ncol(sce_dropseq),ncol(sce_dropseq_rep2))) cluster <- c(sce_10x$cell_line,sce_10x_rep2$cell_line,sce_celseq$cell_line, sce_celseq_rep2$cell_line,sce_dropseq$cell_line, sce_dropseq_rep2$cell_line) group <- rep(c("10x","celseq","dropseq"), c(2*ncol(sce_10x),2*ncol(sce_celseq),2*ncol(sce_dropseq))) allcounts <- cbind(counts(sce_10x),counts(sce_10x_rep2), counts(sce_celseq), counts(sce_celseq_rep2), counts(sce_dropseq), counts(sce_dropseq_rep2)) sce_all <- SingleCellExperiment(assays = list(counts = allcounts)) sce_all$sample <- sample sce_all$group <- group sce_all$cluster <- cluster ## ----------------------------------------------------------------------------- sce_all <- scater::logNormCounts(sce_all) sce_all <- scater::runPCA(sce_all) sce_all <- scater::runUMAP(sce_all) ## ----fig.width=12, fig.height=6----------------------------------------------- pca1 <- scater::plotReducedDim(sce_all, dimred = "PCA", colour_by = "cluster") + ggtitle("Cell line") pca2 <- scater::plotReducedDim(sce_all, dimred = "PCA", colour_by = "group") + ggtitle("Technology") pca1 + pca2 ## ----fig.width=12, fig.height=6----------------------------------------------- umap1 <- scater::plotReducedDim(sce_all, dimred = "UMAP", colour_by = "cluster") + ggtitle("Cell line") umap2 <- scater::plotReducedDim(sce_all, dimred = "UMAP", colour_by = "group") + ggtitle("Technology") umap1 + umap2 ## ----------------------------------------------------------------------------- # Perform logit transformation propeller(sce_all) ## ----------------------------------------------------------------------------- # Perform arcsin square root transformation propeller(sce_all, transform="asin") ## ----------------------------------------------------------------------------- propeller(clusters=sce_all$cluster, sample=sce_all$sample, group=sce_all$group) ## ----fig.height=4,fig.width=7------------------------------------------------- plotCellTypeProps(sce_all) ## ----fig.height=5,fig.width=7------------------------------------------------- props <- getTransformedProps(sce_all$cluster, sce_all$sample, transform="logit") barplot(props$Proportions, col = c("orange","purple","dark green"),legend=TRUE, ylab="Proportions") ## ----fig.height=4,fig.width=10------------------------------------------------ par(mfrow=c(1,3)) for(i in seq(1,3,1)){ stripchart(props$Proportions[i,]~rep(c("10x","celseq","dropseq"),each=2), vertical=TRUE, pch=16, method="jitter", col = c("orange","purple","dark green"),cex=2, ylab="Proportions") title(rownames(props$Proportions)[i]) } ## ----------------------------------------------------------------------------- par(mfrow=c(1,1)) plotCellTypeMeanVar(props$Counts) plotCellTypePropsMeanVar(props$Counts) ## ----------------------------------------------------------------------------- names(props) props$TransformedProps ## ----------------------------------------------------------------------------- group <- rep(c("10x","celseq","dropseq"),each=2) pair <- rep(c(1,2),3) data.frame(group,pair) ## ----------------------------------------------------------------------------- design <- model.matrix(~ 0 + group + pair) design ## ----------------------------------------------------------------------------- propeller.anova(prop.list=props, design=design, coef = c(1,2,3), robust=TRUE, trend=FALSE, sort=TRUE) ## ----------------------------------------------------------------------------- design mycontr <- makeContrasts(group10x-groupdropseq, levels=design) propeller.ttest(props, design, contrasts = mycontr, robust=TRUE, trend=FALSE, sort=TRUE) ## ----------------------------------------------------------------------------- group dose <- rep(c(1,2,3), each=2) des.dose <- model.matrix(~dose) des.dose fit <- lmFit(props$TransformedProps,des.dose) fit <- eBayes(fit, robust=TRUE) topTable(fit) ## ----------------------------------------------------------------------------- fit.prop <- lmFit(props$Proportions,des.dose) fit.prop <- eBayes(fit.prop, robust=TRUE) topTable(fit.prop) ## ----fig.height=4,fig.width=10------------------------------------------------ fit.prop$coefficients par(mfrow=c(1,3)) for(i in seq(1,3,1)){ plot(dose, props$Proportions[i,], main = rownames(props$Proportions)[i], pch=16, cex=2, ylab="Proportions", cex.lab=1.5, cex.axis=1.5, cex.main=2) abline(a=fit.prop$coefficients[i,1], b=fit.prop$coefficients[i,2], col=4, lwd=2) } ## ----------------------------------------------------------------------------- des.tech <- model.matrix(~group) dupcor <- duplicateCorrelation(props$TransformedProps, design=des.tech, block=pair) dupcor ## ----------------------------------------------------------------------------- # Fitting the linear model accounting for pair as a random effect fit1 <- lmFit(props$TransformedProps, design=des.tech, block=pair, correlation=dupcor$consensus) fit1 <- eBayes(fit1) summary(decideTests(fit1)) # Differences between celseq vs 10x topTable(fit1,coef=2) # Differences between dropseq vs 10x topTable(fit1, coef=3) ## ----------------------------------------------------------------------------- topTable(fit1, coef=2:3) ## ----------------------------------------------------------------------------- data("pbmc_props") pbmc.props <- pbmc_props$proportions pbmc.sample.info <- pbmc_props$sample_info tot.cells <- pbmc_props$total_cells ## ----fig.width=10, fig.height=6----------------------------------------------- par(mfrow=c(1,1)) barplot(as.matrix(pbmc.props),col=ggplotColors(nrow(pbmc.props)), ylab="Cell type proportions", xlab="Samples") ## ----------------------------------------------------------------------------- prop.list <- convertDataToList(pbmc.props,data.type="proportions", transform="logit", scale.fac=tot.cells/20) ## ----------------------------------------------------------------------------- designAS <- model.matrix(~0+pbmc.sample.info$age + pbmc.sample.info$sex) colnames(designAS) <- c("old","young","MvsF") # Young vs old mycontr <- makeContrasts(young-old, levels=designAS) propeller.ttest(prop.list = prop.list,design = designAS, contrasts = mycontr, robust=TRUE,trend=FALSE,sort=TRUE) ## ----------------------------------------------------------------------------- group.immune <- paste(pbmc.sample.info$sex, pbmc.sample.info$age, sep=".") par(mfrow=c(1,2)) stripchart(as.numeric(pbmc.props["CD8.Naive",])~group.immune, vertical=TRUE, pch=c(8,16), method="jitter", col = c(ggplotColors(20)[20],"hotpink",4, "darkblue"),cex=2, ylab="Proportions", cex.axis=1.25, cex.lab=1.5, group.names=c("F.Old","F.Young","M.Old","M.Young")) title("CD8.Naive: Young Vs Old", cex.main=1.5, adj=0) stripchart(as.numeric(pbmc.props["CD16",])~group.immune, vertical=TRUE, pch=c(8,16), method="jitter", col = c(ggplotColors(20)[20],"hotpink",4, "darkblue"),cex=2, ylab="Proportions", cex.axis=1.25, cex.lab=1.5, group.names=c("F.Old","F.Young","M.Old","M.Young")) title("CD16: Young Vs Old", cex.main=1.5, adj=0) ## ----------------------------------------------------------------------------- sessionInfo()