demo_dreamlet_sims.Rmd

---
title: "Null simulations subsets of NPS/AD"
subtitle: 'STAR-solo R1-4 merged h5ad'
author: "Developed by [Gabriel Hoffman](http://gabrielhoffman.github.io/)"
date: "Run on `r Sys.time()`"
documentclass: article
output: 
  html_document:
  toc: true
  smart: false
vignette: >
  %\VignetteIndexEntry{Decorrelate}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
  %\usepackage[utf8]{inputenc}
---


<!--- 

# results 
# /sc/arion/projects/CommonMind/hoffman/NPS-AD/work/nps_ad/demo/results
cd /sc/arion/projects/CommonMind/hoffman/NPS-AD/work/nps_ad/
ml python git
git pull
R
# rm -rf demo_dreamlet_sims_cache

system("git pull"); rmarkdown::render("demo_dreamlet_sims.Rmd");


# https://hoffmg01.u.hpc.mssm.edu/nps_ad/demo_dreamlet_sims.html


ml git
cd ~/build2/dreamlet
git pull
R CMD INSTALL .


--->


- Simulate only 1 cohort
  - problem with user vs elapsed times
  - balanced vs unbalanced
  - use *exact* same filtering criteria
  - scale up simulations



*****Idea for creating confounding: need to random flip probs*******
- simulate with 1 gene



```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = FALSE,
  warning=FALSE,
  message=FALSE,
  error = FALSE,
  tidy = FALSE,
  dev = c("png", "pdf"),
  cache = TRUE,
  cache.lazy = FALSE)
```

```{r load.packages, cache=FALSE}
# Use cache=FALSE so that package are fully loaded each time
# This ensures that forks within mclapply() have these loaded
# Othewise, mclapply() not have access to these libraries and will fail 
#   unless the libraries are manually loaded within each fork
suppressPackageStartupMessages({
library(SingleCellExperiment)
library(dreamlet)
library(muscat)
library(scater)
library(tidyverse)
library(ggplot2)
library(gtools)
library(DESeq2)
library(variancePartition)
library(parallel)
library(edgeR)
library(lme4)
library(kableExtra)
library(RhpcBLASctl)
omp_set_num_threads(1)
blas_set_num_threads(1)
})
```

```{r load.data}
# Completely collapsed
######################
pbObj.Donor = readRDS("/sc/arion/projects/CommonMind/hoffman/NPS-AD/work/demo_dreamlet_pbobj_Donor.RDS")
pbObj.Donor$SubID = pbObj.Donor$Row.names 

# only MSSM and only batches of >=4 samples
pbObj.Donor = pbObj.Donor[,pbObj.Donor$cohort == "MSSM"]
tab = table(pbObj.Donor$BatchID)
pbObj.Donor = pbObj.Donor[,pbObj.Donor$BatchID %in% names(tab)[tab >=4]]
pbObj.Donor$BatchID = droplevels(pbObj.Donor$BatchID)

# includes replicates
pb = readRDS("/sc/arion/projects/CommonMind/hoffman/NPS-AD/work/demo_dreamlet_pbobj.RDS")
pb = pb[,pb$cohort == "MSSM"]
pb = pb[,pb$BatchID %in% names(tab)[tab >=4]]
pb$BatchID = droplevels(pb$BatchID)
pb$SubID = droplevels(pb$SubID)

# Full Data
outfolder = '/sc/arion/projects/CommonMind/hoffman/NPS-AD/work/nps_ad/demo/results/demo_dreamlet_sce_mod.RDS'
sceFull = readRDS(outfolder)
sceFull = sceFull[,sceFull$SubID %in% pb$SubID]
sceFull$BatchID = droplevels(sceFull$BatchID)
sceFull$SubID = droplevels(sceFull$SubID)

sim.method = "unbal"
# sim.method = "balanced"

# simulate balanced phenotype
set.seed(2)
ids_ad = c()

for( batch in levels(pb$BatchID)){

  # get Subject IDs in this batch
  ids = strsplit(batch, ',')[[1]]

  # keep only ids in the dataset
  ids = ids[ids %in% pb$SubID]

  if(sim.method == "balanced"){
    # select number of samples to be AD
    n = length(ids)

    # if not even, randomly add -0.5 or 0.5
    if( length(ids) %% 2 != 0){
      n = n/2 + sample(c(-0.5, 0.5), 1)
    }else{
      n = n/2
    }

    ids.select = ids[1:n]
  }else{
    # p = 0.9
    # prob = sample(c(p, 1-p), 1)
    # prob = c(prob, 1-prob)

    # ids.select = ids[which(sample(0:1, length(ids), replace=TRUE, prob=prob)==1)]

    n = sample(c(1, length(ids)-1), 1)

    ids.select = ids[1:n]
  }

  ids_ad = append(ids_ad, ids.select)
}

pb$Dxsim = rep("Control", ncol(pb))
pb$Dxsim[pb$SubID %in% ids_ad] = "AD"
pb$Dxsim = factor(pb$Dxsim, c("Control", "AD"))

pbObj.Donor$Dxsim = rep("Control", ncol(pbObj.Donor))
pbObj.Donor$Dxsim[pbObj.Donor$SubID %in% ids_ad] = "AD"
pbObj.Donor$Dxsim = factor(pbObj.Donor$Dxsim, c("Control", "AD"))

sceFull$Dxsim = rep("Control", ncol(sceFull))
sceFull$Dxsim[sceFull$SubID %in% ids_ad] = "AD"
sceFull$Dxsim = factor(sceFull$Dxsim, c("Control", "AD"))

# set cohort
sceFull$cohort = c()
for(CHRT in unique(pbObj.Donor$cohort)){

  ids = pbObj.Donor$SubID[pbObj.Donor$cohort == CHRT]

  sceFull$cohort[sceFull$SubID %in% ids] = CHRT
}
```

