Skip to content

Commit

Permalink
The association tests for correlation between different traits now pr…
Browse files Browse the repository at this point in the history 
…oduce p-values irrespective of whether principal components were computed or not.

git-svn-id: https://hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/RnBeads@114889 bc3139a8-67e5-0310-9ffc-ced21a209358
  • Loading branch information
y.assenov committed Mar 16, 2016
1 parent e01ef7e commit 3486836
Show file tree
Hide file tree
Showing 4 changed files with 42 additions and 41 deletions.
48 changes: 24 additions & 24 deletions R/batch.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -184,7 +184,7 @@ test.traits <- function(x, y, perm.matrix = NULL) {
"test" = as.character(NA),
"correlation" = as.double(NA),
"pvalue" = as.double(NA))

## Focus on common values
if (class(x) == "Date") { x <- as.integer(x) }
if (class(y) == "Date") { y <- as.integer(y) }
Expand All @@ -211,7 +211,7 @@ test.traits <- function(x, y, perm.matrix = NULL) {
return(result)
}
}

## Perform a test or compute correlation
get.p <- function(expr) { tryCatch(suppressWarnings(expr$p.value), error = function(er) { as.double(NA) }) }
if (is.factor(x)) {
Expand Down Expand Up @@ -248,7 +248,7 @@ test.traits <- function(x, y, perm.matrix = NULL) {
result[["pvalue"]] <- mean(values[1] <= values)
}
}

