inst/doc/conditionsVignette.R

In slingshot: Tools for ordering single-cell sequencing

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
require(slingshot)
require(RColorBrewer)
set.seed(1)

## ----dataSetup----------------------------------------------------------------
data('slingshotExample')
rd <- slingshotExample$rd
cl <- slingshotExample$cl
condition <- factor(rep(c('A','B'), length.out = nrow(rd)))
condition[110:140] <- 'A'
ls()

## ---- echo=FALSE--------------------------------------------------------------
plot(rd, asp = 1, pch = 16, col = brewer.pal(3,'Set1')[condition], las=1)
legend('topleft','(x,y)',legend = c('A','B'), title = 'Condition', pch=16, col = brewer.pal(3,'Set1')[1:2])

## -----------------------------------------------------------------------------
sds <- slingshot(rd, cl)

## ---- echo=FALSE--------------------------------------------------------------
plot(rd, asp = 1, pch = 16, col = brewer.pal(3,'Set1')[condition], las=1)
lines(sds, lwd=3)
legend('topleft','(x,y)',legend = c('A','B'), title = 'Condition', pch=16, col = brewer.pal(3,'Set1')[1:2])

## ---- fig.height=4, fig.width=8, echo = FALSE---------------------------------
n <- nrow(rd); L <- ncol(slingPseudotime(sds))
noise <- runif(n, -.1,.1)
plot(as.numeric(slingPseudotime(sds)), rep(1:L, each=n)+noise,pch=16, col = brewer.pal(9,'Set1')[condition], axes=FALSE, xlab='Pseudotime', ylab='Lineage', las=1)
axis(1); axis(2, at=1:L, las=1)

## ---- echo=FALSE--------------------------------------------------------------
# tA1 <- slingPseudotime(sds, na=FALSE)[condition=='A',1]
# wA1 <- slingCurveWeights(sds)[condition=='A',1]; wA1 <- wA1/sum(wA1)
# dA1 <- density(tA1, weights = wA1)
# tB1 <- slingPseudotime(sds, na=FALSE)[condition=='B',1]
# wB1 <- slingCurveWeights(sds)[condition=='B',1]; wB1 <- wB1/sum(wB1)
# dB1 <- density(tB1, weights = wB1)
# tA2 <- slingPseudotime(sds, na=FALSE)[condition=='A',2]
# wA2 <- slingCurveWeights(sds)[condition=='A',2]; wA2 <- wA2/sum(wA2)
# dA2 <- density(tA2, weights = wA2)
# tB2 <- slingPseudotime(sds, na=FALSE)[condition=='B',2]
# wB2 <- slingCurveWeights(sds)[condition=='B',2]; wB2 <- wB2/sum(wB2)
# dB2 <- density(tB2, weights = wB2)
# 
# plot(range(slingPseudotime(sds),na.rm=TRUE), c(1,2+7*max(c(dA2$y,dB2$y))), col='white', xlab='Pseudotime', ylab='Lineage', axes = FALSE, las=1)
# axis(1); axis(2, at=1:2)
# polygon(c(min(dA1$x),dA1$x,max(dA1$x)), 1+7*c(0,dA1$y,0), col=rgb(1,0,0,.5))
# polygon(c(min(dB1$x),dB1$x,max(dB1$x)), 1+7*c(0,dB1$y,0), col=rgb(0,0,1,.5))
# polygon(c(min(dA2$x),dA2$x,max(dA2$x)), 2+7*c(0,dA2$y,0), col=rgb(1,0,0,.5))
# polygon(c(min(dB2$x),dB2$x,max(dB2$x)), 2+7*c(0,dB2$y,0), col=rgb(0,0,1,.5))

layout(matrix(1:2,nrow=1))
boxplot(slingPseudotime(sds)[,1] ~ condition, col = brewer.pal(3,'Set1')[1:2], main = 'Lineage 1', xlab='Condition', ylab='Pseudotime', las=1, pch = 16)
boxplot(slingPseudotime(sds)[,2] ~ condition, col = brewer.pal(3,'Set1')[1:2], main = 'Lineage 2', xlab='Condition', ylab='Pseudotime', las=1, pch = 16)
layout(1)

## ---- eval=FALSE--------------------------------------------------------------
#  # Permutation test
#  t1 <- slingPseudotime(sds, na=FALSE)[,1]
#  w1 <- slingCurveWeights(sds)[,1]
#  d1 <- weighted.mean(t1[condition=='A'], w1[condition=='A']) -
#      weighted.mean(t1[condition=='B'], w1[condition=='B'])
#  dist1 <- replicate(1e4, {
#      condition.i <- sample(condition)
#      weighted.mean(t1[condition.i=='A'], w1[condition.i=='A']) -
#          weighted.mean(t1[condition.i=='B'], w1[condition.i=='B'])
#  })

## ---- echo=FALSE, fig.height=4, fig.width=9-----------------------------------
t1 <- slingPseudotime(sds, na=FALSE)[,1]
w1 <- slingCurveWeights(sds)[,1]
d1 <- weighted.mean(t1[condition=='A'], w1[condition=='A']) - 
    weighted.mean(t1[condition=='B'], w1[condition=='B'])
dist1 <- replicate(1e4, {
    condition.i <- sample(condition)
    weighted.mean(t1[condition.i=='A'], w1[condition.i=='A']) - 
        weighted.mean(t1[condition.i=='B'], w1[condition.i=='B'])
})
t2 <- slingPseudotime(sds, na=FALSE)[,2]
w2 <- slingCurveWeights(sds)[,2]
d2 <- weighted.mean(t2[condition=='A'], w2[condition=='A']) - 
    weighted.mean(t2[condition=='B'], w2[condition=='B'])
