Gene labels to peak markers.

[wilfridricher]

Hello,
We were wondering if it is possible to annotate the peaks using the nearest gene to each regulatory element as in the Figure 2c of Satpathy, A.T. et al*.
We have tried using biomaRt and the script is working, but the BiomaRt server is not enough stable to sustain the annotation of many peaks as we got from our dataset.

Thanks a lot in advance
We like ArchR
Jimena and Wilfrid

• Satpathy, A.T., Granja, J.M., Yost, K.E. et al. Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion. Nat Biotechnol 37, 925–936 (2019).

[rcorces]

Thanks for using ArchR. This information should already be present in the GRanges object of your peakSet which can be accessed using getPeakSet(). For example:

> getPeakSet(projHeme5)
GRanges object with 143400 ranges and 12 metadata columns:
        seqnames              ranges strand |     score replicateScoreQuantile groupScoreQuantile Reproducibility
           <Rle>           <IRanges>  <Rle> | <numeric>              <numeric>          <numeric>       <numeric>
   Mono     chr1       752502-753002      * |    15.359                  0.723              0.509               2
      B     chr1       762685-763185      * |   48.8824                  0.933              0.845               2
  CD4.M     chr1       779917-780417      * |   7.02613                  0.335              0.066               2
    GMP     chr1       793283-793783      * |   3.48882                  0.384              0.098               2
      B     chr1       801006-801506      * |   14.6654                  0.697              0.421               2
    ...      ...                 ...    ... .       ...                    ...                ...             ...
   PreB     chrX 154840785-154841285      * |   3.21394                  0.664              0.278               2
  CD4.N     chrX 154841356-154841856      * |   1.35161                  0.282               0.05               2
   PreB     chrX 154842430-154842930      * |   5.32226                  0.866              0.599               2
     NK     chrX 154862053-154862553      * |   12.3754                  0.494                0.2               2
   Mono     chrX 154997960-154998460      * |    8.4037                   0.59              0.308               2
                   GroupReplicate distToGeneStart nearestGene    peakType distToTSS  nearestTSS        GC       idx
                      <character>       <integer> <character> <character> <integer> <character> <numeric> <integer>
   Mono     Mono._.scATAC_PBMC_R1           10149   LINC00115      Distal     10149  uc001aau.3     0.483         1
      B        B._.scATAC_PBMC_R1              32   LINC00115    Promoter        32  uc001aau.3    0.6926         2
  CD4.M    CD4.M._.scATAC_PBMC_R1           17195   LINC01128    Intronic     16988  uc021oeh.1     0.489         3
    GMP GMP._.scATAC_CD34_BMMC_R1           18648      FAM41C      Exonic     18648  uc021oei.1    0.3852         4
      B        B._.scATAC_PBMC_R1           10925      FAM41C      Distal     10925  uc021oei.1    0.4371         5
    ...                       ...             ...         ...         ...       ...         ...       ...       ...
   PreB               PreB._.Rep2            1586       TMLHE    Intronic      1586  uc004cin.3     0.485      3413
  CD4.N    CD4.N._.scATAC_PBMC_R1            1015       TMLHE    Intronic      1015  uc004cin.3    0.3533      3414
   PreB               PreB._.Rep2              57       TMLHE    Promoter        57  uc004cin.3    0.5749      3415
     NK       NK._.scATAC_PBMC_R1           19680       TMLHE      Distal     19680  uc004cin.3    0.4311      3416
   Mono     Mono._.scATAC_PBMC_R1          155587       TMLHE      Distal       758  uc004cin.3    0.4591      3417
  -------
  seqinfo: 23 sequences from an unspecified genome; no seqlengths

[wilfridricher]

Aaah, it's perfect!
I hadn't seen this function! It should solve our problem!
Thank you very much for your quick response and sorry for the inconvenience!

[Austin-s-h]

Hi! I have a follow-up on a similar topic. I did not use the built-in ArchR function to call peaks, so I believe I am missing the helpful annotation it generates. However, I am specifically interested in knowing The set of peaks used to create each genescore. For example, I want to know the locations of Peak1, Peak2, and Peak3 (along with their corresponding weight/strength) and the gene they link to. Would it be possible to extend the functionality of addGeneScoreMatrix() to accomplish this? Are there any functions that may help already? I am interested in something similar to getPeak2GeneLinks() to compare individual coAccessible elements' contributions to GeneScores.

[rcorces]

Probably best to post this type of "follow up" as a separate question since only one answer can be selected per thread.

To answer your question in brief - gene scores arent based on peaks but rather bins. So there is no information on how much a specific peak contributes. You might be able to extract this information by reading the .addGeneScoreMat() function but we dont provide any output at the moment that will help you accomplish what you want. If you think this should be added, feel free to suggest it in a new thread under the Ideas grouping.