PBMC examples

library(pagoda2)

Loading required package: Matrix
Loading required package: igraph
package ‘igraph’ was built under R version 3.5.3
Attaching package: ‘igraph’

The following objects are masked from ‘package:stats’:

    decompose, spectrum

The following object is masked from ‘package:base’:

    union

library(Matrix)
library(ggplot2)
library(cowplot)

Attaching package: ‘cowplot’

The following object is masked from ‘package:ggplot2’:

    ggsave

Read in data

library(hdf5r)
h5_data <- H5File$new("5k_pbmc_NGSC3_aggr_filtered_feature_bc_matrix.h5")
cd <- Matrix::sparseMatrix(
  i = h5_data[['matrix/indices']][],
  p = h5_data[['matrix/indptr']][],
  x = h5_data[['matrix/data']][],
  dimnames = list(
    h5_data[['matrix/features/name']][],
    h5_data[['matrix/barcodes']][]
  ),
  dims = h5_data[['matrix/shape']][],
  index1 = FALSE
)
h5_data$close_all()

hist(log10(rowSums(counts)+1),main='Molecules per gene',xlab='molecules (log10)',col='cornsilk')
abline(v=1,lty=2,col=2)

counts <- counts[rowSums(counts)>=100,]
dim(counts)

[1] 15612 65833

Doublets

ds <- get.scrublet.scores(counts)

Written 5.2% of 65833 rows in 2 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 37 secs.      
Written 8.6% of 65833 rows in 3 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 32 secs.      
Written 12.0% of 65833 rows in 4 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 30 secs.      
Written 15.4% of 65833 rows in 5 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 28 secs.      
Written 18.8% of 65833 rows in 6 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 26 secs.      
Written 21.5% of 65833 rows in 7 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 26 secs.      
Written 24.7% of 65833 rows in 8 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 25 secs.      
Written 27.9% of 65833 rows in 9 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 24 secs.      
Written 31.1% of 65833 rows in 10 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 22 secs.      
Written 34.5% of 65833 rows in 11 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 21 secs.      
Written 37.9% of 65833 rows in 12 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 20 secs.      
Written 41.1% of 65833 rows in 13 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 19 secs.      
Written 44.5% of 65833 rows in 14 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 18 secs.      
Written 47.9% of 65833 rows in 15 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 16 secs.      
Written 51.5% of 65833 rows in 16 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 15 secs.      
Written 54.9% of 65833 rows in 17 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 14 secs.      
Written 58.3% of 65833 rows in 18 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 13 secs.      
Written 61.8% of 65833 rows in 19 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 12 secs.      
Written 65.2% of 65833 rows in 20 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 11 secs.      
Written 68.6% of 65833 rows in 21 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 9 secs.      
Written 72.0% of 65833 rows in 22 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 8 secs.      
Written 75.4% of 65833 rows in 23 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 7 secs.      
Written 78.8% of 65833 rows in 24 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 6 secs.      
Written 82.2% of 65833 rows in 25 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 5 secs.      
Written 85.6% of 65833 rows in 26 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 4 secs.      
Written 89.0% of 65833 rows in 27 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 3 secs.      
Written 92.4% of 65833 rows in 28 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 2 secs.      
Written 96.0% of 65833 rows in 29 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 1 secs.      
Written 99.7% of 65833 rows in 30 secs using 1 thread. anyBufferGrown=no; maxBuffUsed=49%. ETA 0 secs.      
                                                                                                                                     

Annotation

ann <- read.delim('68k_pbmc_barcodes_annotation.tsv',header=T,sep='\t')

cannot open file '68k_pbmc_barcodes_annotation.tsv': No such file or directoryError in file(file, "rt") : cannot open the connection

Conos processing together with an annotated PBMC dataset

library(pagoda2)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:igraph’:

    as_data_frame, groups, union

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(conos)
library(parallel)
library(cowplot)
library(ggrepel)
library(Matrix)

Co-align with Rahul’s 10k PBMCs

load("../pbmc10k_v3/seurat.RData") # cd, cd.var,meta, umap

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

Pre-process, align and annotate

