Median-balanced quantile normalization

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----bc-installation, eval = FALSE--------------------------------------------
#  # if (!requireNamespace("BiocManager", quietly = TRUE))
#  #    install.packages("BiocManager")
#  BiocManager::install("MBQN")

## ----dependencies1, eval = FALSE----------------------------------------------
#  # if (!requireNamespace("BiocManager", quietly = TRUE))
#  #    install.packages("BiocManager")
#  # BiocManager::install(pkgs = c("preprocessCore","limma","rpx","SummarizedExperiment"))

## ----example1, eval = TRUE----------------------------------------------------
## basic example
library("MBQN")
set.seed(1234)
# data generation 
mtx <- mbqnSimuData("omics.dep")
# data distortion
mtx <- mbqnSimuDistortion(mtx)$x.mod

## ----figure1, fig.height = 5, fig.width = 6, fig.align = "left", fig.cap = "Fig. 1 Boxplot of the unnormalized, distorted intensity data matrix. The first feature is an RI feature (red line). It has maximum intensity for each sample!"----
plot.new()
mbqnBoxplot(mtx, irow = 1, main = "Unnormalized")

## ---- eval = TRUE-------------------------------------------------------------
res <- mbqnGetNRIfeatures(mtx, low_thr = 0.5)

## ----figure2, fig.height = 5, fig.width = 6, fig.align = 'left', fig.cap = "Fig. 2 Quantile normalized intensities with balanced and unbalanced normalized RI feature. Classical quantile normalization suppresses any intensity variation of the RI feature, while the MBQN preserves its variation while reducing systematic batch effects!"----
plot.new()
mbqn.mtx <- mbqnNRI(x = mtx, FUN = median, verbose = FALSE) # MBQN
qn.mtx <- mbqnNRI(x = mtx, FUN = NULL, verbose = FALSE) # QN
mbqnBoxplot(mbqn.mtx, irow = res$ip, vals = data.frame(QN = qn.mtx[res$ip,]), main = "Normalized")

## ----example2, eval = TRUE----------------------------------------------------
## basic example
library("MBQN")
set.seed(1234)
mtx <- mbqnSimuData("omics.dep") 
# Alternatively: mtx <- matrix(
#            c(5,2,3,NA,2,4,1,4,2,3,1,4,6,NA,1,3,NA,1,4,3,NA,1,2,3),ncol=4)

## ---- eval = TRUE-------------------------------------------------------------
qn.mtx <- mbqn(mtx,FUN=NULL, verbose = FALSE)
mbqn.mtx <- mbqn(mtx,FUN = "median", verbose = FALSE)
qn.nri.mtx <- mbqnNRI(mtx,FUN = "median", low_thr = 0.5, verbose = FALSE)

## ---- eval = TRUE-------------------------------------------------------------
res <- mbqnGetNRIfeatures(mtx, low_thr = 0.5)
# Maximum frequency of RI/NRI feature(s):  100 %

## ----example3, eval = TRUE----------------------------------------------------
#plot.new()
mtx <- mbqnSimuData("omics.dep", show.fig = FALSE)
mod.mtx <- mbqnSimuDistortion(mtx, s.mean = 0.05, s.scale = 0.01)
mtx2 <- mod.mtx
mod.mtx <- mod.mtx$x.mod

res <- mbqnGetNRIfeatures(mod.mtx, low_thr = 0.5)

# undistorted feature
feature1 <- mtx[1,]
# distorted feature
feature1mod = mod.mtx[1,]
# feature after normalization
qn.feature1 = mbqn(mod.mtx, verbose = FALSE)[1,]
qn.mtx = mbqn(mod.mtx,verbose = FALSE)

mbqn.mtx = mbqn(mod.mtx, FUN = "mean",verbose = FALSE)
mbqn.feature1 = mbqn(mod.mtx, FUN = "mean",verbose = FALSE)[1,]

## ---- eval = TRUE-------------------------------------------------------------
# undistorted feature
ttest.res0 <- t.test(feature1[seq_len(9)], feature1[c(10:18)],
                    var.equal =TRUE)
# distorted feature
ttest.res1 <- t.test(feature1mod[seq_len(9)], feature1mod[c(10:18)],
                    var.equal =TRUE)
# mbqn normalized distorted feature
ttest.res <- t.test(mbqn.feature1[seq_len(9)], mbqn.feature1[c(10:18)],
                    var.equal =TRUE)

## ---- eval = TRUE-------------------------------------------------------------

