Decidetest results summary different from results written as .csv file
1
0
Entering edit mode
Ankita • 0
@12419918
Last seen 11 months ago
India

Hi, I am trying to identify my DEG from the following code. However, my decide test summary results are different from the one written as csv file, i.e. different number of up and downregulated genes obtained. Please explain the difference.

#Fit linear model for each gene
GSE84054_fit <- lmFit(GSE84054_voomData, GSE84054_design) %>%
  contrasts.fit(GSE84054_contrasts) %>%
  treat(lfc = 1)

#identify DEGs
> GSE84054_results <- decideTests(GSE84054_fit, 
+                                 p.value = 0.05,
+                                 adjust.method = "fdr")
> summary(GSE84054_results)
       SphereVSTumor
Down             705
NotSig         12931
Up                48
#finally write results as
GSE84054_allDEresults <- topTreat(GSE84054_fit, 
                         coef = "SphereVSTumor", 
                         number = Inf, 
                         adjust.method = "fdr") %>%
  as.data.frame() 

#We  filter the allDEresults table to just have the ones which we define to be significantly differentially expressed. In this case, significant differential expression means that genes must have FDR-adjusted p-value < 0.05 and absolute log2 fold change greater than 1 (i.e. either logFC < -1 or logFC > 1)

>GSE84054_allDEresults <- allDEresults %>%
  dplyr::mutate(isSignificant = case_when(
    adj.P.Val < 0.05 & abs(logFC) > 1 ~ TRUE, 
    TRUE ~ FALSE # If conditions in the line above are not met, gene is not DE. 
  ))

>GSE84054_sigDEresults <- GSE84054_allDEresults %>%
  dplyr::filter(isSignificant == TRUE)

#We will export those two tables above into the output folder as CSV spreadsheets (this is useful for sharing data with collaborators)
>GSE84054_allDEresults %>% 
  readr::write_csv("GSE84054_all_DE_results.csv")

>GSE84054_sigDEresults %>%
  readr::write_csv("GSE84054_significant_DE_genes.csv")


sessionInfo( )
edgeR voom decidetests DEG • 553 views
ADD COMMENT
0
Entering edit mode
@gordon-smyth
Last seen 30 minutes ago
WEHI, Melbourne, Australia

Dear Ankita,

It is a bit unfair to post such a long block of code without proper explanation and expect us to decipher it. It would be better to ask a more focused and better explained question.

Anyway, looking over your code, it seems in the first half of your code that you have applied treat() to prioritise differentially expressed genes. That is a highly recommended limma pipeline, although the choice of lfc=1 is much higher than I recommend. Why are you not following the advice given in the treat help page, which is to set lfc to a low value?

On the second half of your code, you seem to have instead done your own ad hoc filtering of the DE genes with a logFC cutoff, something that the limma documentation advises you not to do.

The first approach is good. The second is bad. Why would you expect the two approaches to be the same?

Instead of posting a question here, you would learn more by reading the documentation about treat, including the reference that is cited on the treat help page.

ADD COMMENT

Login before adding your answer.

Traffic: 732 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6