Append peakType column to P2G output

[RegnerM2015]

Hi @rcorces and @jgranja24 ,

I want to append the peakType to the output of Peak2Gene links. Is this okay to use? Please let me know if there is a more efficient way of doing this. Thanks in advance!

> p2geneDF <- metadata(proj@peakSet)$Peak2GeneLinks
> p2geneDF$geneName <- mcols(metadata(p2geneDF)$geneSet)$name[p2geneDF$idxRNA]
> p2geneDF$peakName <- (metadata(p2geneDF)$peakSet %>% {paste0(seqnames(.), ":", start(.), "-", end(.))})[p2geneDF$idxATAC]
> annot <- readRDS(metadata(p2geneDF)$seATAC)
> p2geneDF$peakType <- annot@rowRanges$peakType[p2geneDF$idxATAC]
> p2g.df.obs <- as.data.frame(p2geneDF)
> p2g.df.obs <- p2g.df.obs[complete.cases(p2g.df.obs),]
> head(p2g.df.obs)
  idxATAC idxRNA Correlation      RawPVal          FDR  VarQATAC   VarQRNA   EmpPval    EmpFDR geneName           peakName peakType
1       1      2   0.1196724 7.567003e-03 1.482268e-02 0.3052280 0.1450149 0.5659268 0.9999790   FAM87B chr1:794813-795313   Distal
2       2      2   0.3420600 4.351283e-15 3.609684e-14 0.5430366 0.1450149 0.0836577 0.6683863   FAM87B chr1:804672-805172   Distal
3       3      2  -0.1826303 4.208727e-05 1.235637e-04 0.7824071 0.1450149 0.3187364 0.9955482   FAM87B chr1:817102-817602 Promoter
4       4      2  -0.1150224 1.027849e-02 1.954973e-02 0.3068608 0.1450149 0.5184814 0.9999790   FAM87B chr1:817782-818282   Exonic
5       5      2  -0.1957482 1.105089e-05 3.535143e-05 0.3405486 0.1450149 0.2867701 0.9775311   FAM87B chr1:818808-819308   Exonic
6       6      2   0.1795370 5.693157e-05 1.638087e-04 0.5402933 0.1450149 0.3760558 0.9999790   FAM87B chr1:826658-827158   Exonic

Note that annot@rowRanges is the same set of ranges as metadata(p2geneDF)$peakSet:

> test <- metadata(p2geneDF)$peakSet
> length(test)
[1] 409733
> test <- paste0(test@seqnames,"-",test@ranges)
> length(annot)
[1] 409733
> length(which(test == paste0(annot@rowRanges@seqnames,"-",annot@rowRanges@ranges)))
[1] 409733

[rcorces]

seems fine to me

[RegnerM2015]

Thanks for the quick reply!

[RegnerM2015]

Hi @rcorces and @jgranja24 ,

I am trying to perform a similar operation for ArchR's CoAccessibility implementation.

After running addCoAccessibility on the ArchR project, I use the following code below to extract the corresponding peakNames based on the indices provided byqueryHits and subjectHits:

coaDF <- metadata(proj@peakSet)$CoAccessibility
coaDF$queryPeak <- (metadata(coaDF)$peakSet %>% {paste0(seqnames(.), ":", start(.), "-", end(.))})[coaDF$queryHits]
coaDF$subjectPeak <- (metadata(coaDF)$peakSet %>% {paste0(seqnames(.), ":", start(.), "-", end(.))})[coaDF$subjectHits]
coaDF <- as.data.frame(coaDF)

head(coaDF)
  queryHits subjectHits seqnames  correlation Variability1 Variability2       TStat         Pval          FDR VarQuantile1 VarQuantile2
1         1           2     chr1 -0.400501677 0.0008446586 0.0027217335 -9.75401023 1.096224e-20 1.652183e-19    0.4542218    0.7424770
2         1           3     chr1 -0.060927063 0.0008446586 0.0004418009 -1.36217369 1.737589e-01 2.488891e-01    0.4542218    0.2711830
3         1           4     chr1  0.039121375 0.0008446586 0.0047801634  0.87369808 3.827038e-01 4.774914e-01    0.4542218    0.8474390
4         1           5     chr1  0.494381301 0.0008446586 0.0017332506 12.69212766 3.542041e-32 1.042517e-30    0.4542218    0.6403246
5         2           1     chr1 -0.400501677 0.0027217335 0.0008446586 -9.75401023 1.096224e-20 1.652183e-19    0.7424759    0.4542225
6         2           3     chr1  0.001323007 0.0027217335 0.0004418009  0.02952414 9.764584e-01 9.831759e-01    0.7424759    0.2711830
           queryPeak        subjectPeak
1 chr1:804681-805181 chr1:817112-817612
2 chr1:804681-805181 chr1:826480-826980
3 chr1:804681-805181 chr1:827270-827770
4 chr1:804681-805181 chr1:832050-832550
5 chr1:817112-817612 chr1:804681-805181
6 chr1:817112-817612 chr1:826480-826980

Please let me know if this a valid approach for linking up the peakNames to the indices. Note that this approach was inspired by #368 (comment). Thanks!

[rcorces]

I dont want to get into the business of checking other people's code. You'll have to do your own spot checks to make sure its correct.

[RegnerM2015]

Sounds good! Thanks for your response.