R Notebook

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(ggplot2)
library(cowplot)

Attaching package: ‘cowplot’

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

    ggsave

Load data

mat <- readRDS("E12.5.rds")

source("~/m/pavan/DLI/conp2.r")

data.table 1.11.8  Latest news: r-datatable.com

doublet.f <- get.scrublet.scores(mat)

Basic processing

r <- basicP2proc(mat[,doublet.f[colnames(mat)]<0.2],n.cores=20,nPcs=30,make.geneknn=F)

11590 cells, 27998 genes; normalizing ... using plain model winsorizing ... log scale ... done.
calculating variance fit ... using gam 1373 overdispersed genes ... 1373persisting ... done.
running PCA using 3000 OD genes .... done
creating space of type angular done
adding data ... done
building index ... done
querying ... done
Estimating embeddings.

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

running tSNE using 20 cores:
 - point 10000 of 11590

r$makeKnnGraph(k = 50, type = "PCA", center = TRUE, weight.type = "none", n.cores = 30, distance = "cosine")

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

r$getKnnClusters(type='PCA',method=function(x) conos::leiden.community(x,r=0.3),name='leiden')

set.seed(3) # 4/1e7  3/1e7   4/3e7/M2/a0.8   5/M2/a0.5/3e7   5/M2/a0.5/3e7  5/M3/a0.2/5e7
M <- 2; r$getEmbedding(type='PCA', embeddingType='largeVis', M=M,  gamma = 1/M,sgd_batches=5e7,perplexity = NA)

Estimating embeddings.

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

set.seed(2)
r$getEmbedding(type='PCA',embeddingType='tSNE',perplexity = 400)

calculating distance ... pearson ...running tSNE using 20 cores:
 - point 10000 of 11590

r$plotEmbedding(type = 'PCA', embedding = 'tSNE', colors=r$counts[,'Mki67'],alpha=0.5,cex=0.2)

treating colors as a gradient with zlim: 0 0.2466913 

r$plotEmbedding(type = 'PCA', embedding = 'tSNE', colors=r$counts[,'Id2'],alpha=0.5,cex=0.2)

treating colors as a gradient with zlim: 0 0.6343238 

r$plotEmbedding(type = 'PCA', embedding = 'tSNE', clusterType = 'leiden', mark.clusters = T,alpha=0.2,cex=0.1,mark.cluster.cex = 1)

par(mfrow=c(3,3))
gns <- c('Klk13','Pitx2','Sparcl1','Col1a1','Myl9','Stmn2','Acta2','Cav1','Wt1')

lapply(gns,function(gene) r$plotEmbedding(type = 'PCA', embedding = 'tSNE', colors=r$counts[,gene],alpha=0.5,cex=0.2,main=gene))

treating colors as a gradient with zlim: 0 0.3347215 
treating colors as a gradient with zlim: 0 0.3414924 
treating colors as a gradient with zlim: 0 0.454883 
treating colors as a gradient with zlim: 0 0.3963445 
treating colors as a gradient with zlim: 0 0.2912576 
treating colors as a gradient with zlim: 0 0.03126894 
treating colors as a gradient with zlim: 0 1.99891 
treating colors as a gradient with zlim: 0 0.4640115 
treating colors as a gradient with zlim: 0 0.341166 
[[1]]
NULL

[[2]]
NULL

[[3]]
NULL

[[4]]
NULL

[[5]]
NULL

[[6]]
NULL

[[7]]
NULL

[[8]]
NULL

[[9]]
NULL

par(mfrow=c(3,3))
gns <- c('Cpa1','Spp1','Neurog3','Col1a1','Fabp7','Tlx2','Pecam1','Col3a1','Hbb-bt')

lapply(gns,function(gene) r$plotEmbedding(type = 'PCA', embedding = 'tSNE', colors=r$counts[,gene],alpha=0.5,cex=0.2,main=gene))

