Question: limma roast results not in accordance with barcodeplot
1
gravatar for b.nota
12 months ago by
b.nota320
Netherlands
b.nota320 wrote:

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?

limma roast • 287 views
ADD COMMENTlink modified 12 months ago • written 12 months ago by b.nota320
Answer: C: limma roast results not in accordance with barcodeplot
1
gravatar for Gordon Smyth
12 months ago by
Gordon Smyth37k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth37k wrote:

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 COMMENTlink modified 12 months ago • written 12 months ago by Gordon Smyth37k

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 REPLYlink modified 12 months ago • written 12 months ago by b.nota320
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 16.09
Traffic: 299 users visited in the last hour