Limma: Extracting results when the global option of the decideTests is used

Question

[Ridha]

Greetings! Hope you all doing well :)

I have a question regarding extracting differentially expressed genes identified with the global method of the decideTests function from limma. I have 7 contrasts that are related to each other. From the user guide, it is advised to adjust for BETWEEN cotnrasts as well in cases like mine using method="global". I found an answer in other posts that the results from decideTests can be extracted using the write.fit function from limma as follows :

write.fit(fit3, adjust="BH", method="global", file="globalDEresults.txt")

However, I have 2 questions about this method:

1) Using write.fit function outputs a large file for ALL contrasts with their log fold change, p values, t statistcs and adjusted p-values in ONE file. This requires more work to make it tidy per contrast. Normally speaking(if method=separate), topTable function would allow me to extract results per contrast by changing the coef option as follows:

dpe_ckd_alone<-topTable(fit3, coef="CKD_alonevsNoCOM",number = nrow(fit3))

So, my question is: is there a similar method to toptable when you are correcting for between contrasts using method=global from decidetests? Otherwise, how would you extract results(including log fold changes, p values etc) per contrasts easily and assign them to an object?

2) If I don't want to extract the results directly to excel or txt file, but would rather want to manipulate it in R before exporting it, Is it possible to assign the data to an object? because if I write the follwoing code :

my_df<-write.fit(fit, results=resultsG, adjust="BH", method="global", file="globalDEresults.txt")

I get an empty object, which is to be expected as write.fit is intended to export results and not to assign them to an object?

Thank you very much in advance! Best, Ridha

[Melisa]

I have exactly the same problem! I just want to get the genes from the decideTest using method=global!!

[Melisa]

I've got it! https://bioconductor.org/packages/devel/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow.html

here's the answer The top DE genes can be listed using topTreat for results using treat (or topTable for results using eBayes). By default topTreat arranges genes from smallest to largest adjusted p-value with associated gene information, log-FC, average log-CPM, moderated t-statistic, raw and adjusted p-value for each gene. The number of top genes displayed can be specified, where n=Inf includes all genes. Genes Cldn7 and Rasef are amongst the top DE genes for both basal versus LP and basal versus ML.

Answer 1

[Gordon Smyth]

Just save decideTests() to an object:

dt <- decideTests(fit3, method="global")
summary(dt)

Then each column of dt selects the genes that are DE for that contrast.

[Ridha]

Dear professor Gordon Smyth, thank you very much for your answer. I would like to use this opportunity to express my gratitude for all your efforts and scientific contributions in this field, we are so grateful for all of your efforts!

Regarding my question, perhaps I haven't formulated it correctly or clear enough. My question is about extracting NOT only the number of differentially expressed genes, but actually the results of differnential expression. I mean the log fold change of each gene, with the p adjusted value etc per contrast. So, the code you suggsetd, would give me the following in my case:

        dt <- decideTests(fit3, method="global")

            summary(dt)

            contrast 1   contrast 2   contrast 3   contrast 4    contrast 5

Down              72          72           52        110              53                  
NotSig            84         104          115        142             227                 
Up               201         181          190        105              77          

However, what I get using the toptable function as I wrote in the question would be similar to this(using one contrast as an example):

  dpe_craids<-topTable(fit3, coef="CRAIDSvsNoCOM",number = nrow(fit3))
  head(dpe_craids)

       hgnc_name     logFC   AveExpr        t       P.Value     adj.P.Val        B
         WFDC2 0.7802672 10.179460 23.46939 7.709328e-107 2.752230e-104 232.9594
         TNFRSF19 1.2463862  5.994430 22.54186  1.278411e-99  2.281963e-97 216.3990
         TGFBR2 0.8634129  7.544950 22.24706  2.327319e-97  2.769509e-95 211.2150
         TNFRSF10B 1.0812091  7.285260 21.80680  5.132057e-94  4.580361e-92 203.5468
         TNFRSF11A 1.1757852  6.955970 20.93647  1.595200e-87  1.138973e-85 188.6576
          EPHA2 0.7546227  4.609625 20.30856  6.105032e-83  3.632494e-81 178.1491

And this is what I am looking for per contrast instead of the write.fit function which gives me NO log fold changes. It does give me the following colnames p-values,p.adjusted, t statistics and coef values per contrast, and column named A,F, F.p.value, in the same file.

Am I missing something here? Perhaps the idea of usinng decidetests is that I should use the adjusted p values from decidetests when method =global, but the log fold change should be obtained by the toptable function?

[Ridha]

Update: I found answers to my questions in this question. Notes for possibly similar question in the future:

1) the columns that start with coef from the write.fit function ARE ACTUALLY the log fold change that you normally would get using the toptable function.

2) The column named A reffers to the same column named AveExpr(the average expression) from the toptable function.

3) I could not find an answer about how to output results of each contrast alone (like in toptable). So the solution to subset if you have many contrasts is to use the function grepl.

An example would be if your contrast of interest contains the word heart, then you can subset the results from the write.fit function. You will have to export it and then importing it using read.table for instance. like so :

 write.fit(fit3, adjust="BH", method="global", file="globalDEresults.txt")
 globalDEresults<- read.delim("C:/Users/r123k/Desktop/R results and codes/Attempt3/globalDEresults.txt")
 heart_df<-globalDEresults[,grepl("heart",colnames(globalDEresults))] 

Thanks anyways professor Gordon Smyth!

[1769mkc]

This is useful answer as I have used limma not so much but I have another question here you are have 7 contrasts as mentioned in the question so if you want to see the differential gene between 2 condition the way normally done using the contrast function of deseq2 how can i do the same in case of this object ? in my case I have mutation or mutated gene as my design matrix so I want to compare how genes are different between samples which carry two different mutation ?