Question: bumphunter for 3 (or more) groups
0
gravatar for alfonsosaera
11 weeks ago by
IISPV/Tarragona/Spain
alfonsosaera0 wrote:

Dear BioStar Users,

I am using minfi to analyze the results of an EPIC array. To check differential methylation I decided to look for differentially methylated regions (DMRs) using bumphunter as recommended.

bumphunter function returns candidate DMRs with several columns of information including value and area (clarification in this biostars post).

For 2 groups, value column is the average of the coefficients of the linear model calculated for each site within the DMR while area column is the sum of all coefficients. This can be easily checked by calculating the difference, between the 2 groups, of the mean value (M or beta) for each methylation site of the DMR and then averaging (value column) or summing (area column) these differences.

Thus, a positive value indicates more methylation in the non-control group, and the opposite for a negative value. Until here, everything is clear and makes sense.

My problem is that I have 3 groups, say control, drug A and drug B for example. But, bumphunter still returns a unique value and area. How do you interpret the value and area of each candidate DMR? What do they mean in terms methylation levels in each group? How would you filter them to get biological information?

To answer these questions I went to check the function code but got more confused.

bumphunter function can be accessed by typing getMethod("bumphunter", "matrix") and you get:

function (object, ...) 
{
  .local <- function (object, design, chr = NULL, pos, cluster = NULL, 
                      coef = 2, cutoff = NULL, pickCutoff = FALSE, pickCutoffQ = 0.99, 
                      maxGap = 500, nullMethod = c("permutation", "bootstrap"), 
                      smooth = FALSE, smoothFunction = locfitByCluster, useWeights = FALSE, 
                      B = ncol(permutations), permutations = NULL, verbose = TRUE, 
                      ...) 
  {
    nullMethod <- match.arg(nullMethod)
    if (missing(design)) 
      stop("design must be specified")
    if (missing(pos)) 
      stop("If object is a matrix, pos must be specified")
    bumphunterEngine(object, design = design, chr = chr, 
                     pos, cluster = cluster, coef = coef, cutoff = cutoff, 
                     pickCutoff = pickCutoff, pickCutoffQ = pickCutoffQ, 
                     maxGap = maxGap, nullMethod = nullMethod, smooth = smooth, 
                     smoothFunction = smoothFunction, useWeights = useWeights, 
                     B = B, permutations = NULL, verbose = verbose, ...)
  }
  .local(object, ...)
}

So we find that it is, basically, a S4 method to call the bumphunterEngine function that does most of the job.

So, typing bumphunterEngine you get:

function (mat, design, chr = NULL, pos, cluster = NULL, coef = 2, 
          cutoff = NULL, pickCutoff = FALSE, pickCutoffQ = 0.99, maxGap = 500, 
          nullMethod = c("permutation", "bootstrap"), smooth = FALSE, 
          smoothFunction = locfitByCluster, useWeights = FALSE, B = ncol(permutations), 
          permutations = NULL, verbose = TRUE, ...) 
{
  ...

  if (useWeights & smooth) {
    tmp <- .getEstimate(mat = mat, design = design, coef = coef, 
                        full = TRUE)
    rawBeta <- tmp$coef
    weights <- 1/tmp$sigma
    rm(tmp)
  }
  else {
    rawBeta <- .getEstimate(mat = mat, design = design, coef = coef, 
                            full = FALSE)
    weights <- NULL
  }

  ...

  else {
    beta <- rawBeta
    Index <- seq(along = beta)
  }

  ...

  comp <- computation.tots2(tab = tab, A = A)
  rate2 <- comp$rate2
  pvalues2 <- comp$pvalues2
  tab$p.value <- pvalues1
  tab$fwer <- rate1
  tab$p.valueArea <- pvalues2
  tab$fwerArea <- rate2
  tab <- tab[order(tab$fwer, -tab$area), ]
  algorithm <- list(version = packageDescription("bumphunter")$Version, 
                    coef = coef, cutoff = cutoff, pickCutoff = pickCutoff, 
                    pickCutoffQ = pickCutoffQ, smooth = smooth, maxGap = maxGap, 
                    nullMethod = nullMethod, B = B, permutations = permutations, 
                    useWeights = useWeights, smoothFunction = deparse(substitute(smoothFunction)))
  ret <- list(table = tab, coef = rawBeta, fitted = beta, pvaluesMarginal = pvs, 
              null = list(value = V, length = L), algorithm = algorithm)
  class(ret) <- "bumps"
  return(ret)
}

