scripts/Lewy2b_removebias.R

In vladpetyuk/LewyBodies.SN.BottomUp2: Collection of bottom-up LC-MS/MS data for Lewy bodies project

#' ---
#' title: Removing plate column and patient match effects
#' author:
#' date:
#' output:
#'      html_document:
#'          highlight: default
#'          theme: default
#' ---

#+ init, echo=FALSE
set.seed(1)
library(knitr)
opts_chunk$set(message=FALSE, warning=FALSE, echo=FALSE,
               fig.path = "figure/", fig.width = 8, fig.height = 4,
               fig.align='center')

suppressPackageStartupMessages(library(vp.misc))
suppressPackageStartupMessages(library(lme4))
suppressPackageStartupMessages(library(multcomp))
suppressPackageStartupMessages(library(parallel))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(qvalue))
suppressPackageStartupMessages(library(reshape2))
suppressPackageStartupMessages(library(LewyBodies.SN.BottomUp))




#' # Lewy2b
data(Lewy2b_maxquant) # 3480
m <- m[rowSums(!is.na(exprs(m))) > 25,] # > half -> 2963 (15% loss)
m <- normalizeByGlob(m) # median polish normalization 
dSet <- cbind(pData(m)[,c('subject.type','PlateCol','match.group')], t(exprs(m)))
res <- mclapply(featureNames(m), function(x){
    fmla1 <- sprintf("`%s` ~ subject.type + (1|PlateCol) + (1|match.group)",x)
    m1 <- lmer(as.formula(fmla1), data = dSet)
}, mc.cores = 8)



col.eff <- t(sapply(res, function(x){coef(x)$PlateCol[m$PlateCol,1]}))
mch.eff <- t(sapply(res, function(x){coef(x)$match.group[m$match.group,1]}))

exprs(m) <- exprs(m) - col.eff
save(m, file = "Lewy2b_maxquant_no_col_bias.RData")

exprs(m) <- exprs(m) - mch.eff
save(m, file = "Lewy2b_maxquant_no_full_bias.RData")