cdl.p2 <- lapply(list(pbmc60=counts,pbmc10=cd), basicP2proc, n.cores=30, min.cells.per.gene=0, n.odgenes=2e3, get.largevis=FALSE, make.geneknn=FALSE,get.tsne=T)

con <- Conos$new(cdl.p2,n.cores=30)
con$buildGraph(k=15, k.self=5, space='PCA', ncomps=30, n.odgenes=2000, matching.method='mNN', metric='angular', verbose=TRUE)

found 0 out of 1 cached PCA  space pairs ... running 1 additional PCA  space pairs . done
inter-sample links using  mNN  . done
local pairs local pairs  done
building graph ..done

Cluster and generate a 2D ebmedding

con$findCommunities(r=2)
con$embedGraph(method='UMAP'); 

Convert graph to adjacency list...
Done
Estimate nearest neighbors and commute times...
Estimating hitting distances: 17:26:54.

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Done.
Estimating commute distances: 17:28:32.
Hashing adjacency list: 17:28:32.

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**********************

Done.

****************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|

Estimating distances: 17:28:43.

**************************************************|

Done
Done.
All done!: 17:28:56.
Done
Estimate UMAP embedding...

Maximal number of estimated neighbors is 31. Consider increasing min.visited.verts, min.prob or min.prob.lower.
17:28:56 UMAP embedding parameters a = 0.0267 b = 0.7906
17:28:56 Read 75265 rows and found 1 numeric columns
17:28:57 Commencing smooth kNN distance calibration using 30 threads
17:29:00 Initializing from normalized Laplacian + noise
17:29:22 Commencing optimization for 1000 epochs, with 2092332 positive edges using 30 threads
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
17:29:48 Optimization finished

Done

emb <- con$embedding;

con$misc$embeddings <- list(UMAP=emb)

con$embedGraph(method='largeVis'); 

Estimating embeddings.

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

con$misc$embeddings <- list(largeVis=emb)

propagate labels

new.ann <- con$propagateLabels(labels = satija.ann, verbose=T,fixed.initial.labels=TRUE)

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
*********

Stop after 15 iterations. Norm: 0.0237014
Min weight: 1.67017e-05, max weight: 0.367879, fading: (10, 0.1)

*****************************************|

new.ann <- new.ann$labels

Take a look:

alpha <- 0.1; size <- 0.1;
p1 <- con$plotGraph(title='clusters',alpha=alpha,size=size);
p2 <- con$plotGraph(color.by='sample',title='samples',mark.groups=F,alpha=alpha,size=size);
p3 <- con$plotGraph(groups=new.ann,alpha=alpha,size=size,title='annotation',plot.na=F)
p4 <- con$plotGraph(groups=satija.ann,alpha=alpha,size=size,title='old annotation',plot.na=F)
cowplot::plot_grid(plotlist=list(p1,p2,p3,p4),nrow=2)

alpha <- 0.1; size <- 0.1;
p2 <- con$plotGraph(color.by='sample',title='samples',mark.groups=F,alpha=alpha,size=size);
p3 <- con$plotGraph(groups=new.ann,alpha=alpha,size=size,title='annotation',plot.na=F)
p4 <- con$plotGraph(groups=satija.ann,alpha=alpha,size=size,title='old annotation',plot.na=F)
p1 <- con$plotGraph(gene='CD14')
cowplot::plot_grid(plotlist=list(p1,p2,p3,p4),nrow=2)

p2 processing

r <- Pagoda2$new(counts[,names(ds)[ds<0.2]],log.scale=T, n.cores=30)

59525 cells, 15612 genes; normalizing ... using plain model log scale ... done.

#p2 <- con$samples$pbmc60

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

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

r$calculatePcaReduction(nPcs=50,n.odgenes=3e3)

running PCA using 3000 OD genes .... done

r$makeKnnGraph(k=40,type='PCA',center=T,distance='cosine');

creating space of type angular done
adding data ... done
building index ... done
querying ... done

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

Too many cells to pre-calcualte correlation distances, switching to L2. Please, consider using UMAP.

