`scran::findMarkers` always return all available markers as differential, albeit with different p-values

[denvercal1234GitHub]

Hi there,

Thank you for the package.

I applied scran::findMarkers to determine the differentially expressed proteins of my clusters after performing clustering of the spectral flow cytometry data with CATALYST as below.

QUESTION 1. I noticed for any cluster, there are the same number of markers returned in the table (in fact all markers available in the object), although these markers with different p-values. So, does scran::findMarkers simply rank the markers available with p-values being how much it differentiates a given cluster against all other clusters?

QUESTION 2. Is there a way to also show the logFC column in the result instead of pre-setting it with lfc=2 within the function? Because is it true that p-values in these tests are only rough estimates and not quite reliable statistically to define cluster biomarkers; thus using logFC is better?

####### Select only the protein markers of interest to be analysed so that the analysis will only consider these proteins
F37_sce_backboneClustering <- filterSCE(F37_sce_backboneClustering, channel_name %in% markers_to_analyzed_vector)

F37_sce_backboneClustering_DEPwilcoxUplfc2ptypeAllMarkers = scran::findMarkers(assay(F37_sce_backboneClustering), groups= F37_sce_backboneClustering_meta16_clusteringIDperCell, test.type="wilcox", direction="up", lfc=2, pval.type="all", sorted=T, add.summary = T)

Thank you for your help!

[LTLA]

1: findMarkers will report statistics for all genes. It's up to you to decide where you want to draw the line with respect to defining a "marker". For example, you could set a threshold on the FDR, or you could take the top 100 genes, etc.

2: The lfc= argument is a testing threshold (inspired by limma::treat), it only affects the p-value calculation.

In your case, the actual issue is that you've set test.type="wilcox", which is causing the function to report AUCs instead of log-fold changes. This is because the AUC is a more relevant effect size when using the Wilcoxon test. If you want log-fold changes, just call the function again with test.type="t", which reports the log-fold change for each pairwise comparison. (Best to do so with sorted=FALSE to ensure that the order of rows is the same between outputs.)

Alternatively, you may consider using scoreMarkers, which dispenses with p-values altogether and just reports a variety of useful effect sizes, including Cohen's d (i.e., the standardized log-fold change) and the AUC.

[denvercal1234GitHub]

Hi @LTLA -- Thank you so much for the response. Very helpful.

Does it mean it is more accurate to use AUC (instead of FC) when using wilcox? I am unclear because in scRNAseq DEG analysis, Wilcox and not t-test is used with reported logFC and not AUC.

In your experiences, should I sort the markers based on logFC from t-test, or based on AUC from Wilcox, because I am hesitant to rank/sort the markers by "p-values" and "FDR" given that the p-values in these kinds of analyses are inherently biased.

[LTLA]

Does it mean it is more accurate to use AUC (instead of FC) when using wilcox?

I don't think "accurate" is the right word here.

The AUC is just a more relevant statistic for the Wilcoxon test, as the p-value is effectively derived from the AUC. So if you get a low p-value, you can be sure that the AUC for at least one pairwise comparison is high. (Or multiple comparisons, depending on pval.type=.)

You can't say that about the log-fold change, because that statistic is not directly related to the p-value. For example, it is easy to get a very large log-fold change with a high p-value if you have a highly variable distribution with large outliers in one group. This inflates the log-fold change but is penalized by the tests.

In your experiences, should I sort the markers based on logFC from t-test, or based on AUC from Wilcox, because I am hesitant to rank/sort the markers by "p-values" and "FDR" given that the p-values in these kinds of analyses are inherently biased.

While it is true that the p-values are "biased" in the sense that they are too optimistic, they are still okay for ranking. For example, if you rank by AUC, you're basically ranking by the Wilcoxon p-value anyway. And ranking by the t-test p-value is still better than ranking by log-fold change if you want to consider the variance.

But you can basically rank by whatever you want. For example, in our app, users can select a variety of different strategies to rank marker genes, each with their own strengths and weaknesses. As long as you understand what those are, you can choose one that suits your needs.

[denvercal1234GitHub]

This is very helpful, @LTLA! Thank you.