## ----figure3, fig.height = 5, fig.width = 6, fig.align = "left", fig.cap = "Fig. 3 "----
plot.new()
matplot(t(rbind(feature1 = feature1,
    mod.feature1 = (feature1mod-mean(feature1mod))/25+mean(feature1),
    qn.feature1 = (qn.feature1-mean(qn.feature1))+mean(feature1),
    mbqn.feature1 = (
        mbqn.feature1-mean(mbqn.feature1))+mean(feature1))),
    type = "b", lty = c(1,1,1), pch = "o",
    ylab = "intensity",
    xlab = "sample",
    main = "Differentially expressed RI feature",
    ylim = c(34.48,34.85))
legend(x=11,y= 34.86, legend = c("feature","distorted feature/25" ,
                                "QN feature", " MBQN feature"),pch = 1,
        col = c(1,2,3,4), lty= c(1,1,1,1), bty = "n", y.intersp = 1.5,
        x.intersp = 0.2)
legend(x = .1, y = 34.6,
        legend = paste("p-value (t-test) =",round(ttest.res1$p.value,2),
            "\np-value (t-test, mbqn) =", round(ttest.res$p.value,4)),
        bty = "n", x.intersp = 0)


if (ttest.res$p.value<0.05)
    message("H0 (=equal mean) is rejected!")


# print(mtx2$x.mod)
# print(mtx2$mx.offset)
# print(mtx2$mx.scale)
print(paste("ttest.undistorted =",ttest.res0)) 
print(paste("ttest.distorted =", ttest.res1))
print(paste("ttest.mbqndistorted =", ttest.res))

## ----example4, eval = TRUE----------------------------------------------------
library("SummarizedExperiment")
pxd_id <- "PXD001584"   

# Download file from PRIDE to currentdir/PXDxxx
#getPXDfile(pxd_id)

# Load file
out <- mbqnLoadFile(pxd_id, file.pattern = "proteinGroups.txt")

# filter for potential contaminants and identified only by site features
out <- out[!rowData(out)[["ixs"]],] 

# extract data and feature annotation
mtx <- assays(out)[["data"]]
featureAnnotations <- rowData(out)

low_thr <- 0.5
ylim <- NULL
ix <- seq_len(ncol(mtx))

ix <- c(1:9,19:27)
mtx <- mtx[,ix]
ylim.qn <- ylim <- c(22.5,36)

####################################################################################

## ----figure4, fig.height = 3, fig.width = 6, fig.align = "left", fig.cap = "Fig. 4 Correctly detected and identified RI (100 %RI) and a NRI feature (75 %RI) each with a data coverage of 100% across samples. "----
plot.new()
res <- mbqnPlotAll(mtx,
                   FUN = median,
                   low_thr = low_thr,
                   las = 2,
                   type = "l",
                   feature_index = NULL,
                   show_nri_only = TRUE,
                   axis.cex = 0.5,
                   y.intersp= 0.5)
#dev.off()

# get protein name of strongest nri/ri feature
nri_max <- as.numeric(names(which.max(res$nri)))
#featureAnnotations$proteinDescription[nri_max]
#featureAnnotations$proteinName[nri_max]
#featureAnnotations$nbPeptides[nri_max]

colnames(mtx) <- gsub("LFQ intensity","",colnames(mtx))
mbqn.mtx <- mbqn(mtx,FUN = median)
qn.mtx <- mbqn(mtx,FUN = NULL)

# Boxplot of QN intensity features, highlight RI/NRI Features
if(length(ylim)==0) {
  ylim <- c(floor(min(range(mbqn.mtx, na.rm = TRUE))),ceiling(max(range(mbqn.mtx, na.rm = TRUE))))
  ylim.qn <- c(floor(min(range(qn.mtx, na.rm = TRUE))),ceiling(max(range(mbqn.mtx, na.rm = TRUE))))
}

df <- data.frame(qn.mtx[as.numeric(names(res$nri)),])
if(ncol(df)==1) df <- t(df)
rownames(df) <- paste("QN feature",names(res$nri))
df2 <- data.frame(mtx[as.numeric(names(res$nri)),])
if(ncol(df2)==1) df2 <- t(df2)
rownames(df2) <- paste("unnormal. feature",names(res$nri))
colnames(df) <- colnames(df2)
df <- rbind(df,df2)
df <- as.data.frame(t(df))

## ----figure5, fig.height = 5, fig.width = 6, fig.align = "left", fig.cap = "Fig. 5 "----
  mtx.nri <- mbqnNRI(mtx,FUN = median,low_thr = 0.5, verbose = FALSE)

plot.new()
mbqnBoxplot(mtx = mtx.nri,
            irow = as.numeric(names(res$nri)),
            ylim = ylim.qn, xlab= "", las=2,
            ylab = "LFQ intensity",
            vals = df,lwd = 1.,
            main = "QN with RI/NRI balanced",
            cex.axis = 1, cex.lab = .9, cex = .9, y.intersp = 0.5)