It is a long function, but the relevant parts for this problem IMO are coef = 2, in the parameters of the function, and .getEstimate(mat = mat, design = design, coef = coef, ...) which calculates the coefficients.

Regarding coef = 2, we can see that the original bumphunter function already has it defined as defaults. Also, knowing that bumhunter was originally defined to compare 2 groups I went to check the minfi wrapper to see if they implemented a change.

setMethod(
  "bumphunter",
  signature(object = "GenomicRatioSet"),
  function(object, design, cluster = NULL, coef = 2, cutoff = NULL,
           pickCutoff = FALSE, pickCutoffQ = 0.99, maxGap = 500,
           nullMethod = c("permutation", "bootstrap"), smooth = FALSE,
           smoothFunction = locfitByCluster, useWeights = FALSE,
           B = ncol(permutations), permutations = NULL, verbose = TRUE,
           type = c("Beta","M"), ...) {

    .isMatrixBackedOrStop(object, "bumphunter")

    type <- match.arg(type)
    bumphunterEngine(
      mat = getMethSignal(object, type),
      design = design,
      chr = as.factor(seqnames(object)),
      pos = start(object),
      cluster = cluster,
      coef = coef,
      cutoff = cutoff,
      pickCutoff = pickCutoff,
      pickCutoffQ = pickCutoffQ,
      maxGap = maxGap,
      nullMethod = nullMethod,
      smooth = smooth,
      smoothFunction = smoothFunction,
      useWeights = useWeights,
      B = B,
      permutations = permutations,
      verbose = verbose,
      ...)
  }
)

We see that the minfi wrapper basically fits the right matrix of data (mat = getMethSignal(object, type)) to dhe bumphunterEngine function keeping the coef = 2 option.

At this point, I went to check the other function .getEstimate(mat = mat, design = design, coef = coef, ...) which calculates the coefficients, if I understood properly, based on the data matrix (mat = mat), the experimental design (design = design) and the coefficient (coef = coef, which gets set to 2). The code can be found here.

.getEstimate <- function(mat, design, coef, B=NULL, permutations=NULL,
                         full=FALSE) {
  v <- design[,coef]
  A <- design[,-coef, drop=FALSE]
  qa <- qr(A)
  S <- diag(nrow(A)) - tcrossprod(qr.Q(qa))

  vv <- if (is.null(B)) matrix(v, ncol=1) else{
    if(is.null(permutations)){
      replicate(B, sample(v))
    } else{
      apply(permutations,2,function(i) v[i])
    }
  }


  sv <- S %*% vv
  vsv <- diag(crossprod(vv,sv))

  b <- (mat %*% crossprod(S, vv)) / vsv
  if(!is.matrix(b))
    b <- matrix(b, ncol = 1)
  if (full) {
    sy <- mat %*% S
    df.residual <- ncol(mat) - qa$rank - 1
    if (is.null(B)) {
      sigma <- matrix(sqrt(rowSums((sy - tcrossprod(b, sv))^2) / df.residual), ncol=1)
    } else {
      sigma <- b
      tmp <- sy
      for (j in 1:B) {
        tmp <- tcrossprod(b[,j], sv[,j])
        sigma[,j] <- rowSums((sy-tmp)^2)
      }
      sigma <- sqrt(sigma/df.residual)
    }
    out <- list(coef=b,
                sigma=sigma,
                stdev.unscaled=sqrt(1/vsv),
                df.residual=df.residual)
    if (is.null(B)) out$stdev <- as.numeric(out$stdev)
  } else {
    out <- b
  }
  return(out)
}

and my design would be something like:

des.groups <- c(rep("control", 2), rep("drug.A", 2), rep('drug.B', 2))
design <- model.matrix(~ des.groups)
design

  (Intercept) des.groupsdrug.A des.groupsdrug.B
1           1                0                0
2           1                0                0
3           1                1                0
4           1                1                0
5           1                0                1
6           1                0                1

I am worried that, just at the beggining of the .getEstimate function, it subsets the design based on the coef argument

v <- design[,coef]

getting

> v
1 2 3 4 5 6 
0 0 1 1 0 0

However this is as far as I can go with my linear algebra knowledge. So I am wondering if the value and area returned by bumphunter reflect the mehylation changes of "just" 2 groups (in my case control (Intercept) and des.groupsdrug.A).

What do you think? If this is the case, how would you address the comparison among 3 groups? with 3 pairwise comparisons?

If you arrived here, THANKS for reading the long post and please share your thoughts.

ADD COMMENTlink modified 11 weeks ago • written 11 weeks ago by alfonsosaera0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1676 users visited in the last hour