Detection difference example

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(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)

Loading required package: ggplot2

Attaching package: ‘cowplot’

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

    ggsave

library(ggrepel)
library(Matrix)
library(data.table)

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

Attaching package: ‘data.table’

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

    between, first, last

# some extra code for the comparative analysis here
source("/home/pkharchenko/m/pavan/DLI/conp2.r")

Load data

load("normal.RData")

cdl <- lapply(pagoda2:::sn(names(cdl)),function(n) { x <- cdl[[n]]; colnames(x) <- paste(n,colnames(x),sep='_'); x})
cdml <- lapply(pagoda2:::sn(names(cdl)),function(n) { x <- cdml[[n]]; rownames(x) <- paste(n,rownames(x),sep='_'); x})

cell.type.f <- as.factor(setNames(unlist(lapply(cdml,function(x) x$cell_line)),unlist(lapply(cdml,rownames))))

Integrate with Conos

pre-processing

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

Integrate

con <- Conos$new(matl.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)

Cluster and generate a 2D ebmedding

con$findCommunities(r=1) # 1.5 means slightly higher than default resolution 
con$embedGraph(method='UMAP'); 
con$misc$embeddings <- list(UMAP=con$embedding)

con$misc$embeddings$largeVis <- con$embedGraph()

Take a look:

alpha <- 0.2; size <- 0.8;
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=cell.type.f,alpha=alpha,size=size,title='annotation',font.size=c(8,8))
p4 <- con$plotGraph(gene='ENSG00000148773',alpha=0.8,size=size,title='Mki67')
cowplot::plot_grid(plotlist=list(p3,p2,p1,p4),nrow=2)

Look at gene scatter plots

mat <- con$getJointCountMatrix(raw=T)
mat <- mat/rowSums(mat)*1e6

Calculate average magnitudes per cell type per platform

cta <- tapply(rownames(mat),list(cell.type.f[rownames(mat)],sample.f[rownames(mat)]),function(x) colMeans(mat[x,]))
odf <- do.call(rbind,tapply(rownames(mat),list(cell.type.f[rownames(mat)],sample.f[rownames(mat)]),function(x) c(as.character(cell.type.f[x[1]]),as.character(sample.f[x[1]])),simplify = F))
names(cta) <- rownames(odf) <- apply(odf,1,paste,collapse='-')
sample.f <- con$getDatasetPerCell()

For every pair of platforms, draw three scatter plots comparing different techniques:

cell.lines <- unique(odf[,1])
tpairs <- combn(unique(odf[,2]),2)

i <- 2;
par(mfrow=c(1,3), mar = c(2.5,2.5,2.5,0.5), mgp = c(1.5,0.65,0), cex = 0.85)
lapply(cell.lines,function(n) {
  nn <- paste(n,tpairs[,i],sep='-')
  df <- do.call(cbind,cta[nn]);
  colnames(df) <- c('x','y')
  df <- log10(df+1);
  smoothScatter(df[,1],df[,2],xlab=nn[1],ylab=nn[2])
  abline(lm(y~x+0,data=data.frame(df)),col=2,lty=2)
})

Hihghlight outlier genes

i <- 2;
n <- cell.lines[[1]]; 
nn <- paste(n,tpairs[,i],sep='-')
df <- do.call(cbind,cta[nn]);
colnames(df) <- c('x','y')
df <- log10(df+1);
m <- lm(y~x+0,data=data.frame(df))
thr <- 1;
genes.up <- names(residuals(m))[residuals(m)>thr]
genes.down <- names(residuals(m))[residuals(m)< -1*thr]

pdf(file='platforms.pdf',width=9,height=3.2)
par(mfrow=c(1,3), mar = c(2.5,2.5,2.5,0.5), mgp = c(1.5,0.65,0), cex = 0.85)
lapply(cell.lines,function(n) {
  nn <- paste(n,tpairs[,i],sep='-')
  df <- do.call(cbind,cta[nn]);
  colnames(df) <- c('x','y')
  df <- df[rowSums(df)>0,]
  df <- log10(df+1);
  smoothScatter(df[,1],df[,2],xlab=tpairs[1,i],ylab=tpairs[2,i],main=paste(n,'cells'),colramp = function(n) adjustcolor(colorRampPalette(c("white", blues9))(n),alpha=0.7))
  abline(lm(y~x+0,data=data.frame(df)),col=1,lty=3)
  if(n=='H1975') {
    abline(1,coefficients(m)[1],col='red',lty=2)
    abline(-1,coefficients(m)[1],col='darkgreen',lty=2)
  }
  vg <- rownames(df) %in% genes.up
  points(df[vg,1],df[vg,2],col='red',pch=19,cex=0.2)
  vg <- rownames(df) %in% genes.down
  points(df[vg,1],df[vg,2],col='darkgreen',pch=19,cex=0.3)
})

[[1]]
NULL

[[2]]
NULL

[[3]]
NULL

dev.off()

null device 
          1