PBMC examples

Initialization

library(pagoda2)
library(Matrix)
library(ggplot2)
library(cowplot)
library(dplyr)
library(ggrepel)
require(parallel)
library(ggrastr)
library(conos)

Load Rahul’s 2019 integration pbmc dataset, as described here: https://satijalab.org/signac/articles/pbmc_vignette.html

#so <- readRDS("pbmc_10k_v3.rds")
load("seurat.RData") # cd, cd.var,meta, umap

ann <- as.factor(setNames(meta$celltype,rownames(meta)))

p2 <- Pagoda2$new(cd,log.scale=T, n.cores=30)
#p2 <- con$samples$pbmc60

p2$adjustVariance(plot=T,gam.k=10)

p2$calculatePcaReduction(nPcs=50,n.odgenes=3e3)
p2$makeKnnGraph(k=40,type='PCA',center=T,distance='cosine');

p2$getEmbedding(type='PCA',embeddingType='tSNE',n.cores=30)

a1 <- p2$plotEmbedding(type='PCA',embeddingType='tSNE',groups=ann,alpha=0.4, font.size=c(3,4))
#pdf(file='pmbc10k_tsne.pdf',width=3,height=3); print(a1); dev.off()
a1

p2$getKnnClusters(type='PCA',method=leiden.community,r=2,name='leiden')

p2$plotEmbedding(type='PCA',embeddingType='tSNE',clusterType='leiden',title='clusters (tSNE)')

de <- p2$getDifferentialGenes(type='PCA',verbose=T,groups=as.factor(ann),append.auc=T,upregulated.only = T)

Variance ranking and expression patterns

x <- p2$misc$varinfo
x$rank <- rank(-1*x$res)
x <- x[order(x$lp,decreasing=T),]
x <- x[order(x$res,decreasing=T),]
x <- x[order(x$lp,decreasing=F),]
x <- x[x$rank<5e3,]
ggplot(x,aes(x=rank,y=res)) + geom_point(color=adjustcolor(2,alpha=0.4),size=0.2) +theme_bw() + xlab("overdispersed gene rank") + ylab("residual variance")

gns <- c("S100A9",'HES4',"GAS6")
#gns <- rownames(x)[c(1,2000,2003)]
po <- lapply(gns,function(g) {
  p2$plotEmbedding(gene=g,size=0.2,alpha=0.7,raster=T,raster.dpi=300)+
  #conos::embeddingPlot(p2$embeddings$PCA$tSNE,colors=p2$counts[,g],size=0.2,alpha=0.7,raster=T,raster.height = 3,raster.width=3)+
  annotate('text',x=-Inf,y=Inf,vjust=1.2,hjust=0,label=g,size=6)+
  theme(panel.border = element_rect(color = 1, size=0.2,linetype=1),axis.line=element_blank()) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
})
pp <- plot_grid(plotlist=c(list(a1),po),nrow=1)
pp

pdf(file='pmbc10k_expr.pdf',width=6,height=2); print(plot_grid(plotlist = po,nrow=1)); dev.off()

null device 
          1 

y <- de[['pDC']]
y <- y[order(y$AUC,decreasing = T),]
y$rank <- match(rownames(y),rownames(x))
head(y[,c(1,2,6,10)],100)

x <- p2$misc$varinfo
x$rank <- rank(-1*x$res)
#x$rank <- rank(x$lp)
x <- x[order(x$lp,decreasing=F),]
x <- x[x$rank<5e3,]
x$name <- ''
for(g in gns) x[g,'name'] <- g
ggplot(x,aes(x=rank,y=res,label=name)) + geom_point(color=adjustcolor(2,alpha=0.4),size=0.2) +theme_bw() + xlab("overdispersed gene rank") + ylab("residual variance")+geom_text_repel()

x <- p2$misc$varinfo
x$rank <- rank(-1*x$res)
x$rank <- rank(x$lp)
x <- x[order(x$lp,decreasing=F),]
x <- x[x$rank<5e3,]
x$name <- ''
x$signif <- x$lpa<log(0.05)
for(g in gns) x[g,'name'] <- g
pv <- ggplot(x,aes(x=rank,y=res,label=name,color=signif)) + geom_point(size=0.1) +theme_bw() + xlab("gene overdispersion rank") + ylab("residual variance")+geom_text_repel(color=1) + geom_vline(xintercept = sum(x$lpa<log(0.05)),linetype=2,color=8) + guides(color='none') + scale_color_manual(values=adjustcolor(c('gray40','red'),alpha=0.15))
pv

gns <- c("CD14",'FYB1',"CD8A")

