getPeak2geneLinks corCutoff ABS value

[yamihn]

hello community @rcorces , @jgranja24
I would like to know if it was possible to get a negative correlation between peaks and gene using the getPeaks2GeneLinks function.
As I see from the tutorial:

getPeak2GeneLinks(
  ArchRProj = NULL,
  corCutOff = 0.45,
  FDRCutOff = 1e-04,
  varCutOffATAC = 0.25,
  varCutOffRNA = 0.25,
  resolution = 1,
  returnLoops = TRUE)

corCutOff
A numeric describing the minimum numeric peak-to-gene correlation to return
but it does not specify if is is an ABSOLUTE value of correlation.

So looking in the function getPeak2GeneLinks, my suggestion could be useful:

22 p2g <- metadata(ArchRProj@peakSet)$Peak2GeneLinks
23 p2g <- p2g[which(abs(p2g$Correlation) >= corCutOff & p2g$FDR <= 
24                      FDRCutOff), , drop = FALSE]

Adding the abs() to p2g correlation.
If there is another way to get negatives correlation, I'd appreciate if you can show me of to retain this.

Thank you

[rcorces]

This could be changed to absolute value but I'm more in favor of end users just using the function as written and passing a negative value (-2) as the correlation cutoff. My rationale is that most users are not interested in negative correlations.

[rcorces]

well, an alternative would be to adjust the function to only return positive or negative correlations based on the sign of the corCutOff param. that seems like a better idea than abs value to me

[yamihn]

Yes, It sounds good. Could you apply this modify to the function?

[rcorces]

This has been added to the dev_corCutOff branch via 5942b56

Please re-install and confirm that this works to retrieve both positive and negative correlations for you as I have not tested it.

devtools::install_github("GreenleafLab/ArchR", ref="dev_corCutOff", repos = BiocManager::repositories())
#to unload a package and reload
detach("package:ArchR", unload=TRUE)
library(ArchR)

[yamihn]

Thank you @rcorces , it works perfectly:

p2g <- getPeak2GeneLinks(
  ArchRProj = subset,
  corCutOff = -0.4,
  resolution = 1,
  returnLoops = TRUE)

(p2g[[1]])
GRanges object with 21058 ranges and 2 metadata columns:
          seqnames              ranges strand |     value         FDR
             <Rle>           <IRanges>  <Rle> | <numeric>   <numeric>
      [1]     chr1     4623481-4857694      * | -0.433811 9.56655e-23
      [2]     chr1     9545408-9747920      * | -0.402889 1.91737e-19
      [3]     chr1     9700209-9859290      * | -0.441837 1.16090e-23
      [4]     chr1    9944118-10096714      * | -0.449895 1.31574e-24
      [5]     chr1   10038218-10108157      * | -0.464241 2.34221e-26
      ...      ...                 ...    ... .       ...         ...
  [21054]     chrX 163879589-163958666      * | -0.450336 1.16581e-24
  [21055]     chrX 163879589-164076493      * | -0.604597 7.90067e-49
  [21056]     chrX 166440704-166661991      * | -0.516263 1.83496e-33
  [21057]     chrX 168577233-168773631      * | -0.414319 1.26951e-20
  [21058]     chrX 168795099-169037176      * | -0.461691 4.85671e-26
  -------
  seqinfo: 20 sequences from an unspecified genome; no seqlengths

[rcorces]

Thanks for confirming. This has been merged into dev via #1688 and will be included eventually in release_1.0.3