fixes predict.design.matrix
fixes predict.design.matrix

perishky/meffil
#' Infinium HumanMethylation450 BeadChip normalization design matrix
#'
#' Design matrix derived by applying principal components analysis to control probes.
#'
#' @param qc.objects A list of outputs from \code{\link{meffil.create.qc.object}()}.
#' @param number.pcs Number of principal components to include in the design matrix (Default: all).
#' @param fixed.effects Names of columns in samplesheet that should be included as fixed effects
#' along with control matrix principal components (Default: NULL).
#' @param random.effects Names of columns in samplesheet that should be included as random effects
#' (Default: NULL).
#' @return Design matrix with one column for each of the first \code{number.pcs} prinicipal
#' components.
#'
#' @export
meffil.design.matrix <- function(qc.objects, number.pcs, fixed.effects = NULL, random.effects = NULL){
    if (!is.null(fixed.effects) && is.character(fixed.effects))
        fixed.effects <- extract.from.samplesheet(qc.objects, fixed.effects)
    if (!is.null(random.effects) && is.character(random.effects))
        random.effects <- extract.from.samplesheet(qc.objects, random.effects)
    
    pca.ret <- meffil.pcs(qc.objects, fixed.effects, random.effects)
    pca.to.design.matrix(pca.ret, number.pcs, fixed.effects, random.effects)
}

#' Calculate control probe PCs
#'
#' @param qc.objects A list of outputs from \code{\link{meffil.create.qc.object}()}.
#' @param fixed.effects Names of columns in samplesheet that should be included as fixed effects
#' along with control matrix principal components (Default: NULL).
#' @param random.effects Names of columns in samplesheet that should be included as random effects
#' (Default: NULL).
#' @return PCA of control probes
#' 
#' @export
meffil.pcs <- function(qc.objects, fixed.effects=NULL, random.effects=NULL) {
    control.matrix <- meffil.control.matrix(qc.objects, normalize=T,
                                            fixed.effects=fixed.effects,
                                            random.effects=random.effects)
    control.pca <- prcomp(control.matrix)
    return(control.pca)
}



#' Infinium HumanMethylation450 BeadChip control matrix
#'
#' Matrix containing control probe intensities from the Infinium HumanMethylation450 BeadChip.
#'
#' @param qc.objects A list of outputs from \code{\link{meffil.create.qc.object}()}.
#' @param normalize If \code{TRUE}, then control matrix is scaled and specified
#' fixed and random effects removed from the matrix.  Otherwise, the raw control matrix is returned.
#' (Default: \code{FALSE}).
#' @param fixed.effects Names of columns in samplesheet that should be included as fixed effects
#' along with control matrix principal components (Default: NULL).
#' @param random.effects Names of columns in samplesheet that should be included as random effects
#' (Default: NULL).
#' @return Matrix with one row per object consisting of control probe intensities and summaries.
#'
#' @export
meffil.control.matrix <- function(qc.objects, normalize=F,
                                  fixed.effects = NULL, random.effects = NULL) {
    if (!is.null(fixed.effects) && is.character(fixed.effects))
        fixed.effects <- extract.from.samplesheet(qc.objects, fixed.effects)
    if (!is.null(random.effects) && is.character(random.effects))
        random.effects <- extract.from.samplesheet(qc.objects, random.effects)
    
    stopifnot(all(sapply(qc.objects, meffil:::is.qc.object)))
    stopifnot(length(qc.objects) >= 2)
    stopifnot(is.null(fixed.effects) || nrow(fixed.effects) == length(qc.objects))
    stopifnot(is.null(random.effects) || nrow(random.effects) == length(qc.objects))

    control.matrix <- t(sapply(qc.objects, function(object) object$controls))

    if (normalize) {
        control.matrix <- meffil:::impute.matrix(control.matrix, margin = 2)
        control.matrix <- scale(control.matrix)
        control.matrix[control.matrix > 3] <- 3
        control.matrix[control.matrix < -3] <- -3
        if (!is.null(fixed.effects) || !is.null(random.effects)) {
            control.matrix <- meffil:::adjust.columns(control.matrix, fixed.effects, random.effects)
            control.matrix <- scale(control.matrix)
        }
    }
    control.matrix
}

