Different logFC and FDR in pairwise testing using getMarkerFeatures() depending on order of groups

[mireiacodina]

Hello,

Thanks for developing ArchR. I have a question related to the getMarkerfeatures() function to do pairwise testing between groups. I want to compare the accessibility between two cell types and I get different logFC and FDR values depending on the group that I use in useGroups and the one that I set as bgdGroups.

Here is the code that I run:

### Pairwise comparison between groups
group1 <- "wt_ectoderm(50epi)"
group2 <- "wt_ventrolateral mesoderm(50epi)"

markerTest <- getMarkerFeatures(
  ArchRProj = proj50epi7, 
  useMatrix = "PeakMatrix",
  groupBy = "Clusters_RNA_sub",
  testMethod = "wilcoxon",
  bias = c("TSSEnrichment", "log10(nFrags)"),
  useGroups = group1,
  bgdGroups = group2
)


markerTest2 <- getMarkerFeatures(
  ArchRProj = proj50epi7, 
  useMatrix = "PeakMatrix",
  groupBy = "Clusters_RNA_sub",
  testMethod = "wilcoxon",
  bias = c("TSSEnrichment", "log10(nFrags)"),
  useGroups = group2,
  bgdGroups = group1
)

markerList <- getMarkers(markerTest, cutOff = "FDR <= 0.05 & abs(Log2FC) >= 0.5", returnGR = FALSE)
sig_peaks <- as.data.frame(markerList)
head(sig_peaks)

group         group_name seqnames  idx    start      end    Log2FC          FDR   MeanDiff
32640      1 wt_ectoderm(50epi)    chr12 1859  8915090  8915590 -2.546727 2.346360e-32 -0.4123869
213971     1 wt_ectoderm(50epi)     chr5 8909 42280151 42280651 -2.397811 1.243829e-31 -0.4614626
201083     1 wt_ectoderm(50epi)     chr4 6087 21660738 21661238 -3.300984 1.895344e-28 -0.4402070
41075      1 wt_ectoderm(50epi)    chr13 1104  4433536  4434036 -2.909491 2.280093e-27 -0.4121067
26517      1 wt_ectoderm(50epi)    chr11 5323 25231320 25231820 -2.248990 1.307753e-25 -0.4673072
226621     1 wt_ectoderm(50epi)     chr6 6507 32164076 32164576 -1.976865 7.026898e-25 -0.5328489


markerList2 <- getMarkers(markerTest2, cutOff = "FDR <= 0.05 & abs(Log2FC) >= 0.5", returnGR = FALSE)
sig_peaks2 <- as.data.frame(markerList2)
head(sig_peaks2)
group                       group_name seqnames  idx    start      end   Log2FC          FDR  MeanDiff
32640      1 wt_ventrolateral mesoderm(50epi)    chr12 1859  8915090  8915590 3.773655 5.760450e-40 0.4569012
213971     1 wt_ventrolateral mesoderm(50epi)     chr5 8909 42280151 42280651 2.777479 6.533203e-37 0.4658591
226621     1 wt_ventrolateral mesoderm(50epi)     chr6 6507 32164076 32164576 2.183288 4.831842e-36 0.5751346
52253      1 wt_ventrolateral mesoderm(50epi)    chr14 1343  5557140  5557640 3.220382 1.345166e-34 0.3908514
41530      1 wt_ventrolateral mesoderm(50epi)    chr13 1559  6151759  6152259 2.018802 2.807008e-34 0.4459561
26517      1 wt_ventrolateral mesoderm(50epi)    chr11 5323 25231320 25231820 2.404376 1.429788e-33 0.4239108

Also, the total number of significant peaks based on the cutoff that I used is different (3601 for sig_peaks and 4889 for sig_peaks2). I would be very thankful if I could get any feedback on why the results are different. I also tried running the code without the "bias" variable but it did not change the outcome.

Thank you very much.

---

SessionInfo()
R version 4.2.3 (2023-03-15)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Ventura 13.3

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib

