dream_DE.R

# load library
library('variancePartition')
# load simulated data:
# geneExpr: matrix of gene expression values
# info: information/metadata about each sample
data(varPartData)
form <-  ~ progression + (1 | slide)
browser()
varPart <- fitExtractVarPartModel( geneExpr, form, info )

load("C:/Users/edmondsonef/Desktop/DSP GeoMx/KPC_geoMX.RData")


# convert test variables to factors
pData(target_myData)$progression <- 
  factor(pData(target_myData)$prog4)                        
pData(target_myData)[["slide"]] <-                                
  factor(pData(target_myData)[["MHL Number"]])
assayDataElement(object = target_myData, elt = "log_q") <-
  assayDataApply(target_myData, 2, FUN = log, base = 2, elt = "q_norm")

# run LMM:
# formula follows conventions defined by the lme4 package
results <- c()
for(status in c("Full ROI")) {
  ind <- pData(target_myData)$segment == status
  mixedOutmc <-
    mixedModelDE(target_myData[, ind],
                 elt = "log_q",
                 
                 modelFormula = ~ progression + (1 | slide),
                 #modelFormula = ~ testRegion + (1 + testRegion | slide),
                 groupVar = "progression",
                 nCores = parallel::detectCores(),
                 multiCore = FALSE)
  
  # format results as data.frame
  r_test <- do.call(rbind, mixedOutmc["lsmeans", ])
  tests <- rownames(r_test)
  r_test <- as.data.frame(r_test)
  r_test$Contrast <- tests
  
  # use lapply in case you have multiple levels of your test factor to
  # correctly associate gene name with it's row in the results table
  r_test$Gene <- 
    unlist(lapply(colnames(mixedOutmc),
                  rep, nrow(mixedOutmc["lsmeans", ][[1]])))
  r_test$Subset <- status
  r_test$FDR <- p.adjust(r_test$`Pr(>|t|)`, method = "fdr")
  r_test <- r_test[, c("Gene", "Subset", "Contrast", "Estimate", 
                       "Pr(>|t|)", "FDR")]
  results <- rbind(results, r_test)
}