Search
Question: limma decideTests global method with multiple groups
0
gravatar for jamie.gearing
4 weeks ago by
jamie.gearing0 wrote:

Hello,

I have an experiment with six groups: WT and KO cells that are untreated, treated with IFN A or IFN B.
Here are some of the possible contrasts.

cm <- makeContrasts(

    # Treatment effects
      A.WT = WT.A - WT.Un,
      B.WT = WT.B - WT.Un,

      A.KO = KO.A - KO.Un,
      B.KO = KO.B - KO.Un,

    # Differences between treatments within strain
      A.vs.B.WT = WT.A - WT.B,
      A.vs.B.KO = KO.A - KO.B,

    # Basal difference between strains
      KO.vs.WT = KO.Un - WT.Un,

    # Treatment differences between strains
      A.KO.vs.WT = (KO.A - KO.Un) - (WT.A - WT.Un),
      B.KO.vs.WT = (KO.B - KO.Un) - (WT.B - WT.Un),

      levels = design
)

fit2 <- contrasts.fit(fit, cm)
fit2 <- treat(fit2, lfc = log2(1.2))
results <- decideTests(fit2)


One of the main questions of interest is what genes are up-regulated specifically by the treatments. e.g. in WT cells, genes that are significantly induced by A compared to unstimulated and to B, but not induced by B. The common specific genes in WT and KO cells are also of interest.

A.WT.specific <- results[,"A.WT"] == 1 & results[,"A.vs.B.WT"] == 1 & results[,"B.WT"] != 1
# Similarly for A.KO.specific, B.WT.specific and B.KO.specific

A.specific <- A.WT.specific & A.KO.specific
# Similarly for B.specific

I think this is a situation where using decideTests with method="global" would be appropriate, however I want to avoid including unnecessary contrasts.

I was wondering whether anyone could offer any advice about whether the contrast matrix above would be suitable to use with the "global" method or whether I would need to define smaller contrast matrices, such as those shown below, and run contrasts.fit, treat and decideTests again as appropriate.

Thanks very much.

# Keep both WT and KO together
cm.specific <- makeContrasts(

    # Treatment effects
      A.WT = WT.A - WT.Un,
      B.WT = WT.B - WT.Un,

      A.KO = KO.A - KO.Un,
      B.KO = KO.B - KO.Un,

    # Differences between treatments within strain
      A.vs.B.WT = WT.A - WT.B,
      A.vs.B.KO = KO.A - KO.B,

      levels = design
)

# Split even further into WT and KO contrast matrices
cm.wt.specific <- makeContrasts(

    # Treatment effects
      A.WT = WT.A - WT.Un,
      B.WT = WT.B - WT.Un,

    # Differences between treatments within strain
      A.vs.B.WT = WT.A - WT.B,

      levels = design
)
# cm.ko.specific similarly

 

ADD COMMENTlink modified 4 weeks ago by Aaron Lun17k • written 4 weeks ago by jamie.gearing0
1
gravatar for Aaron Lun
4 weeks ago by
Aaron Lun17k
Cambridge, United Kingdom
Aaron Lun17k wrote:

Yes, it would seem that method="global" would be most appropriate here, as this ensures that the same p-value threshold is effectively used in each contrast. Otherwise, the dependence of the FDR calculation on the p-values of other genes may result in different thresholds for significance between contrasts. In particular, genes with larger p-values may be significant in some contrasts if there are many other genes with low p-values; conversely, the same p-values may not be considered significant in other contrasts with few DE genes.

Setting a global correction avoids such counterintuitive behaviour. Obviously, you only want to apply global BH correction to the contrasts that are actually of relevance. There's no point in being affected by the p-values for contrasts that aren't going to be involved in your final selection of genes. A superficial examination of your analysis suggests that cm.specific is probably most suitable here. Note that you don't need to redefine the contrast matrix if you only want a few of its contrasts, you can just subset it to keep the columns of interest.

That being said, there's not much statistical rigour involved in ad hoc combinations of contrast results anyway. The simple task of intersecting two DE lists (each detected at a certain FDR threshold) may not necessarily yield an intersection with a FDR below the specified rate. Selecting genes that are not significant is even more troublesome, as this is not what the vast majority of hypothesis tests are designed for. Unfortunately, the statistically rigorous alternative (two-one sided tests) is so conservative as to be useless, so your current approach is probably the best of a bad lot. Just take the results with a grain of salt.

ADD COMMENTlink modified 4 weeks ago • written 4 weeks ago by Aaron Lun17k

Thanks Aaron. I was leaning towards that particular choice, but it was very helpful to get another opinion. Yes, I shall just subset fit2 for the decideTests and then go from there.

results.specific <- decideTests(fit2[,c("A.WT", "B.WT", "A.vs.B.WT", "A.KO", "B.KO", "A.vs.B.KO")], method = "global")

I had a suspicion that perhaps this was not entirely statistically rigorous, particularly the part about which genes are not significantly different,  but it was best way that I could think of. It seems to get me a few candidate genes and I suppose it is better than nothing.

Thanks again.

ADD REPLYlink written 4 weeks ago by jamie.gearing0
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.2.0
Traffic: 111 users visited in the last hour