inst/doc/cydar.R

In cydar: Using Mass Cytometry for Differential Abundance Analyses

## ---- echo=FALSE, results="hide", message=FALSE-------------------------------
require(knitr)
opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE)

## ----setup, echo=FALSE, message=FALSE-----------------------------------------
library(cydar)
library(BiocParallel)
register(SerialParam())
set.seed(100)

## -----------------------------------------------------------------------------
ncells <- 20000
nda <- 200
nmarkers <- 31
down.pos <- 1.8
up.pos <- 1.2
conditions <- rep(c("A", "B"), each=3)
combined <- rbind(matrix(rnorm(ncells*nmarkers, 1.5, 0.6), ncol=nmarkers),
                  matrix(rnorm(nda*nmarkers, down.pos, 0.3), ncol=nmarkers),
                  matrix(rnorm(nda*nmarkers, up.pos, 0.3), ncol=nmarkers))
combined[,31] <- rnorm(nrow(combined), 1, 0.5) # last marker is a QC marker. 
combined <- 10^combined # raw intensity values                  
sample.id <- c(sample(length(conditions), ncells, replace=TRUE), 
               sample(which(conditions=="A"), nda, replace=TRUE), 
               sample(which(conditions=="B"), nda, replace=TRUE))
colnames(combined) <- paste0("Marker", seq_len(nmarkers))

## -----------------------------------------------------------------------------
library(ncdfFlow)
collected.exprs <- list()
for (i in seq_along(conditions)) {
    stuff <- list(combined[sample.id==i,,drop=FALSE])
    names(stuff) <- paste0("Sample", i)
    collected.exprs[[i]] <- poolCells(stuff)
}
names(collected.exprs) <- paste0("Sample", seq_along(conditions))
collected.exprs <- ncdfFlowSet(as(collected.exprs, "flowSet"))

## -----------------------------------------------------------------------------
pool.ff <- poolCells(collected.exprs)

## -----------------------------------------------------------------------------
library(flowCore)
trans <- estimateLogicle(pool.ff, colnames(pool.ff))
proc.ff <- transform(pool.ff, trans)

## -----------------------------------------------------------------------------
gate.31 <- outlierGate(proc.ff, "Marker31", type="upper")
gate.31
filter.31 <- filter(proc.ff, gate.31)
summary(filter.31@subSet)

## -----------------------------------------------------------------------------
processed.exprs <- transform(collected.exprs, trans)

## -----------------------------------------------------------------------------
processed.exprs <- Subset(processed.exprs, gate.31)

## -----------------------------------------------------------------------------
processed.exprs <- processed.exprs[,1:30]

## -----------------------------------------------------------------------------
cd <- prepareCellData(processed.exprs)

## -----------------------------------------------------------------------------
cd <- countCells(cd, tol=0.5)

## -----------------------------------------------------------------------------
head(assay(cd))

## -----------------------------------------------------------------------------
head(intensities(cd))

## -----------------------------------------------------------------------------
library(edgeR)
y <- DGEList(assay(cd), lib.size=cd$totals)

## -----------------------------------------------------------------------------
keep <- aveLogCPM(y) >= aveLogCPM(5, mean(cd$totals))
cd <- cd[keep,]
y <- y[keep,]

## -----------------------------------------------------------------------------
design <- model.matrix(~factor(conditions))
y <- estimateDisp(y, design)
fit <- glmQLFit(y, design, robust=TRUE)
res <- glmQLFTest(fit, coef=2)

## -----------------------------------------------------------------------------
qvals <- spatialFDR(intensities(cd), res$table$PValue)

## -----------------------------------------------------------------------------
is.sig <- qvals <= 0.05
summary(is.sig)

## -----------------------------------------------------------------------------
sig.coords <- intensities(cd)[is.sig,]
sig.res <- res$table[is.sig,]
coords <- prcomp(sig.coords)

## -----------------------------------------------------------------------------
plotSphereLogFC(coords$x[,1], coords$x[,2], sig.res$logFC)

## ----fig.width=10,fig.height=12-----------------------------------------------
par(mfrow=c(6,5), mar=c(2.1, 1.1, 3.1, 1.1))
limits <- intensityRanges(cd, p=0.05)
all.markers <- markernames(cd)
for (i in order(all.markers)) { 
    plotSphereIntensity(coords$x[,1], coords$x[,2], sig.coords[,i],
        irange=limits[,i], main=all.markers[i])
}

## -----------------------------------------------------------------------------
nonred <- findFirstSphere(intensities(cd), res$table$PValue)
summary(nonred)

## -----------------------------------------------------------------------------
all.coords <- prcomp(intensities(cd))
app <- interpretSpheres(cd, select=nonred, metrics=res$table, run=FALSE,
                        red.coords=all.coords$x[,1:2], red.highlight=is.sig) 
# Set run=TRUE if you want the app to run automatically.

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