if (is.na(result[["pvalue"]])) {
result[["error"]] <- "test failed"
}
Expand Down Expand Up @@ -335,7 +335,7 @@ plot.heatmap.pc.pvalues <- function(report, tbl, fname, width = NULL, height = N
geom_text(size = 3) + scale_x_discrete(expand = c(0, 0)) + scale_y_discrete(expand = c(0, 0)) +
theme(panel.border = element_blank(), panel.grid = element_blank()) +
theme(plot.margin = unit(0.1 + c(0, 0, 0, 0), "in")) +
theme(axis.ticks = element_blank(), legend.position = "none")
theme(axis.ticks = element_blank(), legend.position = "none")
## Fix the areas for x and y axis labels
pp <- suppressWarnings(ggplot_gtable(ggplot_build(pp)))
if (is.null(ylab)) {
Expand Down Expand Up @@ -457,7 +457,7 @@ plot.heatmap.symm <- function(report, tbl.symm, tbl.failures = NULL, fname) {
#' @param target \code{character} singleton specifying the level of DNA methylation infromation. If this is
#' \code{"sites"}, the DNA methylation information for the individual sites or probes is analyzed.
#' Otherwise, this should be one of the supported region types, as returned by
#' \code{\link{rnb.region.types}}.
#' \code{\link{rnb.region.types}}.
#' @return Results of the dimension reduction in the form of a list with the following elements:
#' \describe{
#' \item{\code{pca}}{Results of the PCA as returned by the function \code{\link{prcomp}}.}
Expand Down Expand Up @@ -596,7 +596,7 @@ rnb.section.dreduction <- function(report, pcoordinates, sample.phenotypes = NUL
stop("invalid value for report")
}
nsamples <- validate.pcoordinates.all(pcoordinates)

if (is.character(nsamples)) {
stop(nsamples)
}
Expand Down Expand Up @@ -626,7 +626,7 @@ rnb.section.dreduction.internal <- function(report, pcoordinates, sample.phenoty
col.visual <- names(rnb.sample.groups(sample.phenotypes, rnb.getOption("exploratory.columns")))
use.colors <- (length(col.visual) != 0)
logger.info(c("Mapping", length(col.visual), "traits to point colors and types"))

setting.names <- list("Location type" = targets,
"Principal components" = c("1" = "first and second", "2" = "second and third"),
"Distance" = c("manhattan" = "Manhattan", "euclidean" = "Euclidean"))
Expand All @@ -642,7 +642,7 @@ rnb.section.dreduction.internal <- function(report, pcoordinates, sample.phenoty
setting.names[["Sample color"]] <- snames
}
rm(snames)

create.scatters <- function(dpoints, fprefix) {
report.plots <- list()
dframe <- data.frame(x = dpoints[, 1], y = dpoints[, 2], id = rownames(dpoints))
Expand Down Expand Up @@ -881,7 +881,7 @@ rnb.section.dreduction.internal <- function(report, pcoordinates, sample.phenoty
"components ", txt[2], " available in ", txt[3], "comma-separated values ", txt[4], " accompanying this ",
"report.")
rnb.add.paragraph(report, txt)

colnames(var.tbl) <- c("Location Type", "Number of Components", "Full Table File")
rnb.add.table(report, var.tbl, row.names = FALSE)
}
Expand Down Expand Up @@ -1001,7 +1001,7 @@ rnb.step.dreduction <- function(rnb.set, report, return.coordinates = FALSE) {
return(list(report = report, pcoordinates = pcoordinates))
}
return(report)

}

########################################################################################################################
Expand Down Expand Up @@ -1060,7 +1060,7 @@ rnb.execute.batcheffects <- function(rnb.set, pcoordinates = NULL) {
logger.start(fname = NA) # initialize console logger
}
logger.start("Tests for Associations")

## Get a list of comparable traits
predefined.columns <- rnb.getOption("exploratory.columns")
if (!is.null(predefined.columns)) {
Expand Down Expand Up @@ -1098,13 +1098,13 @@ rnb.execute.batcheffects <- function(rnb.set, pcoordinates = NULL) {
return(NULL)
}
logger.info(c("Testing the following traits for associations:", paste(names(traits), collapse = "; ")))

result <- list()

## Create sample permutations if necessary
perm.matrix <- NULL
perm.count <- rnb.getOption("exploratory.correlation.permutations")
if ((!is.null(pcoordinates)) && perm.count != 0 && sum(!sapply(traits, is.factor)) >= 2) {
if (perm.count != 0 && ((!is.null(pcoordinates)) || sum(!sapply(traits, is.factor)) >= 2)) {
perm.matrix <- mapply(sample, rep(nrow(pheno.table), times = perm.count))
perm.matrix[, 1] <- 1:nrow(perm.matrix)
result[["permutations"]] <- perm.matrix
Expand All @@ -1121,7 +1121,7 @@ rnb.execute.batcheffects <- function(rnb.set, pcoordinates = NULL) {
if (ncol(dpoints) > pc.association.count) {
dpoints <- dpoints[, 1:pc.association.count]
}

init.matrix <- function(x) {
matrix(x, nrow = NT, ncol = ncol(dpoints),
dimnames = list(names(traits), "Principal component" = 1:ncol(dpoints)))
Expand All @@ -1130,7 +1130,7 @@ rnb.execute.batcheffects <- function(rnb.set, pcoordinates = NULL) {
table.test.names <- init.matrix(as.character(NA))
table.correlations <- init.matrix(as.double(NA))
table.pvalues <- init.matrix(as.double(NA))

for (i in 1:NT) {
for (j in 1:ncol(dpoints)) {
t.result <- test.traits(traits[[i]], dpoints[, j], perm.matrix)
Expand All @@ -1151,7 +1151,7 @@ rnb.execute.batcheffects <- function(rnb.set, pcoordinates = NULL) {
rm(i, j, t.result)
logger.status("Computed correlations between principal components and traits.")
}

## Compute correlations and perform tests
if (NT > 1) {
init.matrix <- function(x) {
Expand Down Expand Up @@ -1233,7 +1233,7 @@ rnb.section.batcheffects <- function(report, batcheffects) {
if (rnb.getOption("logging") && logger.isinitialized() == FALSE) {
logger.start(fname = NA) # initialize console logger
}

txt <- c("In this section, different properties of the dataset are tested for significant associations. The ",
"properties can include sample coordinates in the principal component space, phenotype traits and ",
"intensities of control probes. The tests used to calculate a p-value given two properties depend on the ",
Expand Down Expand Up @@ -1262,7 +1262,7 @@ rnb.section.batcheffects <- function(report, batcheffects) {
## Summarize the results of association between principal components and traits

if ("pc" %in% names(batcheffects) && length(batcheffects[["pc"]])>0) {

targets <- names(batcheffects[["pc"]])
names(targets) <- 1:length(targets)
setting.names <- list("Region type" = targets)
Expand Down Expand Up @@ -1374,7 +1374,7 @@ rnb.section.batcheffects <- function(report, batcheffects) {
## Create a triangular heatmap of the calculated p-values for associations
tbl <- trait.tables$pvalues
plots.associations[[2]] <- plot.heatmap.symm(report, tbl, NULL, "heatmap_traits_pvalues")

## Attach the table of the calculated p-values for associations
tbl <- cbind("Traits" = rownames(tbl), as.data.frame(tbl, check.names = FALSE))
fname <- "pvalues_traits.csv"
Expand All @@ -1392,10 +1392,10 @@ rnb.section.batcheffects <- function(report, batcheffects) {
txt <- c(txt, "(2) Table of resulting p-values from the performed tests on pairs of traits. Significant ",
"p-values (less than ", rnb.getOption("exploratory.correlation.pvalue.threshold"), ") are printed in pink boxes ",
"Non-significant values are represented by blue boxes. White cells, if present, denote missing values.")
setting.names <- list("Heatmap of" = c("tests" = "tests performed", "pvalues" = "p-values"))
setting.names <- list("Heatmap of" = c("tests" = "tests performed", "pvalues" = "p-values"))
report <- rnb.add.figure(report, txt, plots.associations, setting.names, selected.image = 2)
rm(plots.associations, txt, tbl, tbl.failures, fname, setting.names)

if (!all(is.na(trait.tables$correlations))) {
## Attach the table of correlations to the report
txt <- "In some cases, a correlation was computed between a pair of traits."
Expand All @@ -1412,7 +1412,7 @@ rnb.section.batcheffects <- function(report, batcheffects) {
rnb.add.paragraph(report, txt)
}
}

return(report)
}

Expand Down Expand Up @@ -1458,7 +1458,7 @@ rnb.step.batcheffects <- function(rnb.set, report, pcoordinates, return.permutat
if (rnb.getOption("logging") && logger.isinitialized() == FALSE) {
logger.start(fname = NA) # initialize console logger
}

if (is.null(pheno(rnb.set))) {
txt <- "Batch effects were not studied because the dataset contains no phenotype data."
report <- rnb.add.section(report, "Batch Effects", txt)
Expand Down
28 changes: 14 additions & 14 deletions R/batch.quality.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,17 +71,17 @@ rnb.execute.batch.qc <- function(rnb.set, pcoordinates, permutations = NULL) {
targets <- names(pcoordinates)
names(targets) <- 1:length(targets)
results<-lapply(1:length(targets), function(target.id){

dpoints <- pcoordinates[[target.id]]$pca$x
if (ncol(dpoints) > pc.association.count) {
dpoints <- dpoints[, 1:pc.association.count]
}

result <- list("correlations" = list())
if (!is.null(permutations)) {
result[["pvalues"]] <- list()
}

get.probes.channel <- function(channel, ids) {
i <- intersect(ids, rownames(channel))
if (length(i) == 0) {
Expand All @@ -98,7 +98,7 @@ rnb.execute.batch.qc <- function(rnb.set, pcoordinates, permutations = NULL) {
matrix(as.double(NA), nrow = N, ncol = ncol(dpoints),
dimnames = list("Probe" = 1:N, "Principal component" = 1:ncol(dpoints)))
}

## Compute correlations between QC probe values and principal components
target2ids <- tapply(control.probe.infos[[id.col]], control.probe.infos[[type.col]], as.character)
for (ci in 1:length(CONTROL.TYPES)) {
Expand All @@ -119,7 +119,7 @@ rnb.execute.batch.qc <- function(rnb.set, pcoordinates, permutations = NULL) {
}
table.correlations <- init.matrix(nrow(matrix.channel))
table.pvalues <- init.matrix(nrow(matrix.channel))

for (i in 1:nrow(matrix.channel)) {
for (j in 1:ncol(dpoints)) {
t.result <- test.traits(matrix.channel[i, ], dpoints[, j], permutations)
Expand Down Expand Up @@ -161,17 +161,17 @@ rnb.section.batch.qc <- function(report, qccorrelations) {
stop("invalid value for report")
}
## TODO: Validate qccorrelations

targets<-sapply(qccorrelations, attr, "Target")
report.descriptions <- c("correlations" = as.character(NA), "pvalues" = as.character(NA))

channels<-NULL
ctypes<-NULL
pvals.present<-NULL
table.pvalue.links<-list()

all.plots<-lapply(1:length(targets), function(target.id){

## Create heatmaps of correlations and p-values
pvals.present <<- ("pvalues" %in% names(qccorrelations[[target.id]]))
ctypes <<- names(qccorrelations[[target.id]][["correlations"]])
Expand All @@ -185,7 +185,7 @@ rnb.section.batch.qc <- function(report, qccorrelations) {
table.pvalue.links[[target.id]] <<- matrix(as.character(NA), nrow = length(ctypes), ncol = length(channels),
dimnames = list(ctypes, channels))
}

heatmap.functions <- list("correlations" = plot.heatmap.pc.correlations, "pvalues" = plot.heatmap.pc.pvalues)
for (tblname in names(qccorrelations[[target.id]])) {
heatmap.function <- heatmap.functions[[tblname]]
Expand All @@ -203,7 +203,7 @@ rnb.section.batch.qc <- function(report, qccorrelations) {
report.descriptions[tblname] <<- rplot$description
}
rplot <- rplot$plot

## Save the table of p-values to a file
if (tblname == "pvalues") {
fname <- paste("pvalues_pc_", target.id, "_", cchannel, "_", i, ".csv", sep = "")
Expand All @@ -213,7 +213,7 @@ rnb.section.batch.qc <- function(report, qccorrelations) {
utils::write.csv(tbl, file = fname.full, row.names = FALSE)
table.pvalue.links[[target.id]][i,grep(cchannel, colnames(table.pvalue.links[[target.id]]))] <<- paste("<a href=\"", rnb.get.directory(report, "data"), "/",
fname, "\">csv</a>", sep = "")

rm(fname.full, tbl)
}
} else {
Expand All @@ -230,7 +230,7 @@ rnb.section.batch.qc <- function(report, qccorrelations) {
names(setting.names[["Location type"]]) <- 1:length(targets)
names(setting.names[["Channel"]]) <- channels
names(setting.names[["Probe group"]]) <- 1:length(ctypes)

txt <- c("The heatmaps below visualize the Pearson correlation coefficients between the principal ",
"components and the signal levels of selected quality control probes.")
rnb.add.paragraph(report, txt)
Expand Down Expand Up @@ -275,7 +275,7 @@ rnb.section.batch.qc <- function(report, qccorrelations) {
#' @param rnb.set HumanMethylation450K dataset as an object of type \code{\linkS4class{RnBeadSet}}.
#' @param report Report to contain the quality-associated section. This must be an object of type
#' \code{\linkS4class{Report}}.
#' @param pcoordinates Coordinates of the samples of \code{rnbSet} in the principal components space, as returned by
#' @param pcoordinates Coordinates of the samples of \code{rnb.set} in the principal components space, as returned by
#' \code{\link{rnb.execute.dreduction}}.
#' @param permutations Matrix of sample index permutations, as returned by \code{\link{rnb.execute.batcheffects}}.
#' @return The modified report.
Expand Down
3 changes: 2 additions & 1 deletion R/main.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1346,7 +1346,8 @@ rnb.run.exploratory <- function(rnb.set, dir.reports,
report <- result$report

## Quality-associated batch effects
if (rnb.getOption("exploratory.correlation.qc") && inherits(rnb.set, "RnBeadSet")) {
if (rnb.getOption("exploratory.principal.components") != 0 && rnb.getOption("exploratory.correlation.qc") &&
inherits(rnb.set, "RnBeadSet")) {
report <- rnb.step.batch.qc(rnb.set, report, pcoordinates, permutations = result$permutations)
}
}
Expand Down
4 changes: 2 additions & 2 deletions man/run-RnBClusterRun-methods.Rd
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

0 comments on commit 3486836

Please sign in to comment.