Ordering of seqnames in objects returned by getPeak2GeneLinks()

[jmmuncie]

Hi all - I've noticed that when I run getPeak2GeneLinks() with returnLoops = TRUE, the GRanges object that is returned is organized by seqnames in proper numerical order (chr1, chr2, chr3...chr9, chr10, chr11..chr19, chrX). However, when I run getPeak2GeneLinks() with returnLoops = FALSE, the dataframe that is returned seems to be ordered by seqnames as follows: chr1, chr11, chr12, chr13...chr19, chr2, chr3, chr4..chr9, chrX.

Is there a way to get this dataframe to return with the seqnames in proper numerical order? This is causing issues for me downstream as I'm attempting to plot peak2gene links on browser tracks, and to only display lops that link to differentially accessible peaks. If the objects were returned in the same order, it would be fairly straightfoward, but as is, I don't have a way to properly match between the two.

[rcorces]

depending on exactly what you mean, I think this is actually straightforward. Just sort by idxATAC

p2g <- getPeak2GeneLinks(ArchRProj = projHeme5, returnLoops = FALSE)
p2g <- p2g[order(p2g$idxATAC),]

That being said the results with and without returnLoops = TRUE are going to be different. There isnt a one-to-one match, even if resolution = NULL.

[jmmuncie]

Hi - sorry, I wasn't clear about what I'm attempting to do in my initial post. I would like to plot a browser track that shows only differentially accessible peaks (DAPs) between two groups using the features input, and also show only the loops that link to those specific DAPs (via the loops input), rather than showing all the p2g loops found in the dataset.

I was attempting to do this by first finding matches between my GRanges object of DAPs obtained by pairwise testing between two groups and the GRanges object of linked peaks to genes:

#Find differential accessible peaks between two groups, e.g. KO vs. WT CMs
markers <- getMarkerFeatures(
  ArchRProj = proj, 
  useMatrix = "PeakMatrix",
  groupBy = "ClusterByGenotype",
  testMethod = "wilcoxon",
  bias = c("TSSEnrichment", "log10(nFrags)"),
  useGroups = "KO_CMs",
  bgdGroups = "WT_CMs"
)
markers_upinWT <- getMarkers(markers, cutOff = "FDR <= 0.1 & Log2FC <= -0.5", returnGR = TRUE)

#Convert this to a GRanges object
markers_upinWT_gr <- markers_upinWT$`KO_CMs`

#Get peak to gene linkages
proj <- addPeak2GeneLinks(
  ArchRProj = proj,
  reducedDims = "Harmony_Combined_sub",
  useMatrix = "GeneExpressionMatrix"
  )

#Retrieve p2g links as dataframe
p2g_df <- getPeak2GeneLinks(
  ArchRProj = proj,
  corCutOff = 0.45,
  resolution = 1,
  returnLoops = FALSE
)

#Get peak ids from p2g_df
peakids <- p2g_df$idxATAC

#Get the GR object that consists of all peaks 
gr_peaks <- metadata(p2g_df)[[1]]

#Subset only the peaks that are linked to genes
linked_peaks <- gr_peaks[peakids]

#Now find matches between our differential peaks up in WT
#and the linked peaks:
match <- findOverlaps(markers_upinWT_gr, linked_peaks)

#Finally, we can pull the differential peaks that are linked to genes from the 
#full linked peaks data frame
p2g_diff_df <- p2g_df[match@to, ]

Then, I was thinking that if the p2g loops GRanges object was in the same order as in the GRanges object associated with the p2g dataframe, I'd be able to use those same indexes ( [match@to, ] ) to pull only the linked loops and feed those into the browser track plot, something like this:

p2g_loops <- getPeak2GeneLinks(
  ArchRProj = proj_Mef2c_v13_E9_subset,
  corCutOff = 0.45,
  resolution = 1,
  returnLoops = TRUE
)
p2g_loops_gr <- p2g_loops$Peak2GeneLinks

p2g_diff_loops <- p2g_loops_gr[match@to, ]

tracks <- plotBrowserTrack(
  ArchRProj = proj, 
  groupBy = "ClusterByGenotype", 
  useGroups = c("KO_CMs", "WT_CMs"),
  geneSymbol = c("Nppa"),
  upstream = 100000,
  downstream = 100000,
  features =  getMarkers(markers, cutOff = "FDR <= 0.1 & Log2FC <= -0.5", returnGR = TRUE)["KO_CMs"],
  loops = p2g_diff_loops
)