load('/sc/arion/scratch/hoffmg01/dreamlet/data_balanced_1.RDATA')
load('/sc/arion/scratch/hoffmg01/dreamlet/data_unbal_1.RDATA')

Evaluate balance
```{r balance}
canCorPairs( ~ Dxsim + BatchID, colData(pbObj.Donor))

fit = glm(Dxsim ~ BatchID, data=colData(pbObj.Donor), binomial())
calcVarPart(fit)
summary(fit)

tab = xtabs(~ BatchID + Dxsim, colData(pbObj.Donor))
frac = tab[,1] / rowSums(tab)
mean(frac)
hist(frac, xlab="AD fraction")
```

#### variancePartition
```{r vp}
res.proc = processAssays( pbObj.Donor, ~ 1, min.samples=6)

vp.lst = fitVarPart( res.proc, ~ (1|Dxsim) + (1|BatchID), BPPARAM=SnowParam(12))
```

```{r vp.plt, fig.height=12, fig.width=8}
# Summarize variance fractions genome-wide for each cell type
plotVarPart(vp.lst, label.angle=60)
```


```{r vp.plt2, fig.height=12, fig.width=8}
vp.lst2 = fitVarPart( res.proc, ~ Dxsim)

# Summarize variance fractions genome-wide for each cell type
plotVarPart(vp.lst2, label.angle=60)
```


```{r define.methods}
run_method = function(pb, formula, method){

  if( method == "dreamlet"){
    res.time = system.time({
      res.proc = processAssays( pb, ~ 1, BPPARAM=SnowParam(4), quiet=TRUE, min.samples=6) 
      fit = dreamlet( res.proc, formula, BPPARAM=SnowParam(4), quiet=TRUE)
      tab <- topTable(fit, coef='DxsimAD', number=Inf)
    })
  }else if( method %in% c("limma-voom", "edgeR", "DESeq2")){

      # convert counts to integer
      for(CT in assayNames(pb)){
        assay(pb, CT) = as.matrix( assay(pb, CT))
      }

      df = data.frame(sample_id = pb$SubID, 
                      group_id = pb$Dxsim) 

      metadata(pb)$experiment_info = df
      pb$group_id = df$group_id
      pb$sample_id = df$sample_id

      design = model.matrix(formula, droplevels(colData(pb)))
      colnames(design) = gsub(',', '_', colnames(design))

      metadata(pb)$agg_pars$by = c("cluster_id", "sample_id" )

      contrast = suppressWarnings(makeContrasts("DxsimAD", levels = design))

      res.time = system.time({
        res = pbDS(pb, method=method, design = design, contrast=contrast, verbose=FALSE)

        tab = lapply( names(res$table$DxsimAD), function(CT){
          tab = res$table$DxsimAD[[CT]]
          tab$assay = CT
          tab$ID = tab$gene
          tab$cluster_id = CT
          tab$AveExpr = tab$logCPM 
          tab$P.Value = tab$p_val 
          tab
          })
        tab = do.call(rbind, tab)
        tab$adj.P.Val = p.adjust(tab$P.Value, "BH")
      })
  }else if( method == "glmer"){

    # filter genes and samples using dreamlet code
    res.proc = processAssays( pb, ~ 1, quiet=TRUE, min.samples=6)

    sce = sceFull[,sceFull$SubID %in% unique(pb$SubID)]
  
    res.time = system.time({
    # DESeq2
    res = lapply(assayNames(res.proc), function(CT){
        
        # get included genes and samples
        tmp = assay(res.proc, CT)

        # cat(CT, dim(tmp))

        # subset based on dreamlet filter
        pb.sub = pb[rownames(tmp), colnames(tmp)]

        sce.CT = sce[,sce$celltype8==CT]
        df = droplevels(as.data.frame(colData(sce.CT)))

        counts = as(assay(sce.CT, "X"), "sparseMatrix")
        rm(sce.CT)

        df$cs = colSums(counts)

        res.glmer = mclapply(rownames(pb.sub), function(gene){

          df$y = c(counts[gene,])
          form = as.formula(paste('y ~ offset(log(cs)) + ', as.character(formula)[2] ,'+ (1|SubID)'))

          fit = suppressMessages(glmer(form, data=df, family = 'poisson'))

          data.frame(ID = gene, t(coef(summary(fit))['DxsimAD',]))
          }, mc.cores=4)
        res.glmer = do.call(rbind, res.glmer)
        rm(counts)

        data.frame(ID = res.glmer$ID,
            cluster_id = CT,
            logFC = res.glmer$Estimate,
            AveExpr = NA, 
            t = res.glmer$z.value, 
            P.Value = res.glmer$Pr...z.., 
            p_adj.loc = p.adjust(res.glmer$Pr...z.., "BH"))
    })
    tab = do.call(rbind, res)
    tab$adj.P.Val = p.adjust(tab$P.Value, "BH")
    })    
  }else{
    stop("method not found:", method)
  }

  tab$n = ncol(pb)
  tab$method = method
  tab$formula = paste(as.character(formula), collapse=' ')
  tab$time.user = res.time[1]
  tab$time.system = res.time[2]
  tab$time.elapsed = res.time[3]

  tab
}
```



