The results of FindMarkers and FoldChange are inconsistent

[caodudu]

I know the value of "avg_log2FC" calculated by the function FindMarker is actually calculated by FoldChange. But when I tried to calculate with default parameters ( I always only set ident.1 and ident.2.), I found that the results of two ways were inconsistent strangely.
When I tried in my other data, the results of two ways were consistent surprisingly.
Can the dear author or anyone tell me the potential reason or give me some advice?

[lucygarner]

I have the same issue, as indicated in #6773.

[caodudu]

I have the same issue, as indicated in #6773.
I notice that these three issues (#6654, #6701 and #6773) are now being discussed together. And now I know more about the problem which will be explained in detail below.

1. Conclusion
   The underlying cause of the inconsistencies and complexity of FoldChange and FindMarkers is the ambiguity of key code within the authors' functions.

2. Preparation
   To really understand and solve this problem, it's a good idea to query the source code for several functions.
   Please confirm the version(update to version 4.3) of Seurat first. This is because there are differences between versions of the source code due to Seurat's constant updates.

library(Seurat)
packageVersion("Seurat")

The functions we need to query are:

getAnywhere(FoldChange.default)
getAnywhere(FoldChange.Assay)
getAnywhere(FoldChange.Seurat)
getAnywhere(FindMarkers.default)
getAnywhere(PerformDE)
getAnywhere(WilcoxDETest)
getAnywhere(FindMarkers.Assay)
getAnywhere(FindMarkers.Seurat)

It seems to be a lot of functions, but please don't be surprised or intimidated by this. Please follow me to understand their nested relationships, the key to understanding the problem.

3. Analysis
   3.1 FoldChange.default is the most important function at the bottom level
   The core script is:

data.1 <- mean.fxn(object[features, cells.1, drop = FALSE])
data.2 <- mean.fxn(object[features, cells.2, drop = FALSE])
fc <- (data.1 - data.2)

It provides a basic model for comparing expression values between two groups, where mean.fxn refers to how to calculate the mean expression for a group of cells. The reason for subtraction is that, by default, the ideal mean.fxn should contain logarithmic processing internally. Anyway, we should realize that mean.fxn doesn't specify an average algorithm right now.

Here is a preliminary response to the question in github Issue #6654.
Mr. nathanhaigh has a legitimate question about calculating group-level mean value, but there seems to be some confusion. A clearer pattern is to consider analysis as two separate steps, including how to compute representative expression values at the group level (step 1), and how to find FC between groups (step 2).

3.2 FoldChange.Assay invoke FoldChange.default
The core script is:

 default.mean.fxn <- function(x) {
        return(log(x = rowMeans(x = x) + pseudocount.use, base = base))
    }
mean.fxn <- mean.fxn %||% switch(EXPR = slot, data = switch(EXPR = norm.method %||% 
        "", LogNormalize = function(x) {
        return(log(x = rowMeans(x = expm1(x = x)) + pseudocount.use, 
            base = base))
    }, default.mean.fxn), scale.data = rowMeans, default.mean.fxn)

The core script of 3.2 implements the solution of step 1 in 3.1 and invokes FoldChange.default to implement the step 2 in 3.1.

Internally, this function provides several ways to compute representative expression values at the group level by case, i.e. different mean.fxn.
The structure is shown as follows:

slot--'data'--norm.method--'LogNormalize'----------strategy1{expm1→mean→log2(x+1)}
      │                    │
      │                    └'others'---------------strategy2{mean→log2(x+1)}
      └'scale.data'--------------------------------strategy3{mean}
      │
      └'others'------------------------------------strategy2{mean→log2(x+1)}

#6654 should think that we should implement strategy3 instead of strategy1 to calculate the group mean value in the case of slot='data',norm.method=LogNormalize'
You can also look for the norm.method parameter, which is the key to causing FoldChange and FindMarkers to conflict. Since norm.method is a non-user-specified built-in variable, FoldChange will default to blank and not prompt, so strategy2 is automatically selected. However, FindMarkers will test the command variable to determine which strategy is suitable for norm.method. We will explain later.

3.3 FoldChange.Seurat invokes FoldChange.Assay
The author should unify FoldChange and FindMarkers by adding the judgement of Command variable at this step.

3.4 FindMarkers.default
The default parameter test.use for this function is wilcox, which is implemented by calling PerformDE.
The PerformDE function invokes WilcoxDETest to compute the p value. The result is called de.results but this step does not involve computing FC or logFC.

3.5 FindMarkers.Assay invokes FindMarkers.default and FoldChange.Assay
The core script is:

fc.results <- FoldChange(object = object, slot = data.slotcells.1 = cells.1, cells.2 = cells.2, features = features,pseudocount.use = pseudocount.use, mean.fxn = mean.fxn,fc.name = fc.name, base = base, norm.method = norm.method)
de.results <- FindMarkers(object = data.use, slot = data.slotcounts = counts, cells.1 = cells.1, cells.2 = cells.2,features = features, logfc.threshold = logfc.threshold,test.use = test.use, min.pct = min.pct,...)

It's confirmed that the fc.results were computed by calling FoldChange.Assay independently.
Note that FoldChange entered the norm.method parameter now, indicating that FindMarkers considers it in greater detail.

3.6 FindMarkers.Seurat
The core script is:

norm.command <- paste0("NormalizeData.assay)
norm.method <- if (norm.command %in% Command(object = object) && is.null(x = reduction)){Command(object = object, command = norm.command, value = "normalization.method")}

This script highlights that FindMarkers.Seurat queries the record from Command to confirm the normalization of the object, while passing parameters to the internally invoked FoldChange function.

FoldChange itself does not carry such similar judgements. Hopefully this explanation will be helpful for #6773, #6701 and similar issues.

4. Suggestions
   4.1 If you are only considering how to match FoldChange to FindMarkers, setting aside the question of whether the FC strategy outlined in Incorrect/Inconsistent Fold Change Calculation #6654 makes sense, it is advisable to specify the mean.fxn parameter like this:

fc=FoldChange(object=data,ident.1=group1,ident.2=group2,mean.fxn = function(x) {
	return(log(x = rowMeans(x = expm1(x = x)) + 1, base =2))
})

4.2 At times, FindMarkers have shown similar errors, especially when using normalized matrix introduced by someone else or other tools, which may not be recorded in command.
It is best to check the command or specify the mean.fxn function.

4.3 Since FoldChange and FindMarkers calculate logarithmic values with a default base of 2, instead of the natural logarithmic base that used in NormalizeData, which can lead to changes in the logFC value. A base of 2 can lead to overestimation of the logFC threshold. The logFC threshold should be chosen carefully, not mechanically as 1 or 2.

4.4 It is worth discussing whether strategy1 should be replaced by strategy3, see #6654 for details.

4.5 The default method of differential gene expression analysis in Seurat is the Wilcox rank sum test. The widespread use of Seurat package may result in the abuse of this rough method, in spite of nine alternative internal methods.
It is emphasized here that the Wilcox rank sum test is suitable for cell-level differential gene expression analysis, but not for sample-level differential gene expression analysis. The reason behind this is related to pseudoreplication, for details see https://www.nature.com/articles/s41467-021-25960-2 and https://www.nature.com/articles/s41467-021-21038-1

[nathanhaigh]

It would be good to have minimal working mock dataset(s) together with some testthat tests to reproduce and expose these issues. This would make it easier to expose the bug(s), determine what the correct behaviour should be and implement fix(es).

[lucygarner]

Thank you @caodudu and @nathanhaigh. I was wondering how fold change is calculated for scRNA-seq data using scanpy and came across the following: scverse/scanpy#673. And the function is here: https://github.com/scverse/scanpy/blob/ada761bc5c76007f0afc5a3516f1d05db99ae9ad/scanpy/tools/_rank_genes_groups.py
np.log2((expm1_func(mean_in_cluster) + 1e-9)/(expm1_func(mean_out_cluster) + 1e-9))

@caodudu, your suggestion to define the mean.fxn for FoldChange led to consistent results from FoldChange and FindMarkers.

@nathanhaigh, how would you recommend to set the mean.fxn. As:
mean.fxn = function(x) {rowMeans(log(expm1(x = x) + 1, base = 2))}

This gave me reduced values for the fold changes e.g. see CCL3 in the two tables below (upper table is with the mean.fxn used by FindMarkers as provided by @caodudu and lower table is using the approach above).

These fold changes are still considerably higher than obtained with the FoldChange function using default mean.fxn (see below).

@caodudu, what is FoldChange using as mean.fxn by default? I am using the "RNA" assay.

[caodudu]

@lucygarner In the 'RNA' assay, the default strategy is present like this:

slot--'data'--norm.method--'LogNormalize'----------strategy1{expm1→mean→log2(x+1)}
      │                    │
      │                    └'others'---------------strategy2{mean→log2(x+1)}
      └'scale.data'--------------------------------strategy3{mean}
      │
      └'others'------------------------------------strategy2{mean→log2(x+1)}

The key point to choose mean.fxn function, is not depend on assay, but on the slot and the methods of normalization. For common users, if not specified, the slot data is selected for DEG analysis.

[caodudu]

@lucygarner Decording to your description, the strategy designed by scanpy is mean→expm1→log2(x+1), which is similar with the strategy1 I mentioned, expm1→mean→log2(x+1). That is why the results between seurat and scanpy is similar but different. At issue is when to take the average within a given group, whether to take the pre-logarithmic data or the post-logarithmic data. I'm sorry it's hard for me to give a recommendation

[saketkc]

Apologies everyone for the long wait.

I have just pushed 992c1a9 to develop branch which should fix this and related issues. Can you pull from the develop branch and re-open this issue if you are still facing issues?

devtools::install_github("satijalab/seurat", ref="develop")

This should also be in the next Seuratv4 release.

[caodudu]

Apologies everyone for the long wait.

I have just pushed 992c1a9 to develop branch which should fix this and related issues. Can you pull from the develop branch and re-open this issue if you are still facing issues?

devtools::install_github("satijalab/seurat", ref="develop")

This should also be in the next Seuratv4 release.

@saketkc Thank you. But the FindMarkers function in Seurat V5 has bugs too. This can be attributed to the fact that the data structure and classes designed by BPCells do not fit the findmarkers function.