Random number generation:
 RNG:     L'Ecuyer-CMRG 
 Normal:  Inversion 
 Sample:  Rejection 
 
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggrepel_0.9.3                       presto_1.0.0                        nabor_0.5.0                         SeuratObject_4.1.3                  sp_1.6-0                           
 [6] BSgenome.Drerio.UCSC.danRer11_1.4.2 BSgenome_1.66.3                     rtracklayer_1.58.0                  Biostrings_2.66.0                   XVector_0.38.0                     
[11] dplyr_1.1.1                         tidyr_1.3.0                         rhdf5_2.42.0                        Rcpp_1.0.10                         Matrix_1.5-3                       
[16] data.table_1.14.8                   stringr_1.5.0                       plyr_1.8.8                          magrittr_2.0.3                      ggplot2_3.4.2                      
[21] gtable_0.3.3                        gtools_3.9.4                        gridExtra_2.3                       ArchR_1.0.2                         SummarizedExperiment_1.28.0        
[26] Biobase_2.58.0                      GenomicRanges_1.50.2                GenomeInfoDb_1.34.9                 IRanges_2.32.0                      S4Vectors_0.36.2                   
[31] BiocGenerics_0.44.0                 MatrixGenerics_1.10.0               matrixStats_0.63.0                 

loaded via a namespace (and not attached):
 [1] BiocManager_1.30.20      GenomeInfoDbData_1.2.9   Rsamtools_2.14.0         yaml_2.3.7               globals_0.16.2           pillar_1.9.0             lattice_0.20-45         
 [8] glue_1.6.2               digest_0.6.31            RColorBrewer_1.1-3       colorspace_2.1-0         htmltools_0.5.5          XML_3.99-0.14            pkgconfig_2.0.3         
[15] pheatmap_1.0.12          listenv_0.9.0            zlibbioc_1.44.0          purrr_1.0.1              scales_1.2.1             BiocParallel_1.32.6      tibble_3.2.1            
[22] farver_2.1.1             generics_0.1.3           withr_2.5.0              ROCR_1.0-11              cli_3.6.1                XLConnect_1.0.7          crayon_1.5.2            
[29] evaluate_0.20            parallelly_1.35.0        future_1.32.0            fansi_1.0.4              progressr_0.13.0         Cairo_1.6-0              tools_4.2.3             
[36] BiocIO_1.8.0             lifecycle_1.0.3          Rhdf5lib_1.20.0          kernlab_0.9-32           munsell_0.5.0            DelayedArray_0.24.0      compiler_4.2.3          
[43] rlang_1.1.0              RCurl_1.98-1.12          rhdf5filters_1.10.1      rstudioapi_0.14          rjson_0.2.21             labeling_0.4.2           bitops_1.0-7            
[50] rmarkdown_2.21           restfulr_0.0.15          codetools_0.2-19         R6_2.5.1                 GenomicAlignments_1.34.1 knitr_1.42               future.apply_1.10.0     
[57] fastmap_1.1.1            utf8_1.2.3               rJava_1.0-6              stringi_1.7.12           parallel_4.2.3           vctrs_0.6.1              tidyselect_1.2.0        
[64] xfun_0.38     

[rcorces]

I believe you are running into the same phenomenon that was described here and fixed a while back.
#1409

Specifically, the closest param

ArchR/R/MarkerFeatures.R

Line 44 in 629cf1d

#' @param closest A boolean value that indicated whether to use closest cells from foreground and background instead of random sampling

[mireiacodina]

Thank you very much for your swift reply.
I added the updated getMarkerFeatures() in my code and I could include the closest parameter. After using it, the results look much more similar, however, the total number of significant peaks is still different depending on the group that I use in useGroups and the one that I set as background. It is a different case form the one reported in #1409, since there it was reported different number of peaks when changing the order in the useGroups but using the exact same background group.

