Peak-to-gene correlation analysis using an FDR emprically-derived via permutation

[RegnerM2015]

Peak-to-gene correlation analysis with an empirically derived FDR via permutation

Hi ArchR users and @rcorces,

In Regner et al. "A multi-omic single-cell landscape of human gynecologic malignancies", the authors use an empirically-derived FDR for assessing statistical significance of peak-to-gene associations. This is similar to the concept implemented by the addEmpiricalPval parameter in addPeak2GeneLinks. The idea is to find background or null peak-to-gene associations and use this background to estimate the false discovery rate for a chosen alpha. Instead of finding null associations between peaks and genes located on different chromosomes, you can also induce the null condition by permuting the scATAC aggregate/metacell labels (Figure 1). In some cases, this may offer increased statistical power to detect peak-to-gene associations.

Figure 1 | Peak-to-gene correlation analysis with an empirically derived FDR.

Adapted from Regner et al., A multi-omic single-cell landscape of human gynecologic malignancies, Molecular Cell (2021), https://doi.org/10.1016/j.molcel.2021.10.013

A) Schematic illustrating the peak-to-gene correlation analysis. Step 1: the peak count
information from scATAC-seq cells is aggregated into metacells via a KNN algorithm before
computing the correlation between every peak and every gene located up to 250kb away. The
distribution of correlation values and raw p-values are visualized in histograms and the number
of observed peak-to-gene link tests <= alpha is recorded. Step 2: the peak-to-gene correlations
are re-computed for a permuted null version of the dataset where the scATAC-seq metacell
labels are shuffled, breaking the link or any potential correlation between peaks and genes. The
process is repeated for a total of 100 permutations. For each permutation, the number of null
peak-to-gene link tests <= alpha is recorded. As an example, the distribution of null correlation
values and null raw p-values are visualized in histograms for one permutation run. Step 3: The
median number of null peak-to-gene link tests <= alpha across 100 permutations is divided by
the number of observed peak-to-gene link tests <= alpha to arrive at an empirically derived FDR
(eFDR).
B) The peak-to-gene correlation analysis, as described in A, but using an alpha threshold equal
to the first quartile in the distribution of raw p-values.
C) The peak-to-gene correlation analysis, as described in A, but using an alpha threshold equal
to 1e-12.
D) Flow chart demonstrating the screening procedure used to identify distal peak-to-gene links
with positive regulatory effects. The pie chart depicts the proportion of positive regulatory peak-
to-gene links that are distal peak-to-gene relationships.

How to implement using ArchR?

The authors developed this approach by making only a slight modification the source code of addPeak2GeneLinks to create a new function addPermPeak2GeneLinks. Before computing the correlations, the ATAC aggregate/metacell labels are randomly shuffled with sample() while the corresponding RNA aggregate/metacell labels are left unmodified. This effectively breaks the link between peak accessibility and gene expression offering a robust set of null peak-to-gene associations:

Code from https://github.com/RegnerM2015/scENDO_scOVAR_2020/blob/main/PeaktoGeneLink_Analysis/Full_Cohort/Archr_Peak_Null_Permute.R#L222-L241

  # #Null Correlations
  # .logDiffTime(main="Computing Background Correlations", t1=tstart, verbose=verbose, logFile=logFile)
  # nullCor <- .getNullCorrelations(seATAC, seRNA, o, 1000)
  # saveRDS(nullCor,"nullCor.rds")
  # .logDiffTime(main="Computing Correlations", t1=tstart, verbose=verbose, logFile=logFile)
  # 
  rm(.Random.seed, envir=globalenv())
  # n <- floor(runif(1, min=0, max=1001))
  # set.seed(n)
  seATAC.mat <- assay(seATAC)
  seATAC.mat <- seATAC.mat[,sample(ncol(seATAC.mat))]#Permute ATAC aggregates
  
  o$Correlation <- rowCorCpp(as.integer(o$A), as.integer(o$B), seATAC.mat,assay(seRNA))
  o$VarAssayA <- .getQuantiles(matrixStats::rowVars(seATAC.mat))[o$A]
  o$VarAssayB <- .getQuantiles(matrixStats::rowVars(assay(seRNA)))[o$B]
  o$TStat <- (o$Correlation / sqrt((1-o$Correlation^2)/(ncol(seATAC.mat)-2))) #T-statistic P-value
  o$Pval <- 2*pt(-abs(o$TStat), ncol(seATAC.mat) - 2)
  o$FDR <- p.adjust(o$Pval, method = "fdr")
  # o$EmpPval <- 2*pnorm(-abs(((o$Correlation - mean(nullCor[[2]])) / sd(nullCor[[2]]))))
  # o$EmpFDR <- p.adjust(o$EmpPval, method = "fdr")

