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( )

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.