In my analysis, I just want to compare 2 groups between them, so I need to specify one as UseGroups and the other one as bgdGroups and, in this situation, the total number of peaks is still different, even when using closest = TRUE. Could it be somehow related to the normalization of the axes when using TSS score and log10(nFrags) that the structure might change depending on the group set as background?

Here is the code that I run with the closest param:

group1 <- "wt_margin(50epi)"
group2 <- "wt_ectoderm(50epi)"

markerTest <- getMarkerFeatures(
  ArchRProj = proj50epi7, 
  useMatrix = "PeakMatrix",
  groupBy = "Clusters_wt_margin",
  testMethod = "wilcoxon",
  bias = c("TSSEnrichment", "log10(nFrags)"),
  useGroups = group1,
  bgdGroups = group2,
  closest = TRUE
)

markerList <- getMarkers(markerTest, cutOff = "FDR <= 0.05 & abs(Log2FC) >= 0.5", returnGR = FALSE)
list1 <- as.data.frame(markerList)
head(list1) ### total number of significant peaks = 34,946

     group       group_name seqnames  idx    start      end   Log2FC           FDR  MeanDiff
201083     1 wt_margin(50epi)     chr4 6087 21660738 21661238 3.426848 7.657309e-167 0.5484994
121750     1 wt_margin(50epi)     chr2 7427 35671231 35671731 3.233131 8.314378e-162 0.5141475
57405      1 wt_margin(50epi)    chr14 6495 31980882 31981382 2.658234 1.290151e-158 0.5741179
1476       1 wt_margin(50epi)     chr1 1476  7616618  7617118 3.408917 1.663105e-154 0.4235256
26517      1 wt_margin(50epi)    chr11 5323 25231320 25231820 2.409923 1.208676e-153 0.5648940
255359     1 wt_margin(50epi)     chr8 8029 40335984 40336484 3.779999 1.070482e-136 0.3762477

markerTest1b <- getMarkerFeatures(
  ArchRProj = proj50epi7, 
  useMatrix = "PeakMatrix",
  groupBy = "Clusters_wt_margin",
  testMethod = "wilcoxon",
  bias = c("TSSEnrichment", "log10(nFrags)"),
  useGroups = group2,
  bgdGroups = group1,
  closest = TRUE
)

markerList1b <- getMarkers(markerTest1b, cutOff = "FDR <= 0.05 & abs(Log2FC) >= 0.5", returnGR = FALSE)
list1b <- as.data.frame(markerList1b)
head(list1b) ### total number of significant peaks = 34,741
group         group_name seqnames  idx    start      end    Log2FC           FDR   MeanDiff
201083     1 wt_ectoderm(50epi)     chr4 6087 21660738 21661238 -3.429666 7.619343e-166 -0.5489502
121750     1 wt_ectoderm(50epi)     chr2 7427 35671231 35671731 -3.224900 2.598890e-160 -0.5140247
57405      1 wt_ectoderm(50epi)    chr14 6495 31980882 31981382 -2.671487 1.985932e-158 -0.5754100
1476       1 wt_ectoderm(50epi)     chr1 1476  7616618  7617118 -3.411755 1.608992e-153 -0.4239928
26517      1 wt_ectoderm(50epi)    chr11 5323 25231320 25231820 -2.405494 6.733647e-152 -0.5651749
255359     1 wt_ectoderm(50epi)     chr8 8029 40335984 40336484 -3.797873 1.534297e-136 -0.3770393

I would highly appreciate any feedback on this. Thanks a lot for your help!

[rcorces]

probably just due to random sampling related to the maxCells param.

[mireiacodina]

Thanks a lot for the feedback.
I had to change both maxCells and k parameters to large numbers so I could obtain the same results when changing the order between useGroups and bgdGroups. Changing only maxCells did not work, and I think that the effect of the order when using k = 100 might still be reflected on the selection of different nearby cells for the bias-matched background.

[dhkong78]

I found the parameter "closest" was removed in the latest version. How can I solve the same problem if I would like to get the stable same logFC results?