treating colors as a gradient with zlim: 0 2.03234 
treating colors as a gradient with zlim: 0 1.351252 
treating colors as a gradient with zlim: 0 0.3198094 
treating colors as a gradient with zlim: 0 0.3963445 
treating colors as a gradient with zlim: 0 0.08929173 
treating colors as a gradient with zlim: 0 0.1353626 
treating colors as a gradient with zlim: 0 0.5760832 
treating colors as a gradient with zlim: 0 0.4638031 
treating colors as a gradient with zlim: 0 1.341928 
[[1]]
NULL

[[2]]
NULL

[[3]]
NULL

[[4]]
NULL

[[5]]
NULL

[[6]]
NULL

[[7]]
NULL

[[8]]
NULL

[[9]]
NULL

par(mfrow=c(3,3))
gns <- c('Krt19','Stmn2','Acta2','Col1a1','Fabp7','Rac2','Cxcl12','Nkx2-5','Sparcl1')

lapply(gns,function(gene) r$plotEmbedding(type = 'PCA', embedding = 'largeVis', colors=r$counts[,gene],alpha=0.5,cex=0.2,main=gene))

treating colors as a gradient with zlim: 0 0.3041981 
treating colors as a gradient with zlim: 0 0.03126894 
treating colors as a gradient with zlim: 0 1.99891 
treating colors as a gradient with zlim: 0 0.3963445 
treating colors as a gradient with zlim: 0 0.08929173 
treating colors as a gradient with zlim: 0 0.2900099 
treating colors as a gradient with zlim: 0 0.2836841 
treating colors as a gradient with zlim: 0 0.3769584 
treating colors as a gradient with zlim: 0 0.454883 
[[1]]
NULL

[[2]]
NULL

[[3]]
NULL

[[4]]
NULL

[[5]]
NULL

[[6]]
NULL

[[7]]
NULL

[[8]]
NULL

[[9]]
NULL

r$plotEmbedding(type = 'PCA', embedding = 'largeVis', colors=r$counts[,'Cd44'],alpha=0.5,cex=0.2)

treating colors as a gradient with zlim: 0 0.2337546 

a2 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,groups=r$clusters$PCA$leiden,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())
a2

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

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

de <- r$diffgenes$PCA[['leiden']];
conos::plotDEheatmap(r,r$clusters$PCA$leiden,de=de,n.genes.per.cluster = 10)

#r$plotGeneHeatmap(genes=rownames(de)[1:30],groups=r$clusters$PCA[[1]])

Can’t tell the identity of the two Hbb-bt populations … both look like RBCs 1/17 - erythroblast/erythroid … 17 looks more primitive 15 seems to be myeloid / macrophages C1qc-high 9,5: neuroadrenergic … actively dividing though? 13 - vascular endothelial

Should align with TabMuris data

ann <- r$clusters$PCA$leiden
ann <- setNames(as.character(ann),names(ann))
ann[ann == '8'] <- 'vascular endothelial'
ann[ann %in% c('7')] <- 'neuroadrenergic'
ann[ann %in% c('5')] <- 'erythroid'
ann[ann %in% c('2')] <- 'mesothelium'
ann[ann %in% c('3')] <- 'VSM'
ann[ann %in% c('1','4')] <- 'epithelial'
ann[ann %in% c('6')] <- 'myeloid'
ann <- as.factor(ann)

a1 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,groups=as.factor(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())
a2 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,groups=r$clusters$PCA$leiden,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())
a3 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,colors=pagoda2:::val2col(r$counts[,"Mki67"]),size=0.1,alpha=0.5,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())
plot_grid(plotlist=list(a1,a2,a3),nrow=1)

set.seed(5)
r$getKnnClusters(type='PCA',method=function(x) conos::leiden.community(x,r=5),name='leiden2')

a1 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,groups=as.factor(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())
a2 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,groups=r$clusters$PCA$leiden2,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())
a3 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,colors=pagoda2:::val2col(r$counts[,"Mki67"]),size=0.1,alpha=0.5,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())
plot_grid(plotlist=list(a1,a2,a3),nrow=1)

Make a coarse clusering for illustrating cluster similarity