gns <- rownames(x)[c(1:16)+1016] # IGHG4
po <- lapply(gns,function(g) {
  p2$plotEmbedding(gene=g,size=0.1,alpha=0.6,raster=T,raster.dpi=300) +
  #conos::embeddingPlot(p2$embeddings$PCA$tSNE,colors=p2$counts[,g],size=0.1,alpha=0.6,raster=T,raster.dpi=300,theme=theme_bw())+
  annotate('text',x=-Inf,y=Inf,vjust=1.2,hjust=0,label=g,size=6)
})
pp <- plot_grid(plotlist=po,nrow=4)
pp

Use a simpler linear scale varnorm:

p2l <- Pagoda2$new(cd,log.scale=F, n.cores=30)

p2l$adjustVariance(plot=T,gam.k=10)

gns <- c("CD14",'FYB1',"CD8A")
x <- p2l$misc$varinfo
x$rank <- rank(-1*x$res)
x$rank <- rank(x$lp)
x <- x[order(x$lp,decreasing=F),]
x <- x[x$rank<5e3,]
x$name <- ''
x$signif <- x$lpa<log(0.05)
for(g in gns) x[g,'name'] <- g
pv <- ggplot(x,aes(x=rank,y=res,label=name,color=signif)) + geom_point_rast(size=0.1,raster.width=3,raster.height=2.8) +theme_bw() + xlab("gene overdispersion rank") + ylab("residual variance")+geom_text_repel(color=1) + geom_vline(xintercept = sum(x$lpa<log(0.05)),linetype=2,color=8) + guides(color='none') + scale_color_manual(values=adjustcolor(c('gray40','red'),alpha=0.15)) +ylim(0,4)

Ignoring unknown parameters: raster.width, raster.height

pv

pdf(file='resvar.pdf',width=3,height=2.8); print(pv); dev.off()

x <- p2l$misc$varinfo
x$rank <- rank(x$lp)
x <- x[order(x$lp,decreasing=F),]
x$name <- ''
x$signif <- x$lpa<log(0.05)
x <- na.omit(x)
for(g in gns) x[g,'name'] <- g
pmv <- ggplot(x,aes(x=m,y=v,label=name,color=signif)) + geom_point_rast(size=0.2,raster.width=3,raster.height=2.8) +theme_bw() + xlab("observed mean (log)") + ylab("observed variance (log)")+geom_text_repel(color=1,box.padding = 0.7) + guides(color='none') + scale_color_manual(values=c(adjustcolor('gray40',alpha=0.1),adjustcolor('red',alpha=0.15))) + geom_smooth(span=0.2,aes(x=m,y=v),color='blue',size=0.2)

Ignoring unknown parameters: raster.width, raster.height

pmv

pdf(file='pbmc10k_varnorm.pdf',width=3,height=2.8); print(pmv); dev.off()

Theoretical mean-variance

Downsample cells to a fixed size

# subsamples a cell 
subsample.cell <- function(m,depth) {
  if(sum(m)<depth) stop('not enough molecules')
  as.numeric(rmultinom(1,depth,m/sum(m)))
}

fixed.depth <- 5e3;

cdf <- apply(cd[,colSums(cd)>fixed.depth],2,subsample.cell,fixed.depth)
rownames(cdf) <- rownames(cd)
cdf <- as(cdf,'dgCMatrix')

vi <- intersect(colnames(cdf),names(ann)[ann=="CD4 Naive"])
cdf.mean <- rowSums(cdf[,vi])/length(vi)
cdf.var <- apply(cdf[,vi],1,var)

df <- data.frame(mean=cdf.mean,var=cdf.var)
p <- ggplot(df[df$mean<2,],aes(x=mean,y=var))+geom_point_rast(size=0.3,color=adjustcolor("gray40",alpha=0.3),raster.width=2,raster.height=2) + stat_smooth(method = "lm", formula = y ~ offset(x) + I(x^2), size = 1,linetype=2,col=2) +theme_bw() + geom_abline(slope=1,intercept = 0,col='darkgreen',linetype=2,size=1) + xlab('mean') + ylab("variance")
p

pdf(file='pbmc10k_stdvar.pdf',width=2,height=2); print(p); dev.off()

Log-transformed distribution

d <- de[['CD14+ Monocytes']]
g <- rownames(d)[115]; # IGSF6; LMO2; PTPRE; GLIPR1
x <- p2l$counts[names(ann)[ann=='CD14+ Monocytes'],g]*1e3
df <- data.frame(expr=x)
ggplot(df,aes(x=expr)) + geom_density(fill='#FF00003F') +theme_bw()

ggplot(df,aes(x=expr)) + geom_histogram(aes(y=..density..), colour="black", fill="wheat") +theme_bw()

