In satabdisaha1288/scBT: What the Package Does (One Line, Title Case)

suppressPackageStartupMessages({
  library("splattdr")
  library("snseq.stats")
  library("dplyr")
  library("DescTools")
  library("ggplot2")
  library("MAST")
})

data(params)

sim.dose = c(0, 0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30)
dose.prob = rep(1/9, 9)

params@nGenes = 900

parameter.updates <- read.table('ParTestsCombos.txt', header = TRUE, sep = '\t')

# Create empty lists to append and save outputs
DGE1 = list()
B.LIST = list()
pHH0 = list()
timing.list = list()
fc.calc = list()
pz.calc = list()
de_idx = list()

for (n in 1:nrow(parameter.updates)){
  print(n)
  params@nBatches = length(sim.dose)*3
  params@batch.facLoc = rep(0.1, params@nBatches)
  params@batch.facScale = rep(0.1, params@nBatches)
  params@de.facLoc = as.numeric(parameter.updates$FacLoc[n])
  params@de.facScale = as.numeric(parameter.updates$FacScale[n])
  cellNums = as.numeric(unlist(strsplit(as.character(parameter.updates$CellN[n]), split = ',')))
  if (length(cellNums) == 1){
    params@nCells = cellNums*9
    dose.prob = rep(1/9, 9)
  } else {
    params@nCells = sum(cellNums)
    dose.prob = cellNums/sum(cellNums)
  }
  params@de.prob = as.numeric(parameter.updates$deProb[n])
  params@de.downProb = as.numeric(parameter.updates$deDown[n])
  params@seed = as.numeric(parameter.updates$seed[n])
  libscale = as.numeric(parameter.updates$libScale[n])

  sim = splatSimulateDR(params, dose.names = sim.dose, dose.prob = dose.prob, lib.scale = libscale, verbose = FALSE)
  print(sim)
  sim <- sim[rowSums(logcounts(sim))>0,]
  print(sim)
  DGE1[[parameter.updates$TestName[n]]] = DETest(sim, method = c("All"), fixed.priors = FALSE, verbose = TRUE)

  ######## MAST #####
  timing <- system.time({
    sce = sim
    scaRaw <- FromMatrix(as.matrix(logcounts(sce)),
                        data.frame(colData(sce)),
                        data.frame(rowData(sce))
    )
    cdr <- colSums(assay(scaRaw) > 0)
    colData(scaRaw)$cngeneson <- scale(cdr)
    rowData(scaRaw)$symbolid <- rownames(rowData(sce))
    Dose <- factor(colData(scaRaw)$Dose)
    Dose <- relevel(Dose,"0")
    colData(scaRaw)$Dose <- Dose
    #zlmDose <- zlm(~Dose + cngeneson, sca = scaRaw)
    zlmDose <- zlm(~Dose, sca = scaRaw)

    summaryDose <- rep(list(list()), times = nlevels(Dose)-1)
    names(summaryDose) <- levels(Dose)[-1]
    for(i in levels(Dose)[-1])
    {
      summaryDose[[i]] <- summary(zlmDose, doLRT=paste0("Dose",i))
    }

    summaryDT <- rep(list(data.table::data.table()),times=nlevels(Dose)-1)
    names(summaryDT) <- levels(Dose)[-1]
    fcHurdle <- rep(list(data.table::data.table()),times=nlevels(Dose)-1)
    names(fcHurdle) <- levels(Dose)[-1]
    fcHurdleSig <- rep(list(data.table::data.table()),times=nlevels(Dose)-1)
    names(fcHurdleSig) <- levels(Dose)[-1]
    fcHurdleNonSig <- rep(list(data.table::data.table()),times=nlevels(Dose)-1)
    names(fcHurdleSig) <- levels(Dose)[-1]
    TP_dose <- rep(list(data.table::data.table()),times=nlevels(Dose)-1)
    names(TP_dose) <- levels(Dose)[-1]
    FP_dose <- rep(list(data.table::data.table()),times=nlevels(Dose)-1)
    names(FP_dose) <- levels(Dose)[-1]
    TN_dose <- rep(list(data.table::data.table()),times=nlevels(Dose)-1)
    names(TN_dose) <- levels(Dose)[-1]
    FN_dose <- rep(list(data.table::data.table()),times=nlevels(Dose)-1)
    names(FN_dose) <- levels(Dose)[-1]

    for(i in levels(Dose)[-1]){
      summaryDT[[i]] <- data.frame(summaryDose[[i]]['datatable'])
      fcHurdle[[i]] <- summaryDT[[i]] %>% 
        filter(datatable.contrast == paste0("Dose",i) & datatable.component == 'H') %>% 
        select(datatable.primerid, datatable.Pr..Chisq.)
      fcHurdle[[i]]$FDR <- p.adjust(fcHurdle[[i]]$datatable.Pr..Chisq., 'fdr')
    }

    mast.out <- do.call(cbind, lapply(fcHurdle, data.frame))
    mast.final <- mast.out[,grepl('FDR', colnames(mast.out))]
    rownames(mast.final) <- mast.out[, 1]
    mast.final$adjusted.p <- apply(mast.final, 1, function(x) min(x))
  })
  DGE1[[parameter.updates$TestName[n]]][['MAST']] <- mast.final
  DGE1[[parameter.updates$TestName[n]]]$timing[['MAST']] <- timing
  ###### END MAST ####

  pHH0[[parameter.updates$TestName[n]]] <- DGE1[[parameter.updates$TestName[n]]]$ppH0
  timing.list[[parameter.updates$TestName[n]]] <- DGE1[[parameter.updates$TestName[n]]]$timing
  de_idx[[parameter.updates$TestName[n]]] <- rowData(sim)$DE_idx
  fc.calc[[parameter.updates$TestName[n]]] <- calcFC(sim)
  pz.calc[[parameter.updates$TestName[n]]] <- calcZeroP(sim)
  saveRDS(DGE1, file = 'DGE1.062221.RData')
  saveRDS(pHH0, file = 'ppH0.062221.RData')
  saveRDS(fc.calc, file = 'fc.calc.062221.RData')
  saveRDS(pz.calc, file = 'pz.calc.062221.RData')
  saveRDS(de_idx, file = 'de_idx.062221.RData')

  benchmark.list = list()
  NamesDGE <- names(DGE1[[parameter.updates$TestName[n]]])
  for (t in NamesDGE[which(NamesDGE != 'timing' & NamesDGE != 'ppH0')]){
    for (fc in c(0)){
      for (pct in c(0, 0.05)){
        name <- paste0(t, '_', fc, '_', pct)
        benchmark <- benchmarkSim_batch(sim, DGE1[[parameter.updates$TestName[n]]][[t]], fc.threshold = fc, pct.expressed = pct)
        summary <- summarizeBenchmark(benchmark)
        summary$test = t
        summary$fc = fc
        summary$pct = pct
        benchmark.list[[name]] = summary
        saveRDS(benchmark.list, file = 'benchmark.list.070221.RData')
      }
    }
  }

  B.LIST[[parameter.updates$TestName[n]]] = do.call('rbind', benchmark.list)
  saveRDS(B.LIST, file = 'B.LIST.062221.RData')
}