#' Given a partition of the QC objects into two sets,
#' use the design matrix from one partition to predict the design
#' matrix for the other partition.
#' @param qc.objects A list of outputs from \code{\link{meffil.create.qc.object}()}.
#' @param number.pcs Number of principal components to include in the design matrix (Default: all).
#' @param new.idx Indices into \code{qc.objects} that identify the objects for which
#' the design matrix will be predicted.  The initial design matrix will be derived from
#' the other objects.
#' @param fixed.effects Names of columns in samplesheet that should be included as fixed effects
#' along with control matrix principal components (Default: NULL).
#' @param random.effects Names of columns in samplesheet that should be included as random effects
#' (Default: NULL).
#' @return Design matrix with one column for the \code{new.idx} partition predicted from the
#' first \code{number.pcs} prinicipal components of the objects. 
predict.design.matrix <- function (qc.objects, number.pcs, new.idx,
                                   fixed.effects = NULL,
                                   random.effects = NULL) {
    stopifnot(all(sapply(qc.objects, meffil:::is.qc.object)))
    stopifnot(all(new.idx %in% 1:length(qc.objects)))
    
    if (!is.null(fixed.effects) && is.character(fixed.effects))
        fixed.effects <- meffil:::extract.from.samplesheet(qc.objects, fixed.effects)
    if (!is.null(random.effects) && is.character(random.effects))
        random.effects <- meffil:::extract.from.samplesheet(qc.objects, random.effects)
    
    old.pca <- meffil.pcs(qc.objects[-new.idx],
                          fixed.effects[-new.idx,,drop = F],
                          random.effects[-new.idx,, drop = F])
    new.controls <- meffil.control.matrix(qc.objects, normalize = T,
                                          fixed.effects = fixed.effects,
                                          random.effects = random.effects)[new.idx,]
    new.pca <- predict(old.pca, newdata = data.frame(new.controls, check.names = F))
    meffil:::pca.to.design.matrix(new.pca, number.pcs,
                                  fixed.effects[new.idx,, drop = F],
                                  random.effects[new.idx,, drop = F])
}




#' From PCA applied to the raw control matrix,
#' derive the corresponding design matrix using the first \code{number.pcs} principal components.
pca.to.design.matrix <- function(pca.ret, number.pcs, fixed.effects=NULL, random.effects=NULL) {
    if (class(pca.ret) != "matrix")
        pca.ret <- pca.ret$x
    if (missing(number.pcs)) 
        number.pcs <- ncol(pca.ret)
    stopifnot(number.pcs >= 1 && number.pcs <= ncol(pca.ret))
    control.components <- pca.ret[, 1:number.pcs, drop = F]
    design.matrix <- model.matrix(~control.components - 1)
    rownames(design.matrix) <- rownames(pca.ret)
    list(fixed = cbind(design.matrix, fixed.effects), random = random.effects)   
}

impute.matrix <- function(x, margin=1, fun=function(x) mean(x, na.rm=T)) {
    if (margin == 2) x <- t(x)
    
    idx <- which(is.na(x) | !is.finite(x), arr.ind=T)
    if (length(idx) > 0) {
        na.idx <- unique(idx[,"row"])
        v <- apply(x[na.idx,],margin,fun) ## v = summary for each row
        v[which(is.na(v))] <- fun(v)      ## if v[i] is NA, v[i] = fun(v)
        x[idx] <- v[match(idx[,"row"],na.idx)] ##

        stopifnot(all(!is.na(x)))
    }

    if (margin == 2) x <- t(x)
    x
}

#' Infinium HumanMethylation450 BeadChip normalization design matrix
#'
#' Design matrix derived by applying principal components analysis to control probes.
#'
#' @param qc.objects A list of outputs from \code{\link{meffil.create.qc.object}()}.
#' @param number.pcs Number of principal components to include in the design matrix (Default: all).
#' @param fixed.effects Names of columns in samplesheet that should be included as fixed effects
#' along with control matrix principal components (Default: NULL).
#' @param random.effects Names of columns in samplesheet that should be included as random effects
#' (Default: NULL).
#' @return Design matrix with one column for each of the first \code{number.pcs} prinicipal
#' components.
#'
#' @export
meffil.design.matrix <- function(qc.objects, number.pcs, fixed.effects = NULL, random.effects = NULL){
    if (!is.null(fixed.effects) && is.character(fixed.effects))
        fixed.effects <- extract.from.samplesheet(qc.objects, fixed.effects)
    if (!is.null(random.effects) && is.character(random.effects))
        random.effects <- extract.from.samplesheet(qc.objects, random.effects)
    
    pca.ret <- meffil.pcs(qc.objects, fixed.effects, random.effects)
    pca.to.design.matrix(pca.ret, number.pcs, fixed.effects, random.effects)
}

#' Calculate control probe PCs
#'
#' @param qc.objects A list of outputs from \code{\link{meffil.create.qc.object}()}.
#' @param fixed.effects Names of columns in samplesheet that should be included as fixed effects
#' along with control matrix principal components (Default: NULL).
#' @param random.effects Names of columns in samplesheet that should be included as random effects
#' (Default: NULL).
#' @return PCA of control probes
#' 
#' @export
meffil.pcs <- function(qc.objects, fixed.effects=NULL, random.effects=NULL) {
    control.matrix <- meffil.control.matrix(qc.objects, normalize=T,
                                            fixed.effects=fixed.effects,
                                            random.effects=random.effects)
    control.pca <- prcomp(control.matrix)
    return(control.pca)
}