However, as I was attempting to do that, I realized the problem that I initially posted - which is that the the GRanges object returned in p2g_loops is organized by seqnames in proper numerical order (chr1, chr2, chr3...chr9, chr10, chr11..chr19, chrX), whereas the GRanges object associated with the p2g_df dataframe is ordered by seqnames as follows: chr1, chr11, chr12, chr13...chr19, chr2, chr3, chr4..chr9, chrX. As a result, the [match@to, ] indexes couldn't be used to pull out the relevant loops.

I didn't realize until I read your response that there isn't necessarily a 1:1 match between the peak2gene loops and the p2g dataframe that matches idxATAC and idxRNA. If that's the case, maybe there isn't really a way to accomplish this.

Hopefully this clarifies a bit what I'm attempting to do. Please let me know if you have any thoughts about how I might be able to subset only the loops that are linked to differential peaks to feed into the browser track plot.

[rcorces]

I'm not sure if I'm fully understanding but I would use idxATAC to extract the peak regions as a GRanges object and then I would forget about indexing and just use overlap (type = "equal") to subset your objects. Something like:

projHeme5 <- addPeak2GeneLinks(
  ArchRProj = projHeme5,
  reducedDims = "IterativeLSI"
)

p2g <- getPeak2GeneLinks(
  ArchRProj = projHeme5,
  corCutOff = 0.45,
  resolution = 1,
  returnLoops = FALSE
)

p2g_gr <- sort(sortSeqlevels(metadata(p2g)$peakSet[p2g$idxATAC]))

peak_gr <- getPeakSet(ArchRProj = projHeme5)

#randomly simulate a subset of the peakset as if they were differential peaks
diff_gr <- peak_gr[sample(x = 1:length(peak_gr), size = 1000, replace = FALSE)]

p2g_subset_gr <- subsetByOverlaps(p2g_gr, diff_gr, type = "equal")

[jmmuncie]

I think what you've suggested is a good alternative for how to subset the overlapping peaks between the differential peak set and the full peak set.

The tougher challenge I've run into is that I'd like to then plot a browser track that shows Peak2GeneLinks under the browser track, but only for links that are anchored to the differential peaks. This is challenging because the loops input for plotBrowserTrack requires the GRanges output from getPeak2GeneLinks with returnLoops = T. So far I haven't been able to come up with a way to match the set of differential peaks to the loops, such that I can then subset the loops object. Please let me know if you can think of a way to do this.

For instance, if the GRanges output from getPeak2GeneLinks with returnLoops = T also included the idxATAC of the peak being linked in each loop, it would be easy to match to my list of differential peaks, and then I'd be able to subset the loops object accordingly before passing it to plotBrowserTrack. But maybe I'm not understanding the way the returnLoops=T output is generated and it's not actually possible to include the associated idxATAC.

[rcorces]

Sorry for not being more explicit in my suggestion. I would take the p2g dataframe, create a GRanges object from the idxATAC position and a GRanges object of the promoter from idxRNA and find overlaps with those. If you keep track of the idxATAC and idxRNA in those overlaps, you can recreate a loops object from the links that overlap. This feels like it should work for what you are trying to accomplish? The loops object is really only designed for plotting. And depending on the resolution there just isnt a one-to-one mapping of idxATAC to idxRNA in the context of loops.

You can see the code to generate the loops from p2g here:

ArchR/R/IntegrativeAnalysis.R

Lines 1320 to 1354 in f6c0388

	if(returnLoops){
	peakSummits <- resize(metadata(p2g)$peakSet, 1, "center")
	geneStarts <- resize(metadata(p2g)$geneSet, 1, "start")
	if(!is.null(resolution)){
	summitTiles <- floor(start(peakSummits) / resolution) * resolution + floor(resolution / 2)
	geneTiles <- floor(start(geneStarts) / resolution) * resolution + floor(resolution / 2)
	}else{
	summitTiles <- start(peakSummits)
	geneTiles <- start(geneStarts)
	}
	loops <- .constructGR(
	seqnames = seqnames(peakSummits[p2g$idxATAC]),
	start = summitTiles[p2g$idxATAC],
	end = geneTiles[p2g$idxRNA]
	)
	mcols(loops)$value <- p2g$Correlation
	mcols(loops)$FDR <- p2g$FDR
	loops <- loops[order(mcols(loops)$value, decreasing=TRUE)]
	loops <- unique(loops)
	loops <- loops[width(loops) > 0]
	loops <- sort(sortSeqlevels(loops))
	loops <- SimpleList(Peak2GeneLinks = loops)
	return(loops)
	}else{
	return(p2g)
	}

hope that helps?

[jmmuncie]

Thanks, that does help a lot! I'll give it a try and let you know how it goes.