ggplot(df,aes(x=log10(expr+1))) + geom_density(fill='#FF00003F') +theme_bw()

ggplot(df,aes(x=log10(expr+1))) + geom_histogram(aes(y=..density..), colour="black", fill="wheat") +theme_bw() 

a1 <- ggplot(df,aes(x=expr)) + geom_histogram( colour="black", fill="wheat") +theme_bw()  + ylab('cell frequency') + xlab('cpm (linear)') 
a2 <- ggplot(df,aes(x=log10(expr+1))) + geom_histogram(colour="black", fill="wheat") +theme_bw() + xlab('log10[ cpm+1 ]') + ylab('cell frequency') 
pp <- plot_grid(plotlist=list(a1,a2),nrow=2)
pp

pdf(file='pbmc10k_dist.pdf',width=3,height=3); print(pp); dev.off()

Poisson posterior

gn <- "GZMA"
gn <- "NPM1"
gn <- "IL32"
#gn <- "GNLY"
rng <- range(p2l$counts[,gn])/1e3
lambda.grid <- seq(0,rng[2]*1.2,length.out=200)

# sample cells
n.cells <- 5
set.seed(0)
cells1 <- sample(names(ann)[ann=='CD8 effector'],n.cells)
cells2 <- sample(names(ann)[ann=='CD8 Naive'],n.cells)

postm <- do.call(cbind,lapply(c(cells1,cells2),function(cell) {
   x<-dpois(cd[gn,cell],sum(cd[,cell])*lambda.grid)
   x/sum(x)
}))
colnames(postm) <- c(cells1,cells2);
rownames(postm) <- lambda.grid
df <- reshape2::melt(postm)
colnames(df) <- c('lambda','cell','p')
df$type <- ann[as.character(df$cell)]
#df$type <- as.factor(gsub("CD8 ","",as.character(df$type) ))

postj <- do.call(cbind,tapply(colnames(postm),ann[colnames(postm)],function(cn) {
  x <- exp(rowSums(log(postm[,cn])))
  x/sum(x)
}))
jdf <- reshape2::melt(postj);
colnames(jdf) <- c('lambda','type','p')

p <- ggplot(df,aes(x=lambda,y=p,group=cell,color=type)) + geom_line(alpha=0.3,size=0.3)+theme_bw() +
  geom_line(data=jdf,aes(x=lambda,y=p,group=type,color=type),size=1,alpha=0.8,linetype='43')+ xlim(0,0.0036) + 
  scale_color_manual(values=adjustcolor(rev(c('dodgerblue','indianred')),alpha=1)) + guides(color='none') + 
  theme(legend.position = c(0.75, 0.8),legend.background = element_rect(fill=alpha('white', 0.1))) + guides(color=guide_legend(title=element_blank())) +
  ylab('probability density') + xlab('expression magnitude') 
p

pdf(file='posteriors.pdf',width=3,height=2); print(p); dev.off()

null device 
          1 

Tree

a1 <- p2$plotEmbedding(type='PCA',embeddingType='tSNE',groups=ann,alpha=0.4, font.size=c(3,4))
#pdf(file='pmbc10k_tsne.pdf',width=3,height=3); print(a1); dev.off()
a1

source("../retina/functions_copy.R")

in 2 dimensions

emb <- p2$embeddings$PCA$tSNE
z2 <- t.ppt.tree(X=t(emb),emb=emb,M=200,lambda=100,sigma=0.02,metrics='euclidean',plot=F,output=F)

the condition has length > 1 and only the first element will be used

a1 <- p2$plotEmbedding(type='PCA',embeddingType='tSNE',groups=ann,alpha=0.4, size=0.1, font.size=c(3,4),mark.groups=T,raster=T, raster.dpi=600)+
  t.ggplot.ppt.tree(z2,emb,size=1,col=adjustcolor(1,alpha=0.3))+
  theme(panel.border = element_rect(color = 1, size=0.2,linetype=1),axis.line=element_blank())+  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
a1

pdf(file='pmbc10k_tree.pdf',width=3,height=3); print(a1); dev.off()

null device 
          1 

Dimensionality reduction and neural net figure

a1 <- p2$plotEmbedding(type='PCA',embeddingType='tSNE',groups=ann,alpha=0.4, raster=T, size=0.1, raster.dpi=600, font.size=c(3,4),title='t-SNE') +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
pdf(file='pmbc10k_tsne.pdf',width=3,height=3); print(a1); dev.off()

null device 
          1 

a1

ICA

library(fastICA)

