getMarkerFeatures's Log2FC changes dramatically by testMethod....

[dmsalsgh97]

Hi, I'm using ArchR's getMarkerFeatures function for finding DARs between two clusters.

First, When I used testMethod=wilcoxson with binarization, the log2FC of my interest is below:

>markerTest_wilcox <- getMarkerFeatures(
        ArchRProj = ArchR_Snap, 
        useMatrix = "PeakMatrix",
        groupBy = "Snap_Celltype",
        maxCells = 2000, 
        binarize = TRUE,
        testMethod = "wilcoxon",
        bias = c("TSSEnrichment", "log10(nFrags)"),
        useGroups = "cluster1", 
        bgdGroups = "cluster2"
        )

>test <- getMarkers(markerTest_wilcox, cutOff = "FDR <= 1 & Log2FC >= 0", returnGR = TRUE)
>TNF <- GRanges(seqnames = c('chr6'), ranges = IRanges(start = c(31573565), end = c(31580336)))
>test[[1]][queryHits(findOverlaps(test[[1]], GenomicRanges::intersect(test[[1]], TNF)))]

GRanges object with 2 ranges and 3 metadata columns:
      seqnames            ranges strand |            Log2FC                 FDR
         <Rle>         <IRanges>  <Rle> |         <numeric>           <numeric>
  [1]     chr6 31578398-31578898      * |  1.14630033094542 0.00107059230232709
  [2]     chr6 31575251-31575751      * | 0.214507084640496   0.132129724391847
                MeanDiff
               <numeric>
  [1] 0.0147868795573871
  [2] 0.0017448363528877
  -------

Second, When I used testMethod=binomial with binarization, the log2FC of my interest is dramatically changed as below:

>markerTest_bi <- getMarkerFeatures(
        ArchRProj = ArchR_Snap, 
        useMatrix = "PeakMatrix",
        groupBy = "Snap_Celltype",
        maxCells = 2000, 
        binarize = TRUE,
        testMethod = "binomial",
        bias = c("TSSEnrichment", "log10(nFrags)"),
        useGroups = "cluster1", 
        bgdGroups = "cluster2"
        )

>rm(test)
>test <- getMarkers(markerTest_bi, cutOff = "FDR <= 1 & Log2FC >= 0", returnGR = TRUE)
>TNF <- GRanges(seqnames = c('chr6'), ranges = IRanges(start = c(31573565), end = c(31580336)))
>test[[1]][queryHits(findOverlaps(test[[1]], GenomicRanges::intersect(test[[1]], TNF)))]

GRanges object with 2 ranges and 3 metadata columns:
      seqnames            ranges strand |            Log2FC                 FDR
         <Rle>         <IRanges>  <Rle> |         <numeric>           <numeric>
  [1]     chr6 31578398-31578898      * |  1.49185309632967 1.0598173894274e-08
  [2]     chr6 31575251-31575751      * | 0.938599455335857  0.0105848574337957
                MeanDiff
               <numeric>
  [1] 0.0190039318479685
  [2] 0.0072083879423329
  -------

My first thought was the different picking of background cells, but When I checked selected background cells, there was no difference...

>length(intersect(metadata(markerTest_wilcox)$MatchInfo$matchbgd$'cluster1'$bgd, metadata(markerTest_bi)$MatchInfo$matchbgd$'cluster1'$bgd))
1516
>length(intersect(metadata(markerTest_wilcox)$MatchInfo$matchbgd$'cluster1'$cells, metadata(markerTest_bi)$MatchInfo$matchbgd$'cluster1'$cells))
1516

Any possible reasons for those differnces?
Thanks.

[rcorces]

This is because there are different normalization factors imposed on the matrix if you use the binomial test. See the code snippet below:

ArchR/R/MarkerFeatures.R

Lines 203 to 223 in 968e442

	if(is.null(normBy)){
	if(tolower(testMethod) == "binomial"){
	normFactors <- NULL
	}else{
	if(tolower(useMatrix) %in% c("tilematrix", "peakmatrix")){
	normBy <- "ReadsInTSS"
	normFactors <- getCellColData(ArchRProj, normBy, drop=FALSE)
	normFactors[,1] <- median(normFactors[,1]) / normFactors[,1]
	}else{
	normFactors <- scaleTo / mColSums
	normFactors <- DataFrame(normFactors)
	}
	}
	}else{
	if(tolower(normBy) == "none"){
	normFactors <- NULL
	}else{
	normFactors <- scaleTo / mColSums
	normFactors <- DataFrame(normFactors)
	}
	}

I believe if you set normBy = "none" instead of NULL, it should no longer show a difference in log2FC. That being said, setting normBy = none is likely not the right way to analyze your data and you should be very familiar with what these different options do it you deviate from the recommended options.

[dmsalsgh97]

Thanks for the detailed reply.
As you said, I tried wilcoxon test without normalization by ( normBy = "none") and it showed the exact same Log2FC.

In this case, I agree that the log2FC calculated by binomial test without normalization is the wrong way.
So, I'll use the default wilcoxon test instead without binarization for further analysis!

And I'm curious why current getMarkerFeatures() function sets normFactors = NULL when testMethod is "binomial".
Is there any reason for that?

Thanks.

[cacaonib]

Hi
I also met a similar issue to find the DARs that compared between WT and KO.
When I used binomial test method, I could get a number of DARs (3617 peaks, FDR =< 0.1, LogFC >=0.5)
but... when I used wilcoxon method, I only got the 16 peaks (FDR =< 0.1, LogFC >=0.5).

So, I'm wonder "Can I set the normFactor = "ReadsInTSS" using binomial Test?

[yangxueqianR]

Hi
I also have a similar issue. I find the DAR that compared celltype.
When I used wilcoxon method , I only got the 68 peaks in one cell type(FDR <= 0.05 & LogFC >= 0.5)
but when I used binomial test method , I get 7897 peaks in one cell type(FDR <= 0.05 & LogFC >= 0.5)

Except for this parameter, everything is the same. So I want to know what should I do choose.