saeyslab / nichenetr

NicheNet: predict active ligand-target links between interacting cells
467 stars 115 forks source link

Different results running. the same. code in different version. of #285

Closed pariaaliour closed 2 weeks ago

pariaaliour commented 1 month ago

Dear Saeyslab Team,

Thank you for the amazing package.

I am encountering an issue where running the same code with identical input data in both RStudio and on a server yields different numbers of geneset_oi. This discrepancy is quite confusing.

`seurat_obj_receiver= subset(seuratObj, idents = receiver)

condition_oi = "als" condition_reference = "con"

DE_table_receiver = FindMarkers(object = seurat_obj_receiver, ident.1 = condition_oi, ident.2 = condition_reference, group.by = "group_id", min.pct = 0.05) %>% rownames_to_column("gene")

geneset_oi = DE_table_receiver %>% dplyr::filter(p_val_adj <= 0.05 & abs(avg_log2FC) >= 0.25) %>% pull(gene)

geneset_oi = geneset_oi %>% .[. %in% rownames(ligand_target_matrix)]

length(geneset_oi)`

In Rstudio I have R version 4.3.3, seurat 4.4.0 and I get the below number > length(geneset_oi) [1] 683 In server, I have R version 4.1.2, seurat 5.1.0. and I get the below number of geneset_oi > length(geneset_oi) [1] 2522

I really appreciate it if you. have any input on this. Best, Paria

pariaaliour commented 1 month ago

I tried to uninstall and then, reinstall R and Rstudio. But the issue exists. I appreciate your help on this! Paria

Eisuan commented 1 month ago

Dear Paria,

the differences in the DEGs you observe are probably due to the different implementations of the Wilcoxon test/fold change computation in Seurat 5 vs 4. Seurat 5 uses Presto to guarantee a faster differential expression analysis. Here you can also find other users making observations of high discrepancies in the number of DEGs called after DEA performed using different Seurat versions. The GitHub issue I linked also pointed to differences in standard parameters. I suggest having a closer look at all the parameters of the FindMarkers function. It may be a pct threshold-related issue.

As for what concerns the NicheNet analysis, in our experience, the optimal number of DEGs ranges from 0.5 to 10% of the number of background-expressed genes. However, this does not guarantee that the main assumption of our method holds (i.e., the DEGs observed reflect the result of cell-cell communication).

Best regards,

Daniele

pariaaliour commented 1 month ago

Thank you so much for your reply. I got the point. So, there are two different changes: first, using Presto; second, the default for min.pct. I think if I set the min.pct = 0.1 with Seurat 5 I still get different results, the difference might be slighter, right? And, do you recommend using the old one? (I plan to use the newer version of Seurat. just wanna make sure I'm not missing anything) Thanks, Paria

Eisuan commented 1 month ago

Dear Paria, I also observed a significantly different number of DEGs when using Seurat V5 when I performed the DEA. However, this has nothing to do with NicheNet, and the only thing I can advise is caution on the criteria to call a gene DE, which can be extremely case-specific. Enrichment procedures could help retrieve an overall overview of the DEG results, as well as the presence of DEGs already known to have a role in the process investigated.

Bests, Daniele

csangara commented 1 month ago

Dear Paria,

It's not just the min.pct but also the function to calculate the mean expression, as seen in this line of code.

The difference in versions 4.4 and 5 is given below. If you use Seurat v5 and give the mean.fxn argument in FindMarkers as the one from v4.4, the result should be the same as when running Seurat v4.4. Also note that the pseudocount.use in v4 is 1 and in v5 is 0.1 (mentioned here).

# 4.4
mean.fxn <- function(x) log(x = rowMeans(x = expm1(x = x)) + pseudocount.use, base = base)

# 5.0.1
mean.fxn <- function(x) log(x = (rowSums(x = expm1(x = x)) + pseudocount.use)/NCOL(x), base = base)
pariaaliour commented 1 month ago

Thank you so much for clarifying this @csangara and @Eisuan! Paria