Search
Question: limma roast results not in accordance with barcodeplot
0
gravatar for b.nota
6 months ago by
b.nota300
Netherlands
b.nota300 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?

ADD COMMENTlink modified 6 months ago • written 6 months ago by b.nota300
1
gravatar for Gordon Smyth
6 months ago by
Gordon Smyth35k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth35k 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 6 months ago • written 6 months ago by Gordon Smyth35k

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 6 months ago • written 6 months ago by b.nota300
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: 360 users visited in the last hour