mcl <- setNames(rep(NA,length(ann)),names(ann))
mcl[ann == 'epithelial'] <- 'E1'
mcl[r$clusters$PCA$leiden2[names(mcl)] %in% c(16,22,17,20)] <- 'E2'
mcl[ann == 'mesothelium'] <- 'M1'
mcl[r$clusters$PCA$leiden2[names(mcl)] %in% c(23,26,31,35)] <- 'M2'
mcl[ann == 'VSM'] <- 'V1'
mcl[r$clusters$PCA$leiden2[names(mcl)] %in% c(32,27)] <- 'V2'
mcl <- as.factor(mcl)

set.seed(0)
a1 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,groups=as.factor(ann),size=0.1,alpha=0.2,raster=T,raster.height = 3,raster.width=3,font.size=c(4,5))+
  theme(panel.border = element_rect(color = 1, size=0.2,linetype=1),axis.line=element_blank())
a2 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,groups=mcl,size=0.1,alpha=0.2,raster=T,raster.height = 3,raster.width=3,font.size=c(5.5),palette = function(n) rainbow(n,s=0.7,v=0.7))+
  theme(panel.border = element_rect(color = 1, size=0.2,linetype=1),axis.line=element_blank())
a3 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,colors=pagoda2:::val2col(r$counts[,"Mki67"]),size=0.1,alpha=0.5,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())
plot_grid(plotlist=list(a1,a2,a3),nrow=1)

Hierarchical clustering

cp <- conos:::collapseCellsByType(r$misc$rawCounts,mcl)
cp <- log10(cp/rowSums(cp)*1e6+1)

hc <- hclust(as.dist(1-cor(t(cp))),method='ward.D2')
plot(hc,main='',sub='',xlab='',yaxt='none',ylab='');

pdf(file='dend.pdf',width=3,height=3); plot(hc,main='',sub='',xlab='',yaxt='none',ylab=''); dev.off();

Clean plots:

set.seed(2)
a1 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,groups=as.factor(ann),size=0.1,alpha=0.2,raster=T,raster.height = 3,raster.width=3,font.size=c(4,5),palette = function(n) sample(rainbow(n,s=0.9,v=0.9)),box.padding=0.5)+
  theme(panel.border = element_rect(color = 1, size=0.2,linetype=1),axis.line=element_blank(),plot.margin = margin(.1, .1, .1, .1, "cm"))
pdf(file='fp_celltypes.pdf',width=3,height=3); print(a1); dev.off()

null device 
          1 

a1

Clusters and markers:

plst <- theme(panel.border = element_rect(color = 1, size=0.2,linetype=1),axis.line=element_blank(),plot.margin = margin(.1, .1, .1, .1, "cm"));
set.seed(0)
a2 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,groups=mcl,size=0.1,alpha=0.2,raster=T,raster.height = 3,raster.width=3,font.size=c(5.5),palette = function(n) rainbow(n,s=0.7,v=0.7))+plst

a3 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,colors=pagoda2:::val2col(r$counts[,"Epcam"]),size=0.1,alpha=0.5,raster=T,raster.height = 3,raster.width=3,font.size=c(3,4))+plst+ annotate('text',x=-Inf,y=Inf,hjust=0,vjust=1.2,label='Epcam',size=6)
a4 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,colors=pagoda2:::val2col(r$counts[,"Col3a1"]),size=0.1,alpha=0.5,raster=T,raster.height = 3,raster.width=3,font.size=c(3,4))+plst+ annotate('text',x=-Inf,y=Inf,hjust=0,vjust=1.2,label='Col3a1',size=6)
a5 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,colors=pagoda2:::val2col(r$counts[,"Mki67"]),size=0.1,alpha=0.5,raster=T,raster.height = 3,raster.width=3,font.size=c(3,4))+plst+ annotate('text',x=-Inf,y=Inf,hjust=0,vjust=1.2,label='Mki67',size=6)
pp <- plot_grid(plotlist=list(a2,a3,a4,a5),nrow=2)
pdf(file='fp_panels.pdf',width=4,height=4); print(pp); dev.off();

png 
  2 

pp

mcl2 <- setNames(as.character(mcl),names(mcl))
mcl2[is.na(mcl2)] <- ' '
mcl2 <- as.factor(mcl2)