```{r run.sims}
run_analysis = function(n_batches, methods){

  resTabl = mclapply(n_batches, function(nbatch){

    idx = pbObj.Donor$BatchID %in% levels(pbObj.Donor$BatchID)[1:nbatch]

    resTabl = lapply(methods, function(method){
      message(method, ' ', nbatch )
      res1 = run_method( pbObj.Donor[1:1000,idx], ~ Dxsim, method)
      res2 = run_method( pbObj.Donor[1:1000,idx], ~ Dxsim + BatchID, method)
      # if( method %in% c("dreamlet", "glmer")){      
      #   res3 = run_method( pbObj.Donor[1:1000,idx], ~ Dxsim + (1|BatchID), method)
      # }else{
        res3 = data.frame()
      # }
      gc()
      smartbind(res1, res2, res3)
    })
    smartbind(list=resTabl)
  }, mc.cores=6)
  gc()
  resTabl = do.call(rbind, resTabl)
}
 
n_batches = c(2, 3, 4)
methods = c('dreamlet', "limma-voom")#, "DESeq2", "edgeR", "glmer")
resTabl1 = run_analysis(n_batches, methods)

n_batches = c(6, 12, 17, 29)#, 41, 50, 68)
methods = c('dreamlet', "limma-voom")#, "DESeq2", "edgeR")
resTabl2 = run_analysis(n_batches, methods)

resTabl = rbind(resTabl1, resTabl2)
```

Note that run time includes disk access to h5ad file
```{r plot.time, cache=FALSE, fig.width=8}
methods = c("dreamlet", "limma-voom", "edgeR", "DESeq2", "glmer")

xmax = 1.02* max(resTabl$n)
form_labels = c("~ Dxsim", "~ Dxsim + BatchID", "~ Dxsim + (1 | BatchID)")

resTabl %>% 
  as_tibble %>% 
  group_by(n, method, formula) %>% 
  summarize( minutes = max(time.elapsed)/60) %>%
  mutate(method = factor(method, unique(methods)),
        formula = factor(formula, form_labels)) %>%
  ggplot(aes(n, minutes, color=method, formula)) + 
    geom_line() + 
    geom_point() + 
    theme_classic() +
    theme(aspect.ratio=1, plot.title = element_text(hjust = 0.5)) + 
    scale_y_log10() + 
    scale_x_continuous(limits=c(0, xmax), expand=c(0,0)) +
    ggtitle("Run time") + 
    xlab("# Subjects") + 
    ylab("Elapsed time (min)") +
    scale_color_brewer(palette="Set1", name="Method") +
    facet_wrap(~formula)
```


