Number of cell figure

library(ggplot2)
library(mgcv)

Loading required package: nlme
This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.

require(purrr)

Loading required package: purrr

library(dplyr)

Attaching package: ‘dplyr’

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

    collapse

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

    filter, lag

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

    intersect, setdiff, setequal, union

library(ggpubr)

Loading required package: magrittr

Attaching package: ‘magrittr’

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

    set_names

library(reshape2)
library(Matrix)

Load simulation results

rl <- readRDS("pc2.rl.rds")

Make a simplified data frame

adf <- do.call(rbind,unlist(c(rl),recursive=F))
adf$auc[is.nan(adf$auc)] <- 0.5;
adf$auc[is.na(adf$auc)] <- 1;

adf %<>% group_by(n,c) %>% summarize(auc=mean(auc)) %>% as.data.frame  # add trim?

view raw data

ggplot(adf,aes(log10(n),log10(c),fill=auc))+geom_point(aes(col=auc)) + scale_color_gradient2(low = "blue", mid = "white", high = "red",midpoint=0.75) + theme_bw() + xlim(1,5)

2D regression fit

x <- adf; x$n <- log10(x$n); x$c <- log10(x$c)
sg <- gam(auc~s(n,c),data=x)

Predict on a regular grid

nbins <- 200;
pdf <- data.frame(n=seq(0,5,length.out=nbins),c=seq(min(x$c),max(x$c),length.out=nbins))
pdf <- data.frame(t(matrix(unlist(cross2(pdf$n,pdf$c)),nrow=2))); colnames(pdf) <- c('n','c')

pdf$auc <- predict(sg,newdata=pdf)
pdf$auc[pdf$auc>1] <- 1; pdf$auc[pdf$auc<0.5] <- 0.5

p <- ggplot(pdf,aes(n,c,z=auc))+geom_raster(aes(fill=auc),interpolate=T)  +
  scale_x_continuous(limits=c(1,5), expand = c(0, 0)) + 
  scale_y_continuous( expand = c(0, 0), limits=c(min(pdf$c),4.5)) +
  #scale_fill_gradientn(colours=adjustcolor(c("#0000FFFF","#0000FFFF","#0000FFFF","#FFFFFFFF","#FF0000FF","#FF0000FF"),alpha=0.65)) +
  scale_fill_gradientn(colours=adjustcolor(c("black","black",'black','white','darkgreen','darkgreen'),alpha=0.75),name='AUC') + 
  xlab('log10[ number of cells ]') +ylab('log10[ UMIs/cell ]') 
p

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

null device 
          1 

Clustering examples

Load sample data

cd <- readRDS(file=paste0('lcd_',1,'.rds'))
cell.sizes <- Matrix::colSums(cd)
gf <- attr(cd,'group')
cell.groups <- tapply(1:ncol(cd),gf,I)

quick.df <- function(n.cells,coverage) {
  # sample cells
  sampled.cells <- as.vector(t(matrix(unlist(lapply(cell.groups,sample,n.cells)),ncol=2)))
  sm <- cd[,sampled.cells,drop=F]
  sub.p <- min(1,coverage/mean(Matrix::colSums(sm)))
  # downsample coverage
  sm <- subsample.matrix(sm,p.sample=sub.p)
  sgf <- gf[colnames(sm)]
  # PCA/lda
  p2 <- Pagoda2$new(sm, n.cores = 1, trim = 0, log.scale =T, min.cells.per.gene = 1,  min.transcripts.per.cell = 0,verbose=F)
  p2$adjustVariance(plot = F, gam.k = 10,verbose=F)
  suppressWarnings(p2$calculatePcaReduction(nPcs = 2, n.odgenes = 1e3, maxit = 1000,verbose=F))
  lr <- lda(p2$reductions$PCA,sgf)
 
  # plot data frame
  df <- data.frame(p2$reductions$PCA); colnames(df) <- c('PC1','PC2'); df$type <- sgf
  df
}

quick.margin.plot <- function(df, alpha=0.5, size=1, trim.outliers=0) {
  # thanks to https://stackoverflow.com/questions/8545035/scatterplot-with-marginal-histograms-in-ggplot2
  require(cowplot)
  require(ggrastr)
  if(trim.outliers>0) {
    pc1t <- rank(abs(df$PC1)) <= nrow(df)*(1-trim.outliers);
    pc2t <- rank(abs(df$PC2)) <= nrow(df)*(1-trim.outliers);
    df <- df[pc1t & pc2t,]
  }
  pal <- scale_color_manual(values=c("blue","red"))
  # Set up scatterplot
  scatterplot <- ggplot(df, aes(x = PC1, y = PC2, color = type)) +
    geom_point_rast(size = size, alpha = alpha, width=1,height=1) +
    guides(color = FALSE) +
    theme_bw() +
    theme(plot.margin = margin()) + 
    scale_x_continuous(labels = NULL,breaks=NULL) +
    scale_y_continuous(labels = NULL,breaks=NULL) +
    pal
  
  
  # Define marginal histogram
  marginal_distribution <- function(x, var, group) {
    ggplot(x, aes_string(x = var, fill = group)) +
      #geom_histogram(bins = 30, alpha = 0.4, position = "identity") +
      geom_density(alpha = 0.4, size = 0.1) +
      guides(fill = FALSE) +
      theme_void() +
      theme(plot.margin = margin()) + scale_fill_manual(values=c("blue","red"))
  }
  
  # Set up marginal histograms
  x_hist <- marginal_distribution(df, "PC1", "type")
  y_hist <- marginal_distribution(df, "PC2", "type") +
    coord_flip()
  
  # Align histograms with scatterplot
  aligned_x_hist <- align_plots(x_hist, scatterplot, align = "v")[[1]]
  aligned_y_hist <- align_plots(y_hist, scatterplot, align = "h")[[1]]
  
  # Arrange plots
  plot_grid(
    aligned_x_hist,
    NULL,
    scatterplot,
    aligned_y_hist,
    ncol = 2,
    nrow = 2,
    rel_heights = c(0.15, 1),
    rel_widths = c(1, 0.15)
  )

}

df1 <- quick.df(round(10^(1.5)),round(10^(4.0)));

p1 <- quick.margin.plot(df1, alpha=0.2, size=1,trim.outliers=0.03)

Loading required package: cowplot

Attaching package: ‘cowplot’

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

    get_legend

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

    ggsave

Loading required package: ggrastr

p1

p2 <- quick.margin.plot(df2, alpha=0.1, size=0.5,trim.outliers=0.03)
p2

p3 <- quick.margin.plot(df3, alpha=0.01, size=0.1,trim.outliers=0.001)
p3

pdf(file='ex1.pdf',width=2,height=2); print(p1); dev.off()

null device 
          1 

pdf(file='ex2.pdf',width=2,height=2); print(p2); dev.off()

null device 
          1 

pdf(file='ex3.pdf',width=2,height=2); print(p3); dev.off()

null device 
          1