LinkPeaks() works but Links don't show in CoveragePlot() using multiome data

[RaghadShu]

Hi,
Thank you for the great tool. I am working with 5 sc-multiome samples of pathogenic plasma cells. I integrated the RNA with rpca, the ATAC with rlsi (by first creating a peak union), then both modalities with WNN.
I called the peaks on the integrated ATAC using group.by the wnn clustering resolution. I now have an assay called "peakunion_MACS2". I am trying to link the peaks to genes to see which peaks regulate the genes.

This is what I did to link peaks: I started with "XBP1" and "SDC1" which are known markers for this data

Idents(seurat) <- "wsnn_res.0.3"
DefaultAssay(seurat) <- "peakunion_MACS2"
# first compute the GC content for each peak
seurat <- RegionStats(seurat, genome = BSgenome.Hsapiens.UCSC.hg38)
# link peaks to genes
seurat <- LinkPeaks(
  object = seurat,
  peak.assay = "peakunion_MACS2",
  expression.assay = "RNA",
  genes.use = c("XBP1", "SDC1"))

and when I do Links(seurat) I get this information:

GRanges object with 2 ranges and 5 metadata columns:
      seqnames            ranges strand |     score        gene
         <Rle>         <IRanges>  <Rle> | <numeric> <character>
  [1]     chr2 20225433-20579145      * | 0.0554166        SDC1
  [2]    chr22 28800597-28823298      * | 0.0522171        XBP1
                        peak    zscore      pvalue
                 <character> <numeric>   <numeric>
  [1] chr2-20578650-20579640   3.73339 9.44593e-05
  [2] chr22-28822785-28823..   3.96908 3.60757e-05
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

But when I plot using this:

CoveragePlot(
  object = seurat,
  region = "SDC1",
  features = "SDC1",
  expression.assay = "RNA",
  extend.upstream = 5000,
  extend.downstream = 5000,
  links = TRUE
)

I get this plot, which has no links.

Can you help by pointing out if I am missing anything?

Best,
Raghad

> sessionInfo()
R version 4.2.2 Patched (2022-11-10 r83330)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.6 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_ES.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=es_ES.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=es_ES.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] BSgenome.Hsapiens.UCSC.hg38_1.4.5 BSgenome_1.66.3                  
 [3] rtracklayer_1.58.0                Biostrings_2.66.0                
 [5] XVector_0.38.0                    EnsDb.Hsapiens.v86_2.99.0        
 [7] ensembldb_2.22.0                  AnnotationFilter_1.22.0          
 [9] GenomicFeatures_1.50.4            AnnotationDbi_1.60.2             
[11] Biobase_2.58.0                    GenomicRanges_1.50.2             
[13] GenomeInfoDb_1.34.9               IRanges_2.32.0                   
[15] S4Vectors_0.36.2                  BiocGenerics_0.44.0              
[17] ggplot2_3.4.4                     Signac_1.12.9003                 
[19] Seurat_5.0.1.9001                 SeuratObject_5.0.1               
[21] sp_2.1-2                         

loaded via a namespace (and not attached):
  [1] utf8_1.2.4                  spatstat.explore_3.2-5     
  [3] reticulate_1.34.0           tidyselect_1.2.0           
  [5] RSQLite_2.3.4               htmlwidgets_1.6.4          
  [7] grid_4.2.2                  BiocParallel_1.32.6        
  [9] Rtsne_0.17                  munsell_0.5.0              
 [11] codetools_0.2-19            ica_1.0-3                  
 [13] future_1.33.1               miniUI_0.1.1.1             
 [15] withr_3.0.0                 spatstat.random_3.2-2      
 [17] colorspace_2.1-0            progressr_0.14.0           
 [19] filelock_1.0.3              ROCR_1.0-11                
 [21] tensor_1.5                  listenv_0.9.0              
 [23] MatrixGenerics_1.10.0       labeling_0.4.3             
 [25] GenomeInfoDbData_1.2.9      polyclip_1.10-6            
 [27] bit64_4.0.5                 farver_2.1.1               
 [29] parallelly_1.36.0           vctrs_0.6.5                
 [31] generics_0.1.3              BiocFileCache_2.6.1        
 [33] R6_2.5.1                    bitops_1.0-7               
 [35] spatstat.utils_3.0-4        cachem_1.0.8               
 [37] DelayedArray_0.24.0         promises_1.2.1             
 [39] BiocIO_1.8.0                scales_1.3.0               
 [41] gtable_0.3.4                globals_0.16.2             
 [43] goftest_1.2-3               spam_2.10-0                
 [45] rlang_1.1.3                 RcppRoll_0.3.0             
 [47] splines_4.2.2               lazyeval_0.2.2             
 [49] spatstat.geom_3.2-7         yaml_2.3.8                 
 [51] reshape2_1.4.4              abind_1.4-5                
 [53] httpuv_1.6.13               tools_4.2.2                
 [55] ellipsis_0.3.2              RColorBrewer_1.1-3         
 [57] ggridges_0.5.5              Rcpp_1.0.12                
 [59] plyr_1.8.9                  progress_1.2.3             
 [61] zlibbioc_1.44.0             purrr_1.0.2                
 [63] RCurl_1.98-1.14             prettyunits_1.2.0          
 [65] deldir_2.0-2                pbapply_1.7-2              
 [67] cowplot_1.1.2               zoo_1.8-12                 
 [69] SummarizedExperiment_1.28.0 ggrepel_0.9.5              
 [71] cluster_2.1.6               magrittr_2.0.3             
 [73] data.table_1.14.10          RSpectra_0.16-1            
 [75] scattermore_1.2             lmtest_0.9-40              
 [77] RANN_2.6.1                  ProtGenerics_1.30.0        
 [79] fitdistrplus_1.1-11         matrixStats_1.2.0          
 [81] hms_1.1.3                   patchwork_1.2.0            
 [83] mime_0.12                   xtable_1.8-4               
 [85] XML_3.99-0.16               fastDummies_1.7.3          
 [87] gridExtra_2.3               compiler_4.2.2             
 [89] biomaRt_2.54.1              tibble_3.2.1               
 [91] KernSmooth_2.23-22          crayon_1.5.2               
 [93] htmltools_0.5.7             later_1.3.2                
 [95] tidyr_1.3.0                 DBI_1.2.1                  
 [97] dbplyr_2.4.0                MASS_7.3-60.0.1            
 [99] rappdirs_0.3.3              Matrix_1.6-5               