odgenes <- rownames(p2$misc[["varinfo"]])[(order(p2$misc[["varinfo"]]$lp, decreasing = FALSE)[1:3e3])]
x <- p2$counts[,odgenes]
x@x <- x@x * rep(p2$misc[["varinfo"]][colnames(x), "gsf"], diff(x@p))

ica <- fastICA(x,2)

a5 <- sccore::embeddingPlot(ica$S,groups=ann,alpha=0.4, font.size=c(3,4),plot.theme=theme_bw(),title='ICA')
a5

NMF

library(NMF)
nmf <- nmf(as.matrix(x),rank=2)

nmf.emb <- nmf@fit@W;
a4 <- sccore::embeddingPlot(nmf.emb,groups=ann,alpha=0.4, font.size=c(3,4),plot.theme=theme_bw(),title='NMF')
a4

PCA

a2 <- sccore::embeddingPlot(p2$reductions$PCA[,c(1,2)],groups=ann,alpha=0.4, font.size=c(3,4),plot.theme=theme_bw(),title='PCA')
a2

load autoencoder results

aer <- readRDS("auto.emb.rds")

Combined figure

pp <- plot_grid(plotlist=lapply(list(
  p2$plotEmbedding(type='PCA',embeddingType='tSNE',groups=ann,alpha=0.4, font.size=c(3,4),title='t-SNE',raster=T),
  sccore::embeddingPlot(p2$reductions$PCA[,c(1,2)],groups=ann,mark.groups=F,alpha=0.4, font.size=c(4,5),plot.theme=theme_bw(),title='PCA',raster=T),
  sccore::embeddingPlot(nmf.emb,groups=ann,alpha=0.4, font.size=c(3,4),mark.groups=F,plot.theme=theme_bw(),title='NMF',raster=T),
  sccore::embeddingPlot(aer[[1]],groups=ann, mark.groups=F, alpha=0.4, font.size=c(3,4),plot.theme=theme_bw(),title='autoencoder',raster=T)
  ),function(g) g+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())),nrow=1)
pdf(file='reductions.pdf',width=12,height=3); print(pp); dev.off();

null device 
          1 

pp

a3 <- sccore::embeddingPlot(aer[[1]],groups=ann,alpha=0.4, font.size=c(3,4),plot.theme=theme_bw(),title='autoencoder')
a3

a3 <- sccore::embeddingPlot(aer[[1]],groups=ann,alpha=0.4, font.size=c(3,4),plot.theme=theme_bw(),title='embedding learning')
a3

a1 <- p2$plotEmbedding(type='PCA',embeddingType='tSNE',groups=ann,alpha=0.4, font.size=c(3,4),title='t-SNE')
#pdf(file='pmbc10k_tsne.pdf',width=3,height=3); print(a1); dev.off()
a1

vc <- setdiff(rownames(emb),rownames(x_test))
a1 <- sccore::embeddingPlot(aer[[3]][vc,],groups=ann,alpha=0.4,size=0.2, font.size=c(3,4),plot.theme=theme_bw(),title='predicted',raster=T)
a2 <- sccore::embeddingPlot(aer[[4]][vc,],groups=ann,alpha=0.4,size=0.2, font.size=c(3,4),plot.theme=theme_bw(),title='actual',raster=T)
pp <- plot_grid(plotlist=list(a2,a1),nrow=1)
pdf(file='ae.tSNE.pdf',width=8,height=4); print(pp); dev.off();

null device 
          1 

pp

Naive/memory T cells only

odgenes <- rownames(p2$misc[["varinfo"]])[(order(p2$misc[["varinfo"]]$lp, decreasing = FALSE)[1:3e3])]
x <- p2$counts[,odgenes]
x@x <- x@x * rep(p2$misc[["varinfo"]][colnames(x), "gsf"], diff(x@p))
x <- x[ann[rownames(x)] %in% c("CD4 Naive","CD4 Memory","CD8 Naive","CD8 effector","Double negative T cell"), ]
p2t <- Pagoda2$new(t(x),log.scale=T, n.cores=30)

3955 cells, 3000 genes; normalizing ... 
using plain model 
log scale ... 
done.

p2t$adjustVariance(plot=F,gam.k=10)

calculating variance fit ...
 using gam 
211 overdispersed genes ... 211
persisting ... 
done.

p2t$calculatePcaReduction(nPcs=10,n.odgenes=3e3)

running PCA using 3000 OD genes .
.
.
.
 done

a1 <- sccore::embeddingPlot(p2t$reductions$PCA[,c(1,2)],groups=ann,mark.groups=T,alpha=0.4, font.size=c(4,5),plot.theme=theme_bw(),title='T cell PCA',raster=T)
pdf(file='tcell.pca.pdf',width=3,height=3); print(a1); dev.off();

null device 
          1 

a1