#' Infinium HumanMethylation450 BeadChip control matrix
#'
#' Matrix containing control probe intensities from the Infinium HumanMethylation450 BeadChip.
#'
#' @param qc.objects A list of outputs from \code{\link{meffil.create.qc.object}()}.
#' @param normalize If \code{TRUE}, then control matrix is scaled and specified
#' fixed and random effects removed from the matrix.  Otherwise, the raw control matrix is returned.
#' (Default: \code{FALSE}).
#' @param fixed.effects Names of columns in samplesheet that should be included as fixed effects
#' along with control matrix principal components (Default: NULL).
#' @param random.effects Names of columns in samplesheet that should be included as random effects
#' (Default: NULL).
#' @return Matrix with one row per object consisting of control probe intensities and summaries.
#'
#' @export
meffil.control.matrix <- function(qc.objects, normalize=F,
                                  fixed.effects = NULL, random.effects = NULL) {
    if (!is.null(fixed.effects) && is.character(fixed.effects))
        fixed.effects <- extract.from.samplesheet(qc.objects, fixed.effects)
    if (!is.null(random.effects) && is.character(random.effects))
        random.effects <- extract.from.samplesheet(qc.objects, random.effects)
    
    stopifnot(all(sapply(qc.objects, meffil:::is.qc.object)))
    stopifnot(length(qc.objects) >= 2)
    stopifnot(is.null(fixed.effects) || nrow(fixed.effects) == length(qc.objects))
    stopifnot(is.null(random.effects) || nrow(random.effects) == length(qc.objects))

    control.matrix <- t(sapply(qc.objects, function(object) object$controls))

    if (normalize) {
        control.matrix <- meffil:::impute.matrix(control.matrix, margin = 2)
        control.matrix <- scale(control.matrix)
        control.matrix[control.matrix > 3] <- 3
        control.matrix[control.matrix < -3] <- -3
        if (!is.null(fixed.effects) || !is.null(random.effects)) {
            control.matrix <- meffil:::adjust.columns(control.matrix, fixed.effects, random.effects)
            control.matrix <- scale(control.matrix)
        }
    }
    control.matrix
}

#' Given a partition of the QC objects into two sets,
#' use the design matrix from one partition to predict the design
#' matrix for the other partition.
#' @param qc.objects A list of outputs from \code{\link{meffil.create.qc.object}()}.
#' @param number.pcs Number of principal components to include in the design matrix (Default: all).
#' @param new.idx Indices into \code{qc.objects} that identify the objects for which
#' the design matrix will be predicted.  The initial design matrix will be derived from
#' the other objects.
#' @param fixed.effects Names of columns in samplesheet that should be included as fixed effects
#' along with control matrix principal components (Default: NULL).
#' @param random.effects Names of columns in samplesheet that should be included as random effects
#' (Default: NULL).
#' @return Design matrix with one column for the \code{new.idx} partition predicted from the
#' first \code{number.pcs} prinicipal components of the objects. 
predict.design.matrix <- function(qc.objects, number.pcs, new.idx, fixed.effects = NULL, random.effects=NULL) {
    stopifnot(all(sapply(qc.objects, meffil:::is.qc.object)))
    stopifnot(all(new.idx %in% 1:length(qc.objects)))
    
    if (!is.null(fixed.effects) && is.character(fixed.effects))
        fixed.effects <- extract.from.samplesheet(qc.objects, fixed.effects)
    if (!is.null(random.effects) && is.character(random.effects))
        random.effects <- extract.from.samplesheet(qc.objects, random.effects)

    old.pca <- meffil.pcs(qc.objects[-new.idx],
                          fixed.effects[,-new.idx,drop=F],
                          random.effects[,-new.idx,drop=F])

    new.controls <- meffil.control.matrix(qc.objects, normalize=T,
                                          fixed.effects=fixed.effects,
                                          random.effects=random.effects)[new.idx,]

    new.pca <- predict(old.pca, newdata=data.frame(new.controls, check.names=F))

    meffil:::pca.to.design.matrix(new.pca, number.pcs,
                                  fixed.effects[,new.idx,drop=F],
                                  random.effects[,new.idx,drop=F])
}

#' From PCA applied to the raw control matrix,
#' derive the corresponding design matrix using the first \code{number.pcs} principal components.
pca.to.design.matrix <- function(pca.ret, number.pcs, fixed.effects=NULL, random.effects=NULL) {
    if (class(pca.ret) != "matrix")
        pca.ret <- pca.ret$x
    if (missing(number.pcs)) 
        number.pcs <- ncol(pca.ret)
    stopifnot(number.pcs >= 1 && number.pcs <= ncol(pca.ret))
    control.components <- pca.ret[, 1:number.pcs, drop = F]
    design.matrix <- model.matrix(~control.components - 1)
    rownames(design.matrix) <- rownames(pca.ret)
    list(fixed = cbind(design.matrix, fixed.effects), random = random.effects)   
}

impute.matrix <- function(x, margin=1, fun=function(x) mean(x, na.rm=T)) {
    if (margin == 2) x <- t(x)
    
    idx <- which(is.na(x) | !is.finite(x), arr.ind=T)
    if (length(idx) > 0) {
        na.idx <- unique(idx[,"row"])
        v <- apply(x[na.idx,],margin,fun) ## v = summary for each row
        v[which(is.na(v))] <- fun(v)      ## if v[i] is NA, v[i] = fun(v)
        x[idx] <- v[match(idx[,"row"],na.idx)] ##

        stopifnot(all(!is.na(x)))
    }

    if (margin == 2) x <- t(x)
    x
}