plst <- theme(panel.border = element_rect(color = 1, size=0.2,linetype=1),axis.line=element_blank(),plot.margin = margin(.1, .1, .1, .1, "cm"));
set.seed(0)
a2 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,groups=mcl2,size=0.2,alpha=0.2,raster=T,mark.groups=F,raster.height = 3,raster.width=3,font.size=c(5.5),palette = function(n) c('gray60',rainbow(n-1,s=0.7,v=0.7)))+plst

a3 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,colors=pagoda2:::val2col(r$counts[,"Epcam"]),size=0.1,alpha=0.5,raster=T,raster.height = 3,raster.width=3,font.size=c(3,4))+plst+ annotate('text',x=-Inf,y=Inf,hjust=0,vjust=1.2,label='Epcam',size=6)
a4 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,colors=pagoda2:::val2col(r$counts[,"Col3a1"]),size=0.1,alpha=0.5,raster=T,raster.height = 3,raster.width=3,font.size=c(3,4))+plst+ annotate('text',x=-Inf,y=Inf,hjust=0,vjust=1.2,label='Col3a1',size=6)
a5 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,colors=pagoda2:::val2col(r$counts[,"Mki67"]),size=0.1,alpha=0.5,raster=T,raster.height = 3,raster.width=3,font.size=c(3,4))+plst+ annotate('text',x=-Inf,y=Inf,hjust=0,vjust=1.2,label='Mki67',size=6)
pp <- plot_grid(plotlist=list(a2,a3,a4,a5),nrow=1)
pdf(file='fp_panels2.pdf',width=8,height=2); print(pp); dev.off();

png 
  2 

pp

Combined figure

plst <- theme(panel.border = element_rect(color = 1, size=0.2,linetype=1),axis.line=element_blank(),plot.margin = margin(.1, .1, .1, .1, "cm"));
ann.size <- 7;
emb <- r$embeddings$PCA$tSNE
a1 <- conos::embeddingPlot(emb,groups=as.factor(ann),size=0.1,alpha=0.2,raster=T,raster.height = 3,raster.width=3,font.size=c(3.5,5.5))+plst
set.seed(0)
a2 <- conos::embeddingPlot(emb,groups=mcl,size=0.1,alpha=0.2,raster=T,raster.height = 3,raster.width=3,font.size=c(5.5),palette = function(n) rainbow(n,s=0.7,v=0.7),plot.na=F) +xlim(range(emb[,1]))+ylim(range(emb[,2]))+plst

a3 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,colors=pagoda2:::val2col(r$counts[,"Epcam"]),size=0.1,alpha=0.5,raster=T,raster.height = 3,raster.width=3,font.size=c(3,4))+plst+ annotate('text',x=-Inf,y=Inf,hjust=-0.1,vjust=1.2,label='Epcam',size=ann.size)
a4 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,colors=pagoda2:::val2col(r$counts[,"Col3a1"]),size=0.1,alpha=0.5,raster=T,raster.height = 3,raster.width=3,font.size=c(3,4))+plst+ annotate('text',x=-Inf,y=Inf,hjust=-0.1,vjust=1.2,label='Col3a1',size=ann.size)
a5 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,colors=pagoda2:::val2col(r$counts[,"Mki67"]),size=0.1,alpha=0.5,raster=T,raster.height = 3,raster.width=3,font.size=c(3,4))+plst+ annotate('text',x=-Inf,y=Inf,hjust=-0.1,vjust=1.2,label='Mki67',size=ann.size)
pp <- plot_grid(plotlist=list(a1,a2,a3,a4,a5),nrow=1)
pdf(file='fp_panels_wide.pdf',width=15,height=3); print(pp); dev.off();

png 
  2 

pp

r$plotEmbedding(type = 'PCA', embedding = 'largeVis', groups=ann, mark.clusters = T,alpha=0.2,cex=0.2, mark.cluster.cex = 1)

a1 <- conos::embeddingPlot(r$embeddings$PCA$largeVis,groups=r$clusters$PCA$multilevel,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())
a1

a1 <- conos::embeddingPlot(r$embeddings$PCA$tSNE,groups=as.factor(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())
a1