Question

limma roast results not in accordance with barcodeplot

1

Entering edit mode

b.nota ▴ 370

@bnota-7379

Last seen 4.2 years ago

Netherlands

I am trying to use limma roast, as advised in a previous post A: Compare groups of different RNAseq sets. However, I have doubts about if it is correct what I am doing, since the barcodeplot gives a different picture of my results. I have an RNA-seq experiment analyzed with limma trend, and want to see if a signature of 48 genes (from another experiment) is enriched herein. Here some essential parts of the code:

design <- model.matrix(~0 + targets$Cell + targets$Donor)

fit <- lmFit(logCPM, design)

contr <- makeContrasts(
hivslo = hi - lo,
hivsint = hi - int,
intvslo = int - lo,
levels = design)

fit2 <- eBayes(contrasts.fit(fit, contr), trend=TRUE)

results <- decideTests(fit2, method="separate")

summary(results)
       hivslo hivsint intvslo
Down     2461    1179      47
NotSig  11245   13746   15486
Up       1891     672      64

Signature <- "X:/names.txt"

signature <- scan(Signature, what = "character")

set.seed(12)
roast_hi_lo_signature <- roast(logCPM, index = signature, design = design,
contrast = contr[,1], trend = T)

roast_hi_lo_signature

         Active.Prop   P.Value
Down       0.2317073 0.8189095
Up         0.2926829 0.1815908
UpOrDown   0.2926829 0.3630000
Mixed      0.5243902 0.0070000

index.vector <- All_genes %in% signature

barcodeplot(fit2$t[,1], index = index.vector)

So here my question, in the results only the mixed group has significant enrichment (p < 0.05), however the barcode plot suggest a enrichment in down? So I think I am doing something wrong? Did I do the roast analysis wrong or the barcode?

roast limma • 1.8k views

ADD COMMENT • link 6.5 years ago b.nota ▴ 370

score 1 · Answer 1 · 2018-05-09

1

Entering edit mode

Gordon Smyth 51k

@gordon-smyth

Last seen 34 minutes ago

WEHI, Melbourne, Australia

I can't tell yet, from the code given, whether either of the roast or barcodeplot calls are quite correct. This is because you haven't made clear yet the relationship of the signature list to the gene annotation in fit2. It isn't clear that you are creating the index vectors correctly.

What is All_genes and how does it relate to the logCPM matrix? What are the row.names of logCPM? What type of IDs are in the signature list? Do they match the row.names of logCPM?

Normally I would expect to see

index.vector <- row.names(logCPM) %in% signature

It would be a good idea to run roast() and barcodeplot() on exactly the same index argument, for example by:

roast(logCPM, index.vector, design, contrast=contr[,1], trend=TRUE, nrot=10000)

Also you might try:

camera(logCPM, index.vector, design, contrast=contr[,1], trend=TRUE)

ADD COMMENT • link 6.5 years ago Gordon Smyth 51k

0

Entering edit mode

Thanks (again) for your help, it's really appreciated. It was indeed the index that was not correct. I used for roast only the ensemble gene names of interest (which worked, but was not correct). With barcodeplot this index gave an error, which is why I used the index vector like in your example index.vector <- row.names(logCPM) %in% signature, which gave the plot I have shown. Now I redid roast with that index.vector and it produces results which are in accordance with the plot.

index.vector <- All_genes %in% signature

set.seed(12)
roast_hi_lo_signature <- roast(logCPM, index = index.vector, design = design,
contrast = contr[,1], trend = T, nrot = 10000)

roast_hi_lo_signature
         Active.Prop    P.Value
Down       0.7195122 0.00339983
Up         0.1951220 0.99665017
UpOrDown   0.7195122 0.00679932
Mixed      0.9146341 0.00139986

barcodeplot(fit2$t[,1], index = index.vector)

ADD REPLY • link 6.5 years ago b.nota ▴ 370