dist2 <- replicate(1e4, {
    condition.i <- sample(condition)
    weighted.mean(t2[condition.i=='A'], w2[condition.i=='A']) - 
        weighted.mean(t2[condition.i=='B'], w2[condition.i=='B'])
})

layout(matrix(1:2,nrow = 1))
hist(dist1, breaks=50, col='grey50', xlim = range(c(d1,dist1)), probability = TRUE, xlab = 'Difference of Weighted Means', main = 'Lineage 1 Permutation Test', las=1)
abline(v=d1, col=2, lwd=2)
legend('topright','(x,y)',legend = c('Null Distn.','Observed'), fill=c('grey50',NA), border=c(1,NA), lty=c(NA,1), lwd=c(NA,2), col=c(NA,2), merge = TRUE, seg.len = .6)

hist(dist2, breaks=50, col='grey50', xlim = range(c(d2,dist2)), probability = TRUE, xlab = 'Difference of Weighted Means', main = 'Lineage 2 Permutation Test', las=1)
abline(v=d2, col=2, lwd=2)
legend('topright','(x,y)',legend = c('Null Distn.','Observed'), fill=c('grey50',NA), border=c(1,NA), lty=c(NA,1), lwd=c(NA,2), col=c(NA,2), merge = TRUE, seg.len = .6)
layout(1)

## -----------------------------------------------------------------------------
paste0('Lineage 1 p-value: ', mean(abs(dist1) > abs(d1)))
paste0('Lineage 2 p-value: ', mean(abs(dist2) > abs(d2)))

## -----------------------------------------------------------------------------
# Kolmogorov-Smirnov test
ks.test(slingPseudotime(sds)[condition=='A',1], slingPseudotime(sds)[condition=='B',1])
ks.test(slingPseudotime(sds)[condition=='A',2], slingPseudotime(sds)[condition=='B',2])

## -----------------------------------------------------------------------------
pt <- slingPseudotime(sds, na=FALSE)
cw <- slingCurveWeights(sds)
assignment <- apply(cw, 1, which.max)
ptAs <- c() #assigned pseudotime
for(ii in 1:nrow(pt)) ptAs[ii] <- pt[ii,assignment[ii]] 

## ----echo=FALSE,fig.height=3.5------------------------------------------------
layout(matrix(1:2,nrow=1))
noise <- runif(n, -.05,.05)
plot(ptAs[assignment == 1], (as.numeric(condition)+noise)[assignment == 1], 
     col = brewer.pal(9,'Set1')[condition[assignment == 1]],
     xlab="Pseudotime", ylab="Condition", main="Lineage 1", pch=16, axes = FALSE)
axis(1); axis(2, at=seq_along(levels(condition)), labels = levels(condition), las = 1)

plot(ptAs[assignment == 2], (as.numeric(condition)+noise)[assignment == 2], 
     col = brewer.pal(9,'Set1')[condition[assignment == 2]],
     xlab="Pseudotime", ylab="Condition", main="Lineage 2", pch=16, axes = FALSE)
axis(1); axis(2, at=seq_along(levels(condition)), labels = levels(condition), las = 1)
layout(1)

## -----------------------------------------------------------------------------
# model for lineage 1: not significant
cond1 <- factor(condition[assignment == 1])
t1 <- ptAs[assignment == 1]
m1 <- glm(cond1 ~ t1, family = quasibinomial(link = "logit"))
summary(m1)

## -----------------------------------------------------------------------------
# model for lineage 2: significant
cond2 <- factor(condition[assignment == 2])
t2 <- ptAs[assignment == 2]
m2 <- glm(cond2 ~ t2, family = quasibinomial(link = "logit"))
summary(m2)

## -----------------------------------------------------------------------------
### note that logistic regression is monotone hence only allows for increasing or decreasing proportion of cells for a particular condition.
### hence, for complex trajectories, you could consider smoothing the pseudotime, i.e.
m1s <- mgcv::gam(cond1 ~ s(t1), family="quasibinomial")
summary(m1s)

## -----------------------------------------------------------------------------
m2s <- mgcv::gam(cond2 ~ s(t2), family="quasibinomial")
summary(m2s)

## ----echo=FALSE---------------------------------------------------------------
layout(matrix(1:2,nrow=1))
plot(m1s, shade = TRUE, main = "Lineage 1", 
     xlab="Pseudotime", ylab="Logit Prob(B)", scheme=0)
plot(m2s, shade = TRUE, main = "Lineage 2",
     xlab="Pseudotime", ylab="Logit Prob(B)")
layout(1)

## -----------------------------------------------------------------------------
require(mgcv, quietly = TRUE)
t1 <- pt[,1]
t2 <- pt[,2]
l1 <- as.numeric(assignment == 1)
l2 <- as.numeric(assignment == 2)
m <- gam(condition ~ s(t1, by=l1, id=1) + s(t2, by=l2, id=1),
         family = quasibinomial(link = "logit"))
summary(m)
### and then we're back to tradeSeq-like inference ...

## ---- echo=FALSE--------------------------------------------------------------
layout(matrix(1:2,nrow=1))
plot(m, select=1, shade=TRUE, main='Lineage 1',
     xlab="Pseudotime", ylab="Logit Prob(B)")
plot(m, select=2, shade=TRUE, main='Lineage 2',
     xlab="Pseudotime", ylab="Logit Prob(B)")
layout(1)

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