[101] cli_3.6.2                   parallel_4.2.2             
[103] dotCall64_1.1-1             igraph_1.6.0               
[105] pkgconfig_2.0.3             GenomicAlignments_1.34.1   
[107] plotly_4.10.4               spatstat.sparse_3.0-3      
[109] xml2_1.3.6                  stringr_1.5.1              
[111] digest_0.6.34               sctransform_0.4.1          
[113] RcppAnnoy_0.0.21            spatstat.data_3.0-4        
[115] leiden_0.4.3.1              fastmatch_1.1-4            
[117] uwot_0.1.16                 restfulr_0.0.15            
[119] curl_5.2.0                  shiny_1.8.0                
[121] Rsamtools_2.14.0            rjson_0.2.21               
[123] lifecycle_1.0.4             nlme_3.1-164               
[125] jsonlite_1.8.8              viridisLite_0.4.2          
[127] fansi_1.0.6                 pillar_1.9.0               
[129] lattice_0.22-5              KEGGREST_1.38.0            
[131] fastmap_1.1.1               httr_1.4.7                 
[133] survival_3.5-7              glue_1.7.0                 
[135] png_0.1-8                   bit_4.0.5                  
[137] stringi_1.8.3               blob_1.2.4                 
[139] RcppHNSW_0.5.0              memoise_2.0.1              
[141] dplyr_1.1.4                 irlba_2.3.5.1              
[143] future.apply_1.11.1 

[timoast]

It looks like these links fall outside of the plotted region, if you extend the region you should see the links appear

[RaghadShu]

Thanks! One way that works is to use the seqnames from the Links(seurat) as the region, and the name of the gene as the features. This worked out pretty nicely. Thanks!

[RaghadShu]

Hi,
I have some followup questions that I couldn't find online answers for. Also I'm new to scATAC so please bear with me:

1. From documentation, Links(seurat) outputs this table:

GRanges object with 2 ranges and 5 metadata columns:
      seqnames            ranges strand |     score        gene
         <Rle>         <IRanges>  <Rle> | <numeric> <character>
  [1]     chr2 20225433-20579145      * | 0.0554166        SDC1
  [2]    chr22 28800597-28823298      * | 0.0522171        XBP1
                        peak    zscore      pvalue
                 <character> <numeric>   <numeric>
  [1] chr2-20578650-20579640   3.73339 9.44593e-05
  [2] chr22-28822785-28823..   3.96908 3.60757e-05
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

Where "peak" is the linked peak and "gene" is the linked gene. What is the information given by the first two columns? seqnames-ranges? Is it the region that encompasses both, the gene and the linked peak?

2. Now that I have a linked peak, how can I annotate with a motif and get the transcription factor it binds to? I tried something like this:

peak_coordinate <- "chr2-20578650-20579640"
motifs <- FindMotifs(object = seurat, assay = "peakunion_MACS2" ,features = peak_coordinate)
top_motif <- head(motifs, 3)

and get something like this:

            motif observed background percent.observed percent.background
MA0819.1 MA0819.1        1        995              100             2.4875
MA1562.1 MA1562.1        1       1214              100             3.0350
MA0616.2 MA0616.2        1       1570              100             3.9250
         fold.enrichment   pvalue motif.name  p.adjust
MA0819.1        40.20101 0.024875      CLOCK 0.9175514
MA1562.1        32.94893 0.030350      SOX14 0.9175514
MA0616.2        25.47771 0.039250       HES2 0.9175514

But I am not sure if CLOCK is the transcription factor that regulates SDC1, or is CLOCK just the name of the motif but any transcription factor can bind to it?

Thanks in advance,
Raghad