In PDiracDelta/CONSTANd: Data normalization by matrix raking
knitr::opts_chunk$set(echo = TRUE,dev="CairoPNG")
library(knitr)
library(tidyr)
library(ggplot2)
library(gridExtra)
library(magick)
library(limma)
Introduction
CONSTANd normalization has been proven to work on mass spectrometry intensity values (Maes et al.) in shotgun proteomics and on RNA-seq count data (Van Houtven et al.), but is in principle applicable to any kind of quantification matrices, for instance in other -omics disciplines. It normalizes the data matrix data by raking (using the RAS method by Bacharach, see references) the Nrows by Ncols matrix such that the row means and column means equal 1. The result is a normalized data matrix K=RAS, a product of row mulipliers R and column multipliers S with the original matrix A. Missing information needs to be presented as NA values and not as zero values, because CONSTANd is able to ignore missing values when calculating the mean.
Getting started
The CONSTANd package has to be loaded to use the function.
To load CONSTANd, type the following command after launching R.
#install.packages("BiocManager")
#BiocManager::install("CONSTANd")
library(CONSTANd)
Data assumptions {#s:assumptions}
To warrant that the normalization procedure itself is not biased, three assumptions about the data are made (as for any type of data-driven normalization method) which may be verified by making MA-plots (see examples in section \@ref(s:maplots)):
- The majority of features (peptides/genes/...) are not differentially expressed, to avoid a bias in the estimate of the mean value. It is assumed that up-down shifts are not due to biological causes. The reference set used in the normalization step is the set of all peptides identified in the experiment.
MA-plot: the observations form a single 'cloud' with a dense center and less dense edges, as opposed to, for instance, two clouds or a cloud with uniform density.
- The number of up-regulated features is roughly equal to the number of down-regulated ones. If the data were skewed, this would lead to a bias in the normalization result.
MA-plot: the cloud of observations exhibits a bilateral symmetry about some axis (usually horizontal, but inclinations may occur).
- Any systematic bias present is linearly proportional to the magnitude of the quantification values. Only this way it is possible to find one appropriate normalization factor for each quantification sample.
MA-plot: the axis of bilateral symmetry is a straight line (which may be inclined), i.e., the moving average M-values of the central cloud form an approximately straight line.
Example of a bad MA plot
The following MA plot violates all three assumptions. In case your MA plots exhibit any of such characteristics, it is not advised to use CONSTANd - or any other data-driven method - to normalize your data, and to instead use advanced model-based approaches or, if appropriate, clean your data first.
par(mfrow=c(1,1))
include_graphics('bad_MA.png')
Normalizing a data matrix (single assay)
The data matrix is an Nrows by Ncols matrix (e.g. an assay from a SummarizedExperiment), representing quantification values of any kind, as long as they form an identically processed (sub)set (IPS). An IPS consists of measurements from biological samples who have been processed from biological harvesting of the tissue to the detection and quantitation of the features of interest (RNA, proteins, ...) in near-identical fashion. In practice, this means that there exist no factors (except those of biological interest) which may cause a systematic bias in the quantification values. One _counter_example in proteomics is samples measured on a different day or a different instrument, one _counter_example in RNA-seq would be samples measured after using different library preparation protocols. In the next section, we demonstrate how in such such cases the IPSs are normalized separately as if they were individual assays (see this section), after which they may be re-combined.
First, let's load the 'Leishmania' mRNA-level data set from GEO (Gene Expression Omnibus) with identifier GSE95353. You can either download and unzip it yourself, are use the file that comes with this package.
leish <- read.csv("../inst/extdata/GSE95353_SL_SEQ_LOGSTAT_counttable.txt" , sep='\t')
We only need the LOG and STAT samples, which correspond to biological quadruplicates of Leishmania in two lifestages (LOG and STAT).
rownames(leish) <- leish$gene
leish <- leish[, c("SL.LOG1", "SL.LOG2", "SL.LOG3", "SL.LOG4",
"SL.STAT1", "SL.STAT2", "SL.STAT3", "SL.STAT4")]
conditions <- c(rep("LOG", 4), rep("STAT", 4))
Let's remove genes with low count numbers (median > 0) for good practice, and make sure that zero counts are replaced by NA to use CONSTANd.
# first set NA to 0 for calculating the row medians
leish[is.na(leish)] <- 0
leish <- leish[apply(leish, 1, median, na.rm = FALSE) > 0,]
leish.norm <- CONSTANd(leish)$normalized_data
Let's make a PCA plot (optionally impute NA as zero) before and after normalization to see if it was successful:
leish[is.na(leish)] <- 0
# scale raw data
leish.pc.raw <- prcomp(x=t(leish), scale. = TRUE)
leish.norm[is.na(leish.norm)] <- 0
# CONSTANd already scaled the data
leish.pc.norm <- prcomp(x=t(leish.norm))
par(mfrow=c(1,2))
colors = ifelse(conditions == 'LOG', 'blue', 'red')
plot(leish.pc.raw$x[,'PC1'], leish.pc.raw$x[,'PC2'],
xlab='PC1', ylab='PC2', main='raw counts', col=colors)
plot(leish.pc.norm$x[,'PC1'], leish.pc.norm$x[,'PC2'],
xlab='PC1', ylab='PC2', main='normalized counts', col=colors)
Indeed, the PCA plot indicates that the first principal component, which is the direction along which there is most varability in the measurements, successfully defines a separation between the biological conditions. In other words: after normalization the biological differences between the samples dominate any other sources of systematic bias.
Normalization across multiple assays
Here, we demonstrate multiple IPSs involved in one experiment may be normalized separately as if they were individual assays (see previous section), after which they may be re-combined for joint analysis.
We demonstrate CONSTANd normalization across multiple assays using the Organs example data set from QCQuan, a proteomics webtool built around the CONSTANd algorithm.
The Organs data set entails samples from 8 different mouse organs in biological quadruplicate, whose peptides have been measured in 4 different instrument runs (1 per mouse) using TMT-labeled (8-plex) tandem mass spectrometry.
par(mfrow=c(1,1))
include_graphics('organ_experiment_design.png')
You can either download the original data from the web, unzip it and set LOAD_ORIGINAL=TRUE, or load the cleaned data from this package.
LOAD_ORIGINAL = FALSE
# load extra functions to clean the original data and study design files
source('../inst/script/functions.R')
if (LOAD_ORIGINAL) { # load original from web
BR1 <- read.csv('Organs_PSMs_1.txt', sep='\t')
BR2 <- read.csv('Organs_PSMs_2.txt', sep='\t')
BR3 <- read.csv('Organs_PSMs_3.txt', sep='\t')
BR4 <- read.csv('Organs_PSMs_4.txt', sep='\t')
organs <- list(BR1, BR2, BR3, BR4)
organs <- clean_organs(organs)
} else { # load cleaned file that comes with this package
organs <- readRDS('../inst/extdata/organs_cleaned.RDS')
}
# construct the study design from the QCQuan DoE file and arrange it properly
study.design <- read.csv("../inst/extdata/Organs_DoE.tsv", sep='\t', header = FALSE)
study.design <- rearrange_organs_design(study.design)
The data is now cleaned up just enough to demonstrate the effect of normalization. In a proper analysis one may want to do additional checks, cleaning, and summarization steps. However, none of those affect the use of CONSTANd and one can use CONSTANd on the PSM, peptide or protein level, although the outcome may of course differ.
# the data has been summarized to peptide level
head(rownames(organs[[1]]))
# only the quantification columns were kept, which here happen to be all of them
quanCols <- colnames(organs[[1]])
quanCols
organs <- lapply(organs, function(x) x[,quanCols])
# make unique quantification column names
for (i in seq_along(organs)) {
colnames(organs[[i]]) <- paste0('BR', i, '_', colnames(organs[[i]])) }
In MS-based proteomics, anlyzing samples even with the same intrument on the same day but in a different run can cause systematic bias between the runs. Therefore, we normalize each BR run dataframe independently:
organs.norm.obj <- lapply(organs, function(x) CONSTANd(x))
organs.norm <- lapply(organs.norm.obj, function(x) x$normalized_data)
However, after CONSTANd normalization these values become comparable. We demonstrate this again using PCA plots, after merging the non-normalized data frames and merging the normalized data frames. Let's prepare the merge by doing some administration on our column names.
merge_list_of_dataframes <- function(list.dfs) {
return(Reduce(function(x, y) {
tmp <- merge(x, y, by = 0, all = FALSE)
rownames(tmp) <- tmp$Row.names
tmp$Row.names <- NULL
return(tmp)
}, list.dfs))}
organs.df <- merge_list_of_dataframes(organs)
organs.norm.df <- merge_list_of_dataframes(organs.norm)
We can now make PCA plots using data from all 4 runs:
# PCA can't deal with NA values: impute as zero (this makes sense in a multiplex)
organs.df[is.na(organs.df)] <- 0
organs.norm.df[is.na(organs.norm.df)] <- 0
# scale raw data
organs.pc <- prcomp(x=t(organs.df), scale. = TRUE)
# CONSTANd already scaled the data
organs.norm.pc <- prcomp(x=t(organs.norm.df))
The randomized study design and separate design file makes coloring the PCA plot a bit tedious.
# prepare plot color mappings
organs.pcdf <- as.data.frame(organs.pc$x)
organs.pcdf <- merge(organs.pcdf, study.design, by=0)
rownames(organs.pcdf) <- organs.pcdf$Row.names
organs.norm.pcdf <- as.data.frame(organs.norm.pc$x)
organs.norm.pcdf <- merge(organs.norm.pcdf, study.design, by=0)
rownames(organs.norm.pcdf) <- organs.norm.pcdf$Row.names
p1 <- ggplot(organs.pcdf, aes(x=PC1, y=PC2, color=condition)) + ggtitle("raw intensities") + geom_point() + theme_bw() + theme(legend.position = "none") + theme(plot.title = element_text(hjust=0.5))
p2 <- ggplot(organs.norm.pcdf, aes(x=PC1, y=PC2, color=condition)) + ggtitle("normalized intensities") + geom_point() + theme_bw() + theme(legend.position = "right") + theme(plot.title = element_text(hjust=0.5))
shared_legend <- extract_legend(p2)
grid.arrange(arrangeGrob(p1, p2 + theme(legend.position = "none"), ncol = 2),
shared_legend, widths=c(5,1))
Optional arguments
target
The variable target=1 sets the mean quantification value in each row and column during the raking process. After normalization, each row and column mean in the matrix will be equal to target up to a certain precision (see further). Setting a custom value does not affect the normalization procedure, but does scale the output values. Caution: quantification values from matrices normalized with different target values may under normal circumstances not be directly compared. Change this value only when you need to and when you know what you are doing.
maxIterations and precision
The variable maxIterations=50 is an integer value that defines a hard stop on the number of raking cycles, forcing the algorithm to return even if the desired precision was not reached (a warning will be printed). The variable precision=1e-5) defines the stopping criteria based on the L1-norm as defined by Friedrich Pukelsheim, Bruno Simeone in "On the Iterative Proportional Fitting Procedure: Structure of Accumulation Points and L1-Error Analysis".
Making MA plots {#s:maplots}
We recommend building upon the following function to make MA plots when using CONSTANd.
MAplot <- function(x,y,use.order=FALSE, R=NULL, cex=1.6, showavg=TRUE) {
# make an MA plot of y vs. x that shows the rolling average,
M <- log2(y/x)
xlab = 'A'
if (!is.null(R)) {r <- R; xlab = "A (re-scaled)"} else r <- 1
A <- (log2(y/r)+log2(x/r))/2
if (use.order) {
orig.order <- order(A)
A <- orig.order
M <- M[orig.order]
xlab = "original rank of feature magnitude within IPS"
}
# select only finite values
use <- is.finite(M)
A <- A[use]
M <- M[use]
# plot
print(var(M))
plot(A, M, xlab=xlab, cex.lab=cex, cex.axis=cex)
# rolling average
if (showavg) { lines(lowess(M~A), col='red', lwd=5) }
}
When verifying that input data adheres to the assumptions in section \@ref(s:assumptions), one could use this function (preferably on each pair of biological conditions in each IPS) without optional arguments:
# Example for cerebrum and kidney conditions in run (IPS) 1 in the Organs data set.
cerebrum1.colnames <- rownames(study.design %>% dplyr::filter(run=='BR1', condition=='cerebrum'))
kidney1.colnames <- rownames(study.design %>% dplyr::filter(run=='BR1', condition=='kidney'))
cerebrum1.raw <- organs.df[,cerebrum1.colnames]
kidney1.raw <- organs.df[,kidney1.colnames]
# rowMeans in case of multiple columns
if (!is.null(dim(cerebrum1.raw))) cerebrum1.raw <- rowMeans(cerebrum1.raw)
if (!is.null(dim(kidney1.raw))) kidney1.raw <- rowMeans(kidney1.raw)
MAplot(cerebrum1.raw, kidney1.raw)
MA plot trick for normalized data
On the normalized data, then, we recommend using the R vector (be careful to select only non-discarded features) from the CONSTANd output to return a sensible measure of feature magnitude after normalization:
cerebrum1.norm <- organs.norm.df[,cerebrum1.colnames]
kidney1.norm <- organs.norm.df[,kidney1.colnames]
# rowMeans in case of multiple columns
if (!is.null(dim(cerebrum1.norm))) cerebrum1.norm <- rowMeans(cerebrum1.norm)
if (!is.null(dim(kidney1.norm))) kidney1.norm <- rowMeans(kidney1.norm)
BR1.R <- organs.norm.obj[[1]]$R[rownames(organs.norm.df)]
MAplot(cerebrum1.norm, kidney1.norm, R=BR1.R)
Note that this re-scaling only affects the horizontal A axis values, rendering such an MA plot comparable with one representing raw data. If the S vector from the output, then, has values close to 1 (meaning the sample sizes were very similar) then this MA plot will look a lot like the one with raw data.
If, for some reason, you do not want to re-scale the normalized data using the R vector, you can also use use.order=TRUE option to make the A axis (horizontal axis) display the data in the order they appear in the MAplot function inputs. In that case, the R argument has no effect.
Subsequent Differential Expression Analysis (DEA)
You can use CONSTANd-normalized data for DEA just like you usually would (on non-transformed quantification values). For instance, you could perform a moderated t-test:
get_design_matrix <- function(referenceCondition, study.design) {
# ANOVA-like design matrix for use in moderated_ttest, indicating group
# (condition) membership of each entry in all_samples.
otherConditions = setdiff(unique(study.design$condition), referenceCondition)
all_samples = rownames(study.design)
N_samples = length(all_samples)
N_conditions = 1+length(otherConditions)
design = matrix(rep(0,N_samples*N_conditions), c(N_samples, N_conditions))
design[, 1] <- 1 # reference gets 1 everywhere
rownames(design) <- all_samples
colnames(design) <- c(referenceCondition, otherConditions)
for (i in 1:N_samples) { # for each channel in each condition, put a "1" in the design matrix.
design[as.character(rownames(study.design)[i]), as.character(study.design$condition[i])] <- 1 }
return(design)
}
moderated_ttest <- function(dat, design_matrix, scale) {
# inspired by http://www.biostat.jhsph.edu/~kkammers/software/eupa/R_guide.html
design_matrix <- design_matrix[match(colnames(dat), rownames(design_matrix)),] # fix column order
reference_condition <- colnames(design_matrix)[colSums(design_matrix) == nrow(design_matrix)]
nfeatures <- dim(dat)[1]
fit <- limma::eBayes(lmFit(dat, design_matrix))
reference_averages <- fit$coefficients[,reference_condition]
logFC <- fit$coefficients
logFC <- log2((logFC+reference_averages)/reference_averages)
# correct the reference
logFC[,reference_condition] <- logFC[,reference_condition] - 1
# moderated q-value corresponding to the moderated t-statistic
if(nfeatures>1) {
q.mod <- apply(X = fit$p.value, MARGIN = 2, FUN = p.adjust, method='BH')
} else q.mod <- fit$p.value
return(list(logFC=logFC, qval=q.mod))
}
design_matrix <- get_design_matrix('cerebrum', study.design)
result <- moderated_ttest(organs.norm.df, design_matrix)
Session info
sessionInfo()
PDiracDelta/CONSTANd documentation built on March 9, 2021, 11:47 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE,dev="CairoPNG") library(knitr) library(tidyr) library(ggplot2) library(gridExtra) library(magick) library(limma)
Introduction
CONSTANd normalization has been proven to work on mass spectrometry intensity values (Maes et al.) in shotgun proteomics and on RNA-seq count data (Van Houtven et al.), but is in principle applicable to any kind of quantification matrices, for instance in other -omics disciplines. It normalizes the data matrix data by raking (using the RAS method by Bacharach, see references) the Nrows by Ncols matrix such that the row means and column means equal 1. The result is a normalized data matrix K=RAS, a product of row mulipliers R and column multipliers S with the original matrix A. Missing information needs to be presented as NA values and not as zero values, because CONSTANd is able to ignore missing values when calculating the mean.
Getting started
The CONSTANd package has to be loaded to use the function.
To load CONSTANd, type the following command after launching R.
#install.packages("BiocManager") #BiocManager::install("CONSTANd") library(CONSTANd)
Data assumptions {#s:assumptions}
To warrant that the normalization procedure itself is not biased, three assumptions about the data are made (as for any type of data-driven normalization method) which may be verified by making MA-plots (see examples in section \@ref(s:maplots)):
- The majority of features (peptides/genes/...) are not differentially expressed, to avoid a bias in the estimate of the mean value. It is assumed that up-down shifts are not due to biological causes. The reference set used in the normalization step is the set of all peptides identified in the experiment.
MA-plot: the observations form a single 'cloud' with a dense center and less dense edges, as opposed to, for instance, two clouds or a cloud with uniform density. - The number of up-regulated features is roughly equal to the number of down-regulated ones. If the data were skewed, this would lead to a bias in the normalization result.
MA-plot: the cloud of observations exhibits a bilateral symmetry about some axis (usually horizontal, but inclinations may occur). - Any systematic bias present is linearly proportional to the magnitude of the quantification values. Only this way it is possible to find one appropriate normalization factor for each quantification sample.
MA-plot: the axis of bilateral symmetry is a straight line (which may be inclined), i.e., the moving average M-values of the central cloud form an approximately straight line.
Example of a bad MA plot
The following MA plot violates all three assumptions. In case your MA plots exhibit any of such characteristics, it is not advised to use CONSTANd - or any other data-driven method - to normalize your data, and to instead use advanced model-based approaches or, if appropriate, clean your data first.
par(mfrow=c(1,1)) include_graphics('bad_MA.png')
Normalizing a data matrix (single assay)
The data matrix is an Nrows by Ncols matrix (e.g. an assay from a SummarizedExperiment), representing quantification values of any kind, as long as they form an identically processed (sub)set (IPS). An IPS consists of measurements from biological samples who have been processed from biological harvesting of the tissue to the detection and quantitation of the features of interest (RNA, proteins, ...) in near-identical fashion. In practice, this means that there exist no factors (except those of biological interest) which may cause a systematic bias in the quantification values. One _counter_example in proteomics is samples measured on a different day or a different instrument, one _counter_example in RNA-seq would be samples measured after using different library preparation protocols. In the next section, we demonstrate how in such such cases the IPSs are normalized separately as if they were individual assays (see this section), after which they may be re-combined.
First, let's load the 'Leishmania' mRNA-level data set from GEO (Gene Expression Omnibus) with identifier GSE95353. You can either download and unzip it yourself, are use the file that comes with this package.
leish <- read.csv("../inst/extdata/GSE95353_SL_SEQ_LOGSTAT_counttable.txt" , sep='\t')
We only need the LOG and STAT samples, which correspond to biological quadruplicates of Leishmania in two lifestages (LOG and STAT).
rownames(leish) <- leish$gene leish <- leish[, c("SL.LOG1", "SL.LOG2", "SL.LOG3", "SL.LOG4", "SL.STAT1", "SL.STAT2", "SL.STAT3", "SL.STAT4")] conditions <- c(rep("LOG", 4), rep("STAT", 4))
Let's remove genes with low count numbers (median > 0) for good practice, and make sure that zero counts are replaced by NA to use CONSTANd.
# first set NA to 0 for calculating the row medians leish[is.na(leish)] <- 0 leish <- leish[apply(leish, 1, median, na.rm = FALSE) > 0,] leish.norm <- CONSTANd(leish)$normalized_data
Let's make a PCA plot (optionally impute NA as zero) before and after normalization to see if it was successful:
leish[is.na(leish)] <- 0 # scale raw data leish.pc.raw <- prcomp(x=t(leish), scale. = TRUE) leish.norm[is.na(leish.norm)] <- 0 # CONSTANd already scaled the data leish.pc.norm <- prcomp(x=t(leish.norm)) par(mfrow=c(1,2)) colors = ifelse(conditions == 'LOG', 'blue', 'red') plot(leish.pc.raw$x[,'PC1'], leish.pc.raw$x[,'PC2'], xlab='PC1', ylab='PC2', main='raw counts', col=colors) plot(leish.pc.norm$x[,'PC1'], leish.pc.norm$x[,'PC2'], xlab='PC1', ylab='PC2', main='normalized counts', col=colors)
Indeed, the PCA plot indicates that the first principal component, which is the direction along which there is most varability in the measurements, successfully defines a separation between the biological conditions. In other words: after normalization the biological differences between the samples dominate any other sources of systematic bias.
Normalization across multiple assays
Here, we demonstrate multiple IPSs involved in one experiment may be normalized separately as if they were individual assays (see previous section), after which they may be re-combined for joint analysis. We demonstrate CONSTANd normalization across multiple assays using the Organs example data set from QCQuan, a proteomics webtool built around the CONSTANd algorithm. The Organs data set entails samples from 8 different mouse organs in biological quadruplicate, whose peptides have been measured in 4 different instrument runs (1 per mouse) using TMT-labeled (8-plex) tandem mass spectrometry.
par(mfrow=c(1,1)) include_graphics('organ_experiment_design.png')
You can either download the original data from the web, unzip it and set LOAD_ORIGINAL=TRUE, or load the cleaned data from this package.
LOAD_ORIGINAL = FALSE # load extra functions to clean the original data and study design files source('../inst/script/functions.R') if (LOAD_ORIGINAL) { # load original from web BR1 <- read.csv('Organs_PSMs_1.txt', sep='\t') BR2 <- read.csv('Organs_PSMs_2.txt', sep='\t') BR3 <- read.csv('Organs_PSMs_3.txt', sep='\t') BR4 <- read.csv('Organs_PSMs_4.txt', sep='\t') organs <- list(BR1, BR2, BR3, BR4) organs <- clean_organs(organs) } else { # load cleaned file that comes with this package organs <- readRDS('../inst/extdata/organs_cleaned.RDS') } # construct the study design from the QCQuan DoE file and arrange it properly study.design <- read.csv("../inst/extdata/Organs_DoE.tsv", sep='\t', header = FALSE) study.design <- rearrange_organs_design(study.design)
The data is now cleaned up just enough to demonstrate the effect of normalization. In a proper analysis one may want to do additional checks, cleaning, and summarization steps. However, none of those affect the use of CONSTANd and one can use CONSTANd on the PSM, peptide or protein level, although the outcome may of course differ.
# the data has been summarized to peptide level head(rownames(organs[[1]])) # only the quantification columns were kept, which here happen to be all of them quanCols <- colnames(organs[[1]]) quanCols organs <- lapply(organs, function(x) x[,quanCols]) # make unique quantification column names for (i in seq_along(organs)) { colnames(organs[[i]]) <- paste0('BR', i, '_', colnames(organs[[i]])) }
In MS-based proteomics, anlyzing samples even with the same intrument on the same day but in a different run can cause systematic bias between the runs. Therefore, we normalize each BR run dataframe independently:
organs.norm.obj <- lapply(organs, function(x) CONSTANd(x)) organs.norm <- lapply(organs.norm.obj, function(x) x$normalized_data)
However, after CONSTANd normalization these values become comparable. We demonstrate this again using PCA plots, after merging the non-normalized data frames and merging the normalized data frames. Let's prepare the merge by doing some administration on our column names.
merge_list_of_dataframes <- function(list.dfs) { return(Reduce(function(x, y) { tmp <- merge(x, y, by = 0, all = FALSE) rownames(tmp) <- tmp$Row.names tmp$Row.names <- NULL return(tmp) }, list.dfs))} organs.df <- merge_list_of_dataframes(organs) organs.norm.df <- merge_list_of_dataframes(organs.norm)
We can now make PCA plots using data from all 4 runs:
# PCA can't deal with NA values: impute as zero (this makes sense in a multiplex) organs.df[is.na(organs.df)] <- 0 organs.norm.df[is.na(organs.norm.df)] <- 0 # scale raw data organs.pc <- prcomp(x=t(organs.df), scale. = TRUE) # CONSTANd already scaled the data organs.norm.pc <- prcomp(x=t(organs.norm.df))
The randomized study design and separate design file makes coloring the PCA plot a bit tedious.
# prepare plot color mappings organs.pcdf <- as.data.frame(organs.pc$x) organs.pcdf <- merge(organs.pcdf, study.design, by=0) rownames(organs.pcdf) <- organs.pcdf$Row.names organs.norm.pcdf <- as.data.frame(organs.norm.pc$x) organs.norm.pcdf <- merge(organs.norm.pcdf, study.design, by=0) rownames(organs.norm.pcdf) <- organs.norm.pcdf$Row.names
p1 <- ggplot(organs.pcdf, aes(x=PC1, y=PC2, color=condition)) + ggtitle("raw intensities") + geom_point() + theme_bw() + theme(legend.position = "none") + theme(plot.title = element_text(hjust=0.5)) p2 <- ggplot(organs.norm.pcdf, aes(x=PC1, y=PC2, color=condition)) + ggtitle("normalized intensities") + geom_point() + theme_bw() + theme(legend.position = "right") + theme(plot.title = element_text(hjust=0.5)) shared_legend <- extract_legend(p2) grid.arrange(arrangeGrob(p1, p2 + theme(legend.position = "none"), ncol = 2), shared_legend, widths=c(5,1))
Optional arguments
target
The variable target=1 sets the mean quantification value in each row and column during the raking process. After normalization, each row and column mean in the matrix will be equal to target up to a certain precision (see further). Setting a custom value does not affect the normalization procedure, but does scale the output values. Caution: quantification values from matrices normalized with different target values may under normal circumstances not be directly compared. Change this value only when you need to and when you know what you are doing.
maxIterations and precision
The variable maxIterations=50 is an integer value that defines a hard stop on the number of raking cycles, forcing the algorithm to return even if the desired precision was not reached (a warning will be printed). The variable precision=1e-5) defines the stopping criteria based on the L1-norm as defined by Friedrich Pukelsheim, Bruno Simeone in "On the Iterative Proportional Fitting Procedure: Structure of Accumulation Points and L1-Error Analysis".
Making MA plots {#s:maplots}
We recommend building upon the following function to make MA plots when using CONSTANd.
MAplot <- function(x,y,use.order=FALSE, R=NULL, cex=1.6, showavg=TRUE) { # make an MA plot of y vs. x that shows the rolling average, M <- log2(y/x) xlab = 'A' if (!is.null(R)) {r <- R; xlab = "A (re-scaled)"} else r <- 1 A <- (log2(y/r)+log2(x/r))/2 if (use.order) { orig.order <- order(A) A <- orig.order M <- M[orig.order] xlab = "original rank of feature magnitude within IPS" } # select only finite values use <- is.finite(M) A <- A[use] M <- M[use] # plot print(var(M)) plot(A, M, xlab=xlab, cex.lab=cex, cex.axis=cex) # rolling average if (showavg) { lines(lowess(M~A), col='red', lwd=5) } }
When verifying that input data adheres to the assumptions in section \@ref(s:assumptions), one could use this function (preferably on each pair of biological conditions in each IPS) without optional arguments:
# Example for cerebrum and kidney conditions in run (IPS) 1 in the Organs data set. cerebrum1.colnames <- rownames(study.design %>% dplyr::filter(run=='BR1', condition=='cerebrum')) kidney1.colnames <- rownames(study.design %>% dplyr::filter(run=='BR1', condition=='kidney')) cerebrum1.raw <- organs.df[,cerebrum1.colnames] kidney1.raw <- organs.df[,kidney1.colnames] # rowMeans in case of multiple columns if (!is.null(dim(cerebrum1.raw))) cerebrum1.raw <- rowMeans(cerebrum1.raw) if (!is.null(dim(kidney1.raw))) kidney1.raw <- rowMeans(kidney1.raw) MAplot(cerebrum1.raw, kidney1.raw)
MA plot trick for normalized data
On the normalized data, then, we recommend using the R vector (be careful to select only non-discarded features) from the CONSTANd output to return a sensible measure of feature magnitude after normalization:
cerebrum1.norm <- organs.norm.df[,cerebrum1.colnames] kidney1.norm <- organs.norm.df[,kidney1.colnames] # rowMeans in case of multiple columns if (!is.null(dim(cerebrum1.norm))) cerebrum1.norm <- rowMeans(cerebrum1.norm) if (!is.null(dim(kidney1.norm))) kidney1.norm <- rowMeans(kidney1.norm) BR1.R <- organs.norm.obj[[1]]$R[rownames(organs.norm.df)] MAplot(cerebrum1.norm, kidney1.norm, R=BR1.R)
Note that this re-scaling only affects the horizontal A axis values, rendering such an MA plot comparable with one representing raw data. If the S vector from the output, then, has values close to 1 (meaning the sample sizes were very similar) then this MA plot will look a lot like the one with raw data.
If, for some reason, you do not want to re-scale the normalized data using the R vector, you can also use use.order=TRUE option to make the A axis (horizontal axis) display the data in the order they appear in the MAplot function inputs. In that case, the R argument has no effect.
Subsequent Differential Expression Analysis (DEA)
You can use CONSTANd-normalized data for DEA just like you usually would (on non-transformed quantification values). For instance, you could perform a moderated t-test:
get_design_matrix <- function(referenceCondition, study.design) { # ANOVA-like design matrix for use in moderated_ttest, indicating group # (condition) membership of each entry in all_samples. otherConditions = setdiff(unique(study.design$condition), referenceCondition) all_samples = rownames(study.design) N_samples = length(all_samples) N_conditions = 1+length(otherConditions) design = matrix(rep(0,N_samples*N_conditions), c(N_samples, N_conditions)) design[, 1] <- 1 # reference gets 1 everywhere rownames(design) <- all_samples colnames(design) <- c(referenceCondition, otherConditions) for (i in 1:N_samples) { # for each channel in each condition, put a "1" in the design matrix. design[as.character(rownames(study.design)[i]), as.character(study.design$condition[i])] <- 1 } return(design) } moderated_ttest <- function(dat, design_matrix, scale) { # inspired by http://www.biostat.jhsph.edu/~kkammers/software/eupa/R_guide.html design_matrix <- design_matrix[match(colnames(dat), rownames(design_matrix)),] # fix column order reference_condition <- colnames(design_matrix)[colSums(design_matrix) == nrow(design_matrix)] nfeatures <- dim(dat)[1] fit <- limma::eBayes(lmFit(dat, design_matrix)) reference_averages <- fit$coefficients[,reference_condition] logFC <- fit$coefficients logFC <- log2((logFC+reference_averages)/reference_averages) # correct the reference logFC[,reference_condition] <- logFC[,reference_condition] - 1 # moderated q-value corresponding to the moderated t-statistic if(nfeatures>1) { q.mod <- apply(X = fit$p.value, MARGIN = 2, FUN = p.adjust, method='BH') } else q.mod <- fit$p.value return(list(logFC=logFC, qval=q.mod)) } design_matrix <- get_design_matrix('cerebrum', study.design) result <- moderated_ttest(organs.norm.df, design_matrix)
Session info
sessionInfo()
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.