In lulizou/spASE: Detecting Allele-Specific Expression in Spatial Transcriptomics

library(dplyr)
library(cowplot)
library(ggplot2)
library(reshape2)
knitr::opts_chunk$set(echo = F, warning = F, message = F)

plot.ci.coverage <- function(df, title='') {
  df %>%
    group_by(n) %>%
    summarise(bb = sum(bb,na.rm=T)/sum(bb>=0,na.rm=T), glm = sum(glm,na.rm=T)/sum(glm>=0,na.rm=T), quasi = sum(quasi,na.rm=T)/sum(quasi>=0,na.rm=T)) %>%
    melt(id.vars = 'n') %>%
    ggplot(aes(x = n, y = value)) +
    geom_hline(yintercept=0.95, lty='dashed',color='red') +
    geom_line(aes(color = variable)) +
    scale_color_discrete(name='', labels=c('beta-binomial wald','binomial wald','quasibinomial wald')) +
    theme_bw() +
    ylab('CI coverage') +
    xlab('num cells or beads') +
    ggtitle(title)
}

plot.all <- function(df, title='') {
  p1 <- plot.ci.coverage(df$results.smartseq, title=paste0('Smart-seq3, ',title))
  p2 <- plot.ci.coverage(df$results.slideseq.low, title=paste0('Slide-seq Low, ',title))
  p3 <- plot.ci.coverage(df$results.slideseq.high, title=paste0('Slide-seq High, ',title))
  # make all in a row have the same y range
  # range1 <- layer_scales(p1)$y$range$range
  # range2 <- layer_scales(p2)$y$range$range
  # range3 <- layer_scales(p3)$y$range$range
  # ylim1 <- min(c(range1, range2, range3))
  # ylim2 <- max(c(range1, range2, range3))
  prow <- plot_grid(
    p1 + theme(legend.position='none'),
    p2 + theme(legend.position='none'),
    p3 + theme(legend.position='none'),
    nrow = 1
  )
  legend_b <- get_legend(
    p1 + 
      guides(color = guide_legend(nrow = 1)) +
      theme(legend.position = "bottom")
  )
  print(
    plot_grid(prow, legend_b, ncol=1, rel_heights=c(1,.1))
  )
}

for (phi in seq(0,1,0.1)) {
  df <- readRDS(paste0('ci_lowrange2_mean0.5_phi',phi,'_iter5000.rds'))
  plot.all(df, title=paste0('phi=',phi))
}

lulizou/spASE documentation built on May 22, 2024, 5:24 a.m.