Limitations

There are some limitations associated with this approach. First, this approach is computationally expensive as you would have to run the addPermPeak2GeneLinks n times to generate n permutations. Secondly, it may be difficult to implement this framework into the current ArchR suite of tools. We are open to suggestions and new ideas for making this approach accessible to the community!

[rcorces]

better late than never? This has now been incorporated into dev via e3c6d4c
Will become part of release_1.0.3. Thanks for contributing to ArchR!

[RegnerM2015]

Thank you for implementing our approach into the ArchR suite! We appreciate your time and efforts making this work.

[RegnerM2015]

Hi @rcorces,

After looking at this implementation some more, I wanted to ask if it would be reasonable to calculate the empirical p-values similar to .getNullCorrelations (e3c6d4c#diff-8b250bfa8bdda4cf87f0737474b753677e0c74ef4ec24c617e5f6e886b63f984R1286).

For every permutation, or shuffling of metacells, we would store the null correlation for every P2G. After N permutations, we would would use the mean and sd of this null correlation set (for every P2G) to calculate p-values:

out$EmpPval <- 2*pnorm(-abs(((out$Correlation - mean(nullCorFromShuffling[[2]])) / sd(nullCorFromShuffling[[2]]))))
out$EmpFDR <- p.adjust(out$EmpPval, method = "fdr")

Then correct the EmpPvals with the BH method ("fdr").

Does this approach seem reasonable? In some ways, shuffling labels could provide a better set of null correlations relative to peaks and genes not on the same chromosome, especially in the context of disease states where chromosomal translocations can occur.

[RegnerM2015]

I guess .getNullCorrelations is used to estimate empirical pvalues where as the approach we introduced is a method to estimate empirical FDR.

Looks like I may have answered my own question.

[AnjaliC4]

Hi @rcorces , how can I use this approach currently with ArchR?

[AnjaliC4]

Hi @rcorces, I tried using this approcah but ended up with this error:

devtools::install_github("GreenleafLab/ArchR", ref="dev", repos = BiocManager::repositories())
proj <- addPeak2GeneLinks(ArchRProj = proj,maxDist = 500000,addPermutedPval = TRUE,nperm = 100)
ArchR-addPeak2GeneLinks-4efe4aee7bee-Date-2022-12-04_Time-00-12-06.log

ArchR logging to : ArchRLogs/ArchR-addPeak2GeneLinks-4efe4aee7bee-Date-2022-12-04_Time-00-12-06.log
If there is an issue, please report to github with logFile!
2022-12-04 00:12:06 : Getting Available Matrices, 0.005 mins elapsed.
2022-12-04 00:16:18 : Filtered Low Prediction Score Cells (29695 of 201735, 0.147), 0.022 mins elapsed.
2022-12-04 00:16:25 : Computing KNN, 0.14 mins elapsed.
2022-12-04 00:16:30 : Identifying Non-Overlapping KNN pairs, 0.218 mins elapsed.
2022-12-04 00:16:33 : Identified 500 Groupings!, 0.262 mins elapsed.
2022-12-04 00:16:33 : Getting Group RNA Matrix, 0.264 mins elapsed.
2022-12-04 00:28:03 : Getting Group ATAC Matrix, 11.761 mins elapsed.
2022-12-04 00:42:17 : Normalizing Group Matrices, 26.003 mins elapsed.
2022-12-04 00:42:32 : Finding Peak Gene Pairings, 26.25 mins elapsed.
2022-12-04 00:42:34 : Computing Correlations, 26.285 mins elapsed.
Performing Permuted P-values similar to Regner et al., 2021
Running Permutation 1 of 100
Error in [[<-(*tmp*, name, value = numeric(0)) :
0 elements in value to replace 5439471 elements
Calls: addPeak2GeneLinks -> $<- -> [[<- -> [[<-
Execution halted

[rcorces]

@AnjaliC4 - I will eventually get to debugging this but have not had time yet

[jblich870]

@rcorces - This is a great addition in principle and the new code runs fine for me.

Where is the new PermFDR stored? I can't seem to find it in the output from getPeak2GeneLinks. Code below (nperm=2 chosen just for example):

Multiproj3_sub_big3 <- addPeak2GeneLinks(
  ArchRProj = Multiproj3_sub_big3,
  reducedDims = "Harmony_ATAC",
  useMatrix = "GeneExpressionMatrix",
  addEmpiricalPval = T,
  addPermutedPval = T,
  nperm = 2,
  maxDist = 500000)

p2g.links <- getPeak2GeneLinks(
  ArchRProj = Multiproj3_sub_big3,
  corCutOff = 0.25,
  FDRCutOff = 1,
  varCutOffATAC = 0.5,
  varCutOffRNA = 0.5,
  returnLoops = F)

p2g <- as.data.frame(p2g.links)

gives

ArchR logging to : ArchRLogs/ArchR-addPeak2GeneLinks-30346c0247d-Date-2023-01-04_Time-17-03-01.log
If there is an issue, please report to github with logFile!
2023-01-04 17:03:01 : Getting Available Matrices, 0.001 mins elapsed.
2023-01-04 17:03:02 : Filtered Low Prediction Score Cells (0 of 9147, 0), 0.003 mins elapsed.
2023-01-04 17:03:02 : Computing KNN, 0.008 mins elapsed.
2023-01-04 17:03:02 : Identifying Non-Overlapping KNN pairs, 0.01 mins elapsed.
2023-01-04 17:03:04 : Identified 497 Groupings!, 0.042 mins elapsed.
2023-01-04 17:03:04 : Getting Group RNA Matrix, 0.043 mins elapsed.
2023-01-04 17:03:50 : Getting Group ATAC Matrix, 0.809 mins elapsed.
2023-01-04 17:05:07 : Normalizing Group Matrices, 2.096 mins elapsed.
2023-01-04 17:05:12 : Finding Peak Gene Pairings, 2.179 mins elapsed.
2023-01-04 17:05:14 : Computing Correlations, 2.201 mins elapsed.
2023-01-04 17:05:20 : Computing Background Correlations, 2.304 mins elapsed.
Performing Permuted P-values similar to Regner et al., 2021
Running Permutation 1 of 2
Running Permutation 2 of 2
2023-01-04 17:06:05 : Completed Peak2Gene Correlations!, 3.052 mins elapsed.
ArchR logging successful to : ArchRLogs/ArchR-addPeak2GeneLinks-30346c0247d-Date-2023-01-04_Time-17-03-01.log

> head(p2g)
  idxATAC idxRNA Correlation          FDR  VarQATAC   VarQRNA    EmpPval    EmpFDR
1      34     15   0.3699423 3.661258e-16 0.5477846 0.7584939 0.02054640 0.5085431
2     117     15   0.2608162 3.363955e-08 0.6911543 0.7584939 0.10367816 0.9522413
3     121     15   0.2874845 7.769521e-10 0.8503387 0.7584939 0.07255156 0.8462945
4     203     15   0.3981960 8.123123e-19 0.9547550 0.7584939 0.01261027 0.4064423
5     214     15   0.2997554 1.183681e-10 0.8396979 0.7584939 0.06105312 0.7954710
6      47     18   0.2842973 1.247085e-09 0.9919732 0.6733081 0.07581211 0.8591338

[rcorces]

@jblich870 - there is some active development going on here that isnt quite finished. But you would currently need to be on the dev_p2g branch to get the desired output.

[jblich870]

@rcorces - great, thanks - I'll give this a go

[AnjaliC4]

Hi @RegnerM2015 - thanks for this method. I was wondering if you get the above error or not?
In addition, I was wondering two things,

1. Can this method be also applied to peak co-accessibility with the same logic?
2. When using background correlations against true correlations, how do you determine the optimal correlation threshold?
   Thanks :)

[RegnerM2015]

Hi @AnjaliC4,

I have not yet tested out these new parameters on the dev branch.

Regarding your questions:

1. Yes, in theory, one could could apply a similar permutation strategy to estimate the eFDR for peak co-accessibility analysis.
2. If you have a background set of null correlations (meaning that there are no real peak-to-gene associations), you can plot a histogram of the correlation values and summarize with summary in R. This should tell you what the extreme correlation values are in this null condition and give you an idea of an appropriate correlation threshold. For example, if the highest correlation in the null condition is 0.25, you may reason that any correlation in the real data less than 0.25 is probably not a robust peak-to-gene association.

[AnjaliC4]

Hi @RegnerM2015 , that makes a lot of sense. Thank you :)

[AnjaliC4]

Hi again @RegnerM2015 , do you suggest using a correlation cut-off as well after using eFDR cut-off or is the idea is that eveything passing eFDR regardless of correlation could be robust peak-to-gene linkages? It may not be easy to summarise null correlations from the pipeline as there will be 100 null correlations for every gene.

[RegnerM2015]

Depending on your situation and downstream analysis tasks, you may consider uing a correlation cutoff as well after applying the significance threshold. While everything passing the significance threshold (eFDR) could be robust peak-to-gene linkages regardless of correlation, it is still important to check the distribution of correlation values after applying this significance threshold in my opinion. Hope this helps!