```{r plot.time.user, cache=FALSE, fig.width=8}
xmax = 1.02* max(resTabl$n)

resTabl %>% 
  group_by(n, method, formula) %>% 
  summarize( minutes = max(time.user)/60) %>% 
  mutate(method = factor(method, unique(methods)),
        formula = factor(formula, form_labels)) %>%
  ggplot(aes(n, minutes, color=method)) + 
    geom_line() + 
    geom_point() + 
    theme_classic() +
    theme(aspect.ratio=1, plot.title = element_text(hjust = 0.5)) + 
    scale_y_log10() + 
    scale_x_continuous(limits=c(0, xmax), expand=c(0,0)) +
    ggtitle("CPU time") + 
    xlab("# Subjects") + 
    ylab("CPU time (min)") +
    scale_color_brewer(palette="Set1", name="Method") +
    facet_wrap(~formula)
```


```{r plot.FD, cache=FALSE, fig.width=8}
# retain fixed set of assays and genes
df = resTabl %>%
  as_tibble %>% 
  group_by(assay, ID) %>% 
  count()
maxValue = max(df$n)
df_keepPairs = df %>%
  filter(n==maxValue) %>%
  mutate(pair = paste(assay, ID)) 

df = resTabl %>%
  as_tibble %>%
  filter(paste(assay, ID) %in% df_keepPairs$pair) %>%
  group_by(n, method, formula) %>%
  summarize( nDE = sum(adj.P.Val < 0.05, na.rm=TRUE), total = sum(!is.na(adj.P.Val))) 
ymax = 1.02* max(df$nDE)
xmax = 1.02* max(df$n)

df %>%
  mutate(method = factor(method, unique(methods)),
        formula = factor(formula, form_labels)) %>%
  ggplot(aes(n, nDE, color=method)) + 
    geom_point() + 
    theme_classic() +
    theme(aspect.ratio=1, plot.title = element_text(hjust = 0.5)) + 
    scale_y_continuous(limits=c(0, ymax)) + 
    scale_x_continuous(limits=c(0, xmax), expand=c(0,0)) +
    ggtitle("False positives from simulations") + 
    xlab("# Subjects") + 
    ylab("False positive genes") +
    scale_color_brewer(palette="Set1", name="Method") +
    facet_wrap(~formula)
```



```{r plot.FPR, cache=FALSE, fig.width=8}
df = resTabl %>%
  as_tibble %>%
  filter(paste(assay, ID) %in% df_keepPairs$pair) %>%
  group_by(n, method, formula) %>%
  summarize( FPR = sum(P.Value < 0.05, na.rm=TRUE) / sum(!is.na(P.Value)))

ymax = 1.02* max(df$FPR)
xmax = 1.02* max(df$n)

df %>%
  mutate(method = factor(method, unique(methods)),
        formula = factor(formula, form_labels)) %>%
  ggplot(aes(n, FPR, color=method)) + 
    geom_hline(yintercept = 0.05, linetype="dashed", color="grey") +
    geom_point() +
    theme_classic() +
    theme(aspect.ratio=1, plot.title = element_text(hjust = 0.5)) + 
    scale_y_continuous(limits=c(0, ymax), expand=c(0,0)) + 
    scale_x_continuous(limits=c(0, xmax), expand=c(0,0)) +
    ggtitle("False positives from simulations") + 
    xlab("# Subjects") + 
    ylab("False positive rate") +
    scale_color_brewer(palette="Set1", name="Method") +
    facet_wrap(~formula)
```



```{r plot.FPR.bar, cache=FALSE, fig.width=8}
df = resTabl %>%
  as_tibble %>%
  # filter(paste(assay, ID) %in% df_keepPairs$pair) %>%
  group_by(n, method, formula) %>%
  summarize( FPR = sum(P.Value < 0.05, na.rm=TRUE) / sum(!is.na(P.Value))) %>%
  as_tibble %>% 
  tidyr::complete(n, method, formula, fill=list(FPR=0))

ymax = 1.02* max(df$FPR)
xmax = 1.02* max(df$n)

df %>%
  mutate(method = factor(method, unique(methods)),
        formula = factor(formula, form_labels)) %>%
  ggplot(aes(factor(n), FPR, fill=method)) + 
    geom_bar(stat="identity", position="dodge2") +
    geom_hline(yintercept = 0.05, linetype="dashed", color="grey") +
    theme_classic() +
    theme(aspect.ratio=1, plot.title = element_text(hjust = 0.5)) + 
    scale_y_continuous(limits=c(0, ymax), expand=c(0,0)) + 
    ggtitle("False positive rate from simulations") + 
    xlab("# Subjects") + 
    ylab("False positive rate") +
    scale_fill_brewer(palette="Set1", name="Method") +
    facet_grid(~formula)
```