How is "nearest gene" determined?

[mitchTwoTimes]

Hi,

When you run getMarkerFeatures, there is a metadata column deeply embedded in the results called nearestGene.
markerTest@metadata$Params$ArchRProj@peakSet$nearestGene

How is this calculated? Most of the time when I plot the peak region, the gene in the browser track is not the nearestGene. For example, the nearestGene according to the getMarkerFeatures output is DYNC1H1, but that's not the case here or in UCSC Genome Browser.

[rcorces]

I think you've gone a little too far into searching for something. What you're looking at is the metadata of the object and specifically the params that were used when calling getMarkerFeatures(). One of those params was the ArchRProject object. By looking at markerTest@metadata$Params$ArchRProj@peakSet$nearestGene, you are effectively just doing
ArchRProj@peakSet$nearestGene which is actually disconnected from the analysis.

> all.equal(ArchRProj@peakSet$nearestGene, markersPeaks@metadata$Params$ArchRProj@peakSet$nearestGene)
[1] TRUE

The nearest gene is not stored in the object explicitly. You would have to match the values in the idx column of the SummarizedExperiment with the values in the idx column of the peakSet but keep in mind that those index values are only unique per chromosome, not across the whole peakSet.

[mitchTwoTimes]

Maybe I should back up and make sure I correctly understand the SummarizedExperiment generated from getMarkerFeatures. Here is an example of a getMarkerFeatures result:

> markerTest
class: SummarizedExperiment 
dim: 396696 1 
metadata(2): MatchInfo Params
assays(7): Log2FC Mean ... AUC MeanBGD
rownames(396696): 1 2 ... 396695 396696
rowData names(4): seqnames idx start end
colnames(1): ulEN_Cortex_GW17
colData names(0):

These are the results from my differential accessibility (DA) comparison, right? Tried to model after this answer.

> DApeaks <- markerTest@elementMetadata
> DApeaks$name <- paste0(DApeaks$seqnames, ":", DApeaks$start, "-", DApeaks$end)
> DApeaks$lfc <- markerTest@assays@data$Log2FC$x
> DApeaks$FDR <- markerTest@assays@data$FDR$x
> DApeaks
DataFrame with 396696 rows and 7 columns
       seqnames     idx     start       end                   name       lfc       FDR
          <Rle> <array>   <array>   <array>            <character> <numeric> <numeric>
1          chr1       1    826629    827129     chr1:826629-827129  1.508038  0.763488
2          chr1       2    827310    827810     chr1:827310-827810  1.585744  0.429330
3          chr1       3    854246    854746     chr1:854246-854746 -0.563896  0.447288
4          chr1       4    859455    859955     chr1:859455-859955  0.000000        NA
5          chr1       5    869696    870196     chr1:869696-870196  1.638633  0.802521
...         ...     ...       ...       ...                    ...       ...       ...
396692     chrX    7671 155881006 155881506 chrX:155881006-15588..  0.826104  0.447288
396693     chrX    7672 155891852 155892352 chrX:155891852-15589..  2.527003  0.429330
396694     chrX    7673 155905377 155905877 chrX:155905377-15590..  0.000000        NA
396695     chrX    7674 155941855 155942355 chrX:155941855-15594..  0.507689  0.447288
396696     chrX    7675 155956624 155957124 chrX:155956624-15595..  0.882480  0.447288

If that's all true, then would this be an appropriate way to get around the per-chromosome idx?

> peakSet <- proj@peakSet
> peakSet$name <- paste0(peakSet@seqnames, ":", start(peakSet@ranges), "-", end(peakSet@ranges))
> all.equal(length(peakSet$name), length(DApeaks$name))
[1] TRUE
> stopifnot(peakSet$name %in% DApeaks$name)
> all.equal(peakSet$name, DApeaks$name)
[1] "352878 string mismatches"

> DApeaks.ordered <- DApeaks[match(peakSet@elementMetadata$name, DApeaks$name), ]
> all.equal(peakSet$name, DApeaks.ordered$name)
[1] TRUE

Please let me know if I misunderstood something @rcorces

[rcorces]

That seems perfectly fine to me.

[mitchTwoTimes]

Thanks for the validation! Can we please circle back to my original question: How is nearest gene determined? Is it just the literal closest annotated gene to the open peak?

[rcorces]

yes. the nearest gene present in your geneAnnotation