Subset Peak2Gene Browser Plot

[simonmfr]

Hi, I wanted to subset the peak2gene links shown in the plotBrowserTrack(), to only show those that are linked to the gene given in geneSymbol (not any other genes).

I did the following, which however results in a different plot compared to plotting with loops = getPeak2GeneLinks(archr_proj):

Code

1. Obtain peak2gene file annotated with gene name (described here getting gene name using idxRNA #368)

p2geneDF <- metadata(archr_proj@peakSet)$Peak2GeneLinks
p2geneDF$geneName <- mcols(metadata(p2geneDF)$geneSet)$name[p2geneDF$idxRNA]
p2geneDF$peakName <- (metadata(p2geneDF)$peakSet %>% {paste0(seqnames(.), "_", start(.), "_", end(.))})[p2geneDF$idxATAC]

2. Get all peaks from target gene

library(data.table)
p2geneDF <- as.data.table(p2geneDF)
target_peaks <- p2geneDF[p2geneDF$geneName == "TARGET_GENE",]

3. Add gene_chr and gene_start values based on getGeneAnnotation(archr_proj)$genes

archr_genes <- as.data.table(getGeneAnnotation(archr_proj)$genes)
target_peaks[, symbol := geneName]
gene_info <- archr_genes[target_peaks[, symbol], 
                         .(symbol, gene_chr=seqnames, gene_start=start, gene_end=end, gene_center=(start+end)/2), 
                         on=.(symbol)]
stopifnot(all(target_peaks$geneName == gene_info$symbol))
target_peaks <- cbind(target_peaks, gene_info)

4. Convert table to GRanges for plotting in plotBrowserTrack(loop = gr_target) (loop from peak center to gene start)

# Convert data.table into a GRanges object
convert_to_granges <- function(dt) {

  # Split Region data into Separate Columns
    dt[, c('peak_chr', 'peak_start', 'peak_end') := tstrsplit(peakName, '_', fixed=TRUE)]
    dt[, peak_center := (as.numeric(peak_start) + as.numeric(peak_end)) / 2]

    granges <- GRanges(seqnames = dt$gene_chr, 
                     ranges = IRanges(start = pmin(dt$gene_start, dt$peak_center), 
                                      end = pmax(dt$gene_start, dt$peak_center)),
                                      # end must always be > start, therefore flip if necessary
                     geneName = dt$geneName, 
                     value = dt$Correlation, 
                     FDR = dt$FDR)
  
  return(granges)
}
gr_target <- convert_to_granges(target_peaks)

5. plot

p <- plotBrowserTrack(
    ArchRProj = archr_proj, 
    groupBy = "CellType", 
    geneSymbol = "TARGET_GENE",
    #loops = getPeak2GeneLinks(archr_proj, 
    #                          resolution = 1, 
    #                          returnLoops = TRUE),
    loops = gr_htra1,
    features = GRangesList(ATAC_Peaks = getPeakSet(archr_proj), 
                           TSS = resize(getTSS(archr_proj), width=101, fix='center')
                           )
)

The peak-center positions are largely the same, I just don't get the loop position of the gene when using standardloops = getPeak2GeneLinks(archr_proj). Loops are ending in the middle of nowhere (no TSS or gene start, see image). Any idea why?

Also, is my processing valid?

Thank you!

loops = getPeak2GeneLinks(archr_proj)

loops = gr_target

I also read #1127 and #1313 where you state that addPeakToGeneLinks() uses the start position of the gene from getGenes(ArchRProj). Shouldn't this result in the plot shown below?