running tSNE using 30 cores:
 - point 10000 of 59525
 - point 20000 of 59525
 - point 30000 of 59525
 - point 40000 of 59525
 - point 50000 of 59525

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

r$plotEmbedding(type='PCA',embeddingType='tSNE',clusterType='leiden',show.legend=F,mark.clusters=T,min.group.size=1,shuffle.colors=F,mark.cluster.cex=1,alpha=0.1,cex=0.01,main='clusters (tSNE)')

r$plotEmbedding(type='PCA',embeddingType='tSNE',groups=new.ann,show.legend=F,mark.clusters=T,min.group.size=1,shuffle.colors=F,mark.cluster.cex=1,alpha=0.01)

using provided groups as a factor

r$plotEmbedding(type='PCA',embeddingType='tSNE',colors=r$counts[,'FOXP3'],show.legend=F,mark.clusters=T,min.group.size=1,shuffle.colors=F,mark.cluster.cex=1,alpha=0.01)

treating colors as a gradient with zlim: 0 1.04557 

de <- r$getDifferentialGenes(type='PCA',verbose=T,clusterType='leiden',append.auc=T,upregulated.only = T)

running differential expression with  20  clusters ... adjusting p-values ... done.

conos:::plotDEheatmap(r,r$clusters$PCA$leiden,de)

Show evidence for cluster 1 and 2

p1 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,groups=as.factor(r$clusters$PCA$leiden),size=0.1,alpha=0.2,raster=T,raster.height = 3,raster.width=3)+
  theme(panel.border = element_rect(color = 1, size=0.2,linetype=1),axis.line=element_blank())
p1

p2 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,groups=new.ann,size=0.1,alpha=0.2,raster=T,raster.height = 3,raster.width=3,font.size=c(3,4))+
  theme(panel.border = element_rect(color = 1, size=0.2,linetype=1),axis.line=element_blank())
p2

gns <- c("CD14",'CD3E',"CD8A")
po <- lapply(gns,function(g) {
  conos::embeddingPlot(r$embeddings$PCA$tSNE,colors=r$counts[,g],size=0.1,alpha=0.6,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())
})
pp <- plot_grid(plotlist=c(list(p1),po),nrow=1)
pp

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

null device 
          1 

Differential expression ladder

cell.groups <- r$clusters$PCA$leiden
cell.groups <- tapply(names(cell.groups),cell.groups,I)
cell.groups <- cell.groups[c('1','2')]
names(cell.groups) <- c("CD4+ T cells","CD14+ Monocytes")

sig.thr <- qnorm(0.01,lower.tail=F);

cell.grid <- 10^(seq(1,log10(1e4),length.out=100))

ngl <- mclapply(cell.grid,function(n.cells) {
  cells <- lapply(cell.groups,sample,n.cells)
  cf <- as.factor(setNames(rep(names(cells),each=n.cells),unlist(cells)))
  suppressWarnings(x <- r$getDifferentialGenes(groups = cf,append.specificity.metrics=F)[[1]])
},mc.cores=30)

ng <- lapply(ngl,function(x) {sum(abs(x$Z)>=sig.thr)})

x <- r$counts[rownames(r$counts)%in% unlist(cell.groups),]
exp5 <- colnames(x)[colSums(x>0) >= nrow(x)*0.05]
exp10 <- colnames(x)[colSums(x>0) >= nrow(x)*0.1]

ng5 <- lapply(ngl,function(x) {sum(abs(x$Z[rownames(x) %in% exp5])>=sig.thr)})

df <- data.frame(n=cell.grid,nde=unlist(ng5),ndf=unlist(ng5)/length(exp5))
p <- ggplot(df,aes(x=log10(n),y=nde))+geom_point(size=0.2)+geom_smooth(color=adjustcolor('red',alpha=0.6),span = 0.5)+theme_bw() +
  xlab("log10[ number of cells ]")+ylab('number of DE genes') +
  scale_y_continuous(sec.axis = sec_axis(~./length(exp5), name = "DE gene fraction")) +
  theme(axis.text.y = element_text(angle = 90, hjust = 0.5),axis.text.y.right = element_text(angle = 90, hjust = 0.5))
p

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

null device 
          1