# Analyzing R-loop data with RLSeq

Henry Miller1,2

1Alex Bishop Laboratory, UT Health San Antonio
2Bioinformatics Research Network

#### 26 October 2021

#### Abstract

This vignette covers basic usage of RLSeq for evaluating data quality and analyzing R-loop locations. RLSeq is part of RLSuite, an R-loop analysis toolchain. RLSuite also includes RLHub, RLBase, and RLPipes.

#### Package

RLSeq 1.0.0

# 1 Introduction

RLSeq is a package for analyzing R-loop mapping data sets, and it is a core component of the RLSuite toolchain. It serves two primary purposes: (1) to facilitate the evaluation of data quality, and (2) to enable R-loop data analysis in the context of genomic annotations and the public data sets in RLBase. The main analysis steps can be conveniently run using the `RLSeq()` function. Then, an HTML report can be generated using the `report()` function. Individual steps of this pipeline are also accessible through separate functions which provide custom analysis capabilities.

This vignette will showcase the primary functionality of RLSeq with data from a publicly-available R-loop data mapping study in Ewing sarcoma cell lines, GSE68845. We have selected two DNA-RNA Immunoprecipitation sequencing (DRIP-seq) samples for demonstration purposes: (1) SRX1025890, a positive R-loop mapping sample (“POS”; condition: S9.6 -RNaseH1), and (2) SRX1025892, a negative control (“NEG”; condition S9.6 +RNaseH1). We will begin by showing a quick-start analysis on SRX1025890, and then we will proceed to discuss, in detail, the specific steps of this analysis with both samples.

# 2 Quick-start

Here, we demonstrate a simple analysis workflow which utilizes a publicly-available data set stored in RLBase (a database of R-loop consensus regions and R-loop-mapping experiments, also part of RLSuite). The commands below download these data, run `RLSeq()`, and generate the HTML report.

    # Peaks and coverage can be found in RLBase
    rlbase <- "https://rlbase-data.s3.amazonaws.com"
    pks <- file.path(rlbase, "peaks", "SRX1025890_hg38.broadPeak")
    cvg <- file.path(rlbase, "coverage", "SRX1025890_hg38.bw")

    # Initialize data in the RLRanges object.
    # Metadata is optional, but improves the interpretability of results
    rlr <- RLRanges(
      peaks = pks,
      coverage = cvg,
      genome = "hg38",
      mode = "DRIP",
      label = "POS",
      sampleName = "TC32 DRIP-Seq"
    )

    # The RLSeq command performs all analyses
    rlr <- RLSeq(rlr)

    # Generate an html report
    report(rlr, reportPath = "rlseq_report_example.html")

The report generated by this code is found here.

# 3 Preliminary

## 3.1 Installation

RLSeq should be installed alongside RLHub to facilitate access to the data required for annotation and analysis. When downloading RLSeq from bioconductor, RLHub is already included.

    if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")

    BiocManager::install("RLSeq")

Both packages can also be installed from github.

    library(remotes)
    install_github("Bishop-Laboratory/RLHub")
    install_github("Bishop-Laboratory/RLSeq")

## 3.2 Obtaining data

RLSeq is compatible with R-loop data generated from a variety of pipelines and tools. However, it is strongly recommended that you use RLPipes, a snakemake-based CLI pipeline tool built specifically for upstream processing of R-loop datasets.

RLPipes can be installed using mamba or conda (slower).

    # conda install -c conda-forge mamba
    mamba create -n rlpipes -c bioconda -c conda-forge rlpipes
    conda activate rlpipes

A typical config file (`CSV`) should be written as such:

experiment
SRX1025890
SRX1025892

And then the pipeline can be run.

    RLPipes build -m DRIP rseq_out/ tests/test_data/samples.csv
    RLPipes run rseq_out/

The resulting directory will contain `peaks/`, `coverage/`, `bam/`, and other processed data sets which are used in downstream analysis.

Note: If you choose to use a different pipeline, use macs2/macs3 for peak calling to ensure compatibility with RLBase.

# 4 End-to-end RLSeq

Here, we describe each step of the analysis pipeline which is run as part of the `RLSeq()` command.

    library(RLSeq)
    library(dplyr)

## 4.1 Data sets

For this example, we will be using DRIP-Seq data from a 2018 Nature paper on R-loops in Ewing sarcoma (Gorthi et al. 2018). The sample has been IP’d for R-loops (S9.6 -RNaseH1; label: “POS”). The data was processed using RLPipes and uploaded to RLBase. Peaks are converted to `GRanges` objects using a helper function from `regioneR`. URLs and file paths for peak files can also be supplied directly to without this step.

    rlbase <- "https://rlbase-data.s3.amazonaws.com"

    # Get peaks and coverage
    s96Pks <- regioneR::toGRanges(file.path(rlbase, "peaks", "SRX1025890_hg38.broadPeak"))
    s96Cvg <- file.path(rlbase, "coverage", "SRX1025890_hg38.bw")

For demonstration purposes, only 10000 ranges are analyzed here.

    # For expediency, peaks we filter and down-sampled to the top 10000 by padj (V9)
    # This is not necessary as part of the typical workflow, however
    s96Pks <- s96Pks[s96Pks$V9 > 2,]
    s96Pks <- s96Pks[sample(names(s96Pks), 10000)]

Finally, `RLRanges` objects were constructed. These are the primary objects used in all `RLSeq` functions. `RLRanges` are an extension of `GRanges` which provide additional metadata and validation functions.

    ## Build RLRanges ##
    # S9.6 -RNaseH1
    rlr <- RLRanges(
      peaks = s96Pks,
      coverage = s96Cvg,
      genome = "hg38",
      mode = "DRIP",
      label = "POS",
      sampleName = "TC32 DRIP-Seq"
    )

## 4.2 Sample quality

Sample quality is assessed by analyzing the association of peaks with R-loop-forming sequences (RLFS). RLFS are genomic sequences that favor the formation of R-loops (Jenjaroenpun et al. 2015). While R-loops can form outside RLFS, there is a strong relationship between them, which provides an unbiased test of whether a set of peaks actually represents successful R-loop mapping.

### 4.2.1 Permutation tests

RLSeq first implements a permutation test to evaluate the enrichment of peaks within RLFS and build a Z-score distribution around RLFS sites.

    # Analyze RLFS for positive sample
    rlr <- analyzeRLFS(rlr, quiet = TRUE)

The resulting objects now contain the permutation test results. These results can be easily visualized with the `plotRLFSRes` function.

    plotRLFSRes(rlr)

Figure 1: Plot of permutation test results (S9.6 -RNaseH1)

### 4.2.2 Quality classification

The quality classifier is an ensemble model based on an online-learning scheme. It predicts “POS” for samples which are predicted to show robust R-loop mapping and “NEG” for samples which are not. The latest version can be accessed via RLHub. For greater detail, please see the `RLHub::modes` reference. To apply the model and predict sample quality, use the `predictCondition()` function.

    # Predict
    rlr <- predictCondition(rlr)

The results from testing our example samples:

    # Access results
    s96_pred <- rlresult(rlr, "predictRes")
    cat("Prediction: ", s96_pred$prediction)

    ## Prediction:  POS

## 4.3 Feature enrichment

The feature enrichment test assesses the enrichment of genomic features within a supplied R-loop dataset. The function queries the RLHub annotation database to retrieve genomic features, and then it performs fisher’s exact test and the relative distance test to assess feature enrichment (Favorov et al. 2012). For a description of all features, see RLHub::annotations.

    # Perform test
    rlr <- featureEnrich(
      object = rlr,
      quiet = TRUE
    )

The results:

    # View Top Results
    annoResS96 <- rlresult(rlr, "featureEnrichment")
    annoResS96 %>%
      relocate(contains("fisher"), .after = type) %>%
      arrange(desc(stat_fisher_rl))

    ## # A tibble: 50 x 13
    ##    db      type  stat_fisher_rl stat_fisher_shuf pval_fisher_rl pval_fisher_shuf
    ##    <chr>   <chr>          <dbl>            <dbl>          <dbl>            <dbl>
    ##  1 Repeat~ SINE          Inf               1.11       2.55e-128           0.128
    ##  2 Transc~ Intr~         Inf               0.994      2.23e-308           0.760
    ##  3 snoRNA~ CDBox          21.6             1.79       2.01e- 23           0.309
    ##  4 Transc~ Exon           16.4             1.00       2.23e-308           0.967
    ##  5 Transc~ TSS            14.8             0.995      2.23e-308           0.921
    ##  6 Transc~ TTS            14.4             0.936      2.23e-308           0.108
    ##  7 knownG~ prot~          12.1             0.970      2.23e-308           0.209
    ##  8 Encode~ enhP           11.5             0.932      2.23e-308           0.0937
    ##  9 PolyA   poly~          10.3             0.946      2.23e-308           0.503
    ## 10 Transc~ thre~           8.95            0.960      2.23e-308           0.536
    ## # ... with 40 more rows, and 7 more variables: num_tested_peaks <int>,
    ## #   num_total_peaks <int>, num_tested_anno_ranges <int>,
    ## #   num_total_anno_ranges <int>, avg_reldist_rl <dbl>, avg_reldist_shuf <dbl>,
    ## #   pval_reldist <dbl>

From the results, we see that there is high enrichment within genic features, such as exons and introns.

### 4.3.1 Visualization of enrichment results

RLSeq provides a helper function, `plotEnrichment`, to facilitate the visualization of enrichment results.

    pltlst <- plotEnrichment(rlr)

This returns a list of plots named according to the corresponding annotation database. For example, Encode cis-regulatory elements (CREs):

    pltlst$Encode_CREs

Note: Caveat on data range

A limitation of this approach is that Fisher’s exact test sometimes returns `Inf` or `-Inf` for the statistic (odds ratio). While these results are useful in demonstrating robust enrichment or non-enrichment, they are difficult to plot in a meaningful way. As a compromise, `plotEnrichment` sets a limited data range of -10 through 15. These values were chosen because they encompass every finite value that can be returned from the implementation of Fisher’s test which RLSeq uses. In the above plots `Inf` results are shown on the y-axis at value `15` and, likewise, `-Inf` is shown at `-10`.

## 4.4 Correlation analysis

Correlation analysis finds inter-sample correlation coefficients of bin-level R-loop signal around gold-standard R-loop sites (sites profiled using ultra-long-read R-loop mapping – “SMRF-Seq”) (Chédin et al. 2021). This analysis helps to answer the question “how well does my data agree with previous results in public datasets?” For greater detail, please refer to the documentation for `corrAnalyze()`.

Note: `corrAnalyze()` does not work on Windows OS systems.

    # corrAnalyze does not work on Windows OS
    if (.Platform$OS.type != "windows") {
      rlr <- corrAnalyze(rlr)
    }

The results of this analysis are visualized using `corrHeatmap`.

    # corrAnalyze does not work on Windows OS
    if (.Platform$OS.type != "windows") {
      corrHeatmap(rlr)
    }

These results demonstrate that our sample correlates well with similar DRIP-Seq data sets.

## 4.5 Gene Annotation

Gene annotations are automatically downloaded using `AnnotationHub()` and then intersected with RLRanges.

    rlr <- geneAnnotation(rlr)

## 4.6 R-Loop Region Test

R-loop regions are consensus R-loop-forming sites discovered from analyzing all high-confidence R-loop mapping samples in RLBase. The `rlRegionTest()` analyzes the enrichment of the ranges in our `RLRanges` object with these consensus R-loop sites, which, like correlation analysis, also helps answer the question “how well does my data agree with previous results?” For greater detail, please refer to the documentation for `rlRegionTest()`.

    rlr <- rlRegionTest(rlr)

The test results can be easily visualized in the following manner.

    # Plot overlap
    plotRLRegionOverlap(
      object = rlr,

      # Arguments for VennDiagram::venn.diagram()
      fill = c("#9ad9ab", "#9aa0d9"),
      main.cex = 2,
      cat.pos = c(-40, 40),
      cat.dist=.05,
      margin = .05
    )

# 5 Accessing RLBase data

For convenience, we also provide pre-analyzed `RLRanges` objects for every sample in RLBase. To access them, you need only provide the ID of the sample which you want to obtain data from. These IDs, along with other metadata, are listed in `RLHub::rlbase_samples()`.

    rlr <- RLRangesFromRLBase(acc = "SRX1025890")
    rlr

    ## GRanges object with 107029 ranges and 6 metadata columns:
    ##          seqnames            ranges strand |                     V4        V5
    ##             <Rle>         <IRanges>  <Rle> |            <character> <integer>
    ##        1     chr1       10034-10345      * | /home/UTHSCSA/miller..        40
    ##        2     chr1     180610-181657      * | /home/UTHSCSA/miller..        53
    ##        3     chr1     182752-182950      * | /home/UTHSCSA/miller..        38
    ##        4     chr1     184149-184628      * | /home/UTHSCSA/miller..        28
    ##        5     chr1     629787-630103      * | /home/UTHSCSA/miller..        42
    ##      ...      ...               ...    ... .                    ...       ...
    ##   107025     chrY 11293162-11294964      * | /home/UTHSCSA/miller..        38
    ##   107026     chrY 11295281-11296131      * | /home/UTHSCSA/miller..        36
    ##   107027     chrY 11297625-11297944      * | /home/UTHSCSA/miller..        23
    ##   107028     chrY 11301521-11301748      * | /home/UTHSCSA/miller..        20
    ##   107029     chrY 26641735-26642007      * | /home/UTHSCSA/miller..        15
    ##                   V6        V7        V8      qval
    ##          <character> <numeric> <numeric> <numeric>
    ##        1           .   4.17431   6.16893   4.04326
    ##        2           .   4.63235   7.59510   5.31977
    ##        3           .   4.21966   5.95157   3.81400
    ##        4           .   3.57860   4.84694   2.87393
    ##        5           .   2.59231   6.38817   4.27687
    ##      ...         ...       ...       ...       ...
    ##   107025           .   2.33249   5.86578   3.81115
    ##   107026           .   2.18689   5.72732   3.67967
    ##   107027           .   2.45115   4.20891   2.32963
    ##   107028           .   2.98060   3.91975   2.09600
    ##   107029           .   2.80910   3.30759   1.56104
    ##   -------
    ##   seqinfo: 595 sequences (1 circular) from hg38 genome
    ##
    ## SRX1025890:
    ##   Mode: DRIP
    ##   Genome: hg38
    ##   Label: POS
    ##
    ## RLSeq Results Available:
    ##   featureEnrichment, correlationMat, rlfsRes, geneAnnoRes, predictRes, rlRegionRes
    ##
    ## prediction: POS

# 6 Session

Session info

    sessionInfo()

    ## R version 4.1.1 (2021-08-10)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server x64 (build 17763)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.1252
    ## [3] LC_MONETARY=English_United States.1252
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.1252
    ##
    ## attached base packages:
    ## [1] stats4    stats     graphics  grDevices utils     datasets  methods
    ## [8] base
    ##
    ## other attached packages:
    ##  [1] GenomicFeatures_1.46.0 AnnotationDbi_1.56.0   Biobase_2.54.0
    ##  [4] GenomicRanges_1.46.0   RLHub_0.99.6           GenomeInfoDb_1.30.0
    ##  [7] IRanges_2.28.0         S4Vectors_0.32.0       BiocGenerics_0.40.0
    ## [10] dplyr_1.0.7            RLSeq_1.0.0            BiocStyle_2.22.0
    ##
    ## loaded via a namespace (and not attached):
    ##   [1] backports_1.2.1               AnnotationHub_3.2.0
    ##   [3] BiocFileCache_2.2.0           plyr_1.8.6
    ##   [5] splines_4.1.1                 BiocParallel_1.28.0
    ##   [7] listenv_0.8.0                 ggplot2_3.3.5
    ##   [9] valr_0.6.3                    digest_0.6.28
    ##  [11] yulab.utils_0.0.4             foreach_1.5.1
    ##  [13] htmltools_0.5.2               magick_2.7.3
    ##  [15] fansi_0.5.0                   magrittr_2.0.1
    ##  [17] memoise_2.0.0                 BSgenome_1.62.0
    ##  [19] tzdb_0.1.2                    aws.signature_0.6.0
    ##  [21] readr_2.0.2                   recipes_0.1.17
    ##  [23] globals_0.14.0                Biostrings_2.62.0
    ##  [25] gower_0.2.2                   matrixStats_0.61.0
    ##  [27] prettyunits_1.1.1             colorspace_2.0-2
    ##  [29] blob_1.2.2                    rappdirs_0.3.3
    ##  [31] xfun_0.27                     callr_3.7.0
    ##  [33] crayon_1.4.1                  RCurl_1.98-1.5
    ##  [35] jsonlite_1.7.2                survival_3.2-13
    ##  [37] iterators_1.0.13              glue_1.4.2
    ##  [39] gtable_0.3.0                  ipred_0.9-12
    ##  [41] zlibbioc_1.40.0               XVector_0.34.0
    ##  [43] DelayedArray_0.20.0           kernlab_0.9-29
    ##  [45] future.apply_1.8.1            scales_1.1.1
    ##  [47] futile.options_1.0.1          DBI_1.1.1
    ##  [49] Rcpp_1.0.7                    xtable_1.8-4
    ##  [51] progress_1.2.2                gridGraphics_0.5-1
    ##  [53] bit_4.0.4                     DT_0.19
    ##  [55] lava_1.6.10                   prodlim_2019.11.13
    ##  [57] htmlwidgets_1.5.4             httr_1.4.2
    ##  [59] ellipsis_0.3.2                farver_2.1.0
    ##  [61] pkgconfig_2.0.3               XML_3.99-0.8
    ##  [63] nnet_7.3-16                   sass_0.4.0
    ##  [65] dbplyr_2.1.1                  utf8_1.2.2
    ##  [67] caret_6.0-90                  ggplotify_0.1.0
    ##  [69] labeling_0.4.2                tidyselect_1.1.1
    ##  [71] rlang_0.4.12                  reshape2_1.4.4
    ##  [73] later_1.3.0                   munsell_0.5.0
    ##  [75] BiocVersion_3.14.0            tools_4.1.1
    ##  [77] cachem_1.0.6                  cli_3.0.1
    ##  [79] ggprism_1.0.3                 generics_0.1.1
    ##  [81] RSQLite_2.2.8                 ExperimentHub_2.2.0
    ##  [83] broom_0.7.9                   aws.s3_0.3.21
    ##  [85] evaluate_0.14                 stringr_1.4.0
    ##  [87] fastmap_1.1.0                 yaml_2.2.1
    ##  [89] processx_3.5.2                ModelMetrics_1.2.2.2
    ##  [91] knitr_1.36                    bit64_4.0.5
    ##  [93] randomForest_4.6-14           purrr_0.3.4
    ##  [95] KEGGREST_1.34.0               pbapply_1.5-0
    ##  [97] future_1.22.1                 nlme_3.1-153
    ##  [99] mime_0.12                     formatR_1.11
    ## [101] xml2_1.3.2                    caretEnsemble_2.0.1
    ## [103] biomaRt_2.50.0                compiler_4.1.1
    ## [105] filelock_1.0.2                curl_4.3.2
    ## [107] png_0.1-7                     interactiveDisplayBase_1.32.0
    ## [109] tibble_3.1.5                  bslib_0.3.1
    ## [111] stringi_1.7.5                 futile.logger_1.4.3
    ## [113] highr_0.9                     ps_1.6.0
    ## [115] lattice_0.20-45               Matrix_1.3-4
    ## [117] vctrs_0.3.8                   pillar_1.6.4
    ## [119] lifecycle_1.0.1               BiocManager_1.30.16
    ## [121] jquerylib_0.1.4               data.table_1.14.2
    ## [123] bitops_1.0-7                  httpuv_1.6.3
    ## [125] rtracklayer_1.54.0            R6_2.5.1
    ## [127] BiocIO_1.4.0                  bookdown_0.24
    ## [129] promises_1.2.0.1              gridExtra_2.3
    ## [131] parallelly_1.28.1             codetools_0.2-18
    ## [133] lambda.r_1.2.4                MASS_7.3-54
    ## [135] assertthat_0.2.1              SummarizedExperiment_1.24.0
    ## [137] rjson_0.2.20                  withr_2.4.2
    ## [139] regioneR_1.26.0               GenomicAlignments_1.30.0
    ## [141] Rsamtools_2.10.0              GenomeInfoDbData_1.2.7
    ## [143] parallel_4.1.1                hms_1.1.1
    ## [145] VennDiagram_1.6.20            grid_4.1.1
    ## [147] rpart_4.1-15                  timeDate_3043.102
    ## [149] tidyr_1.1.4                   class_7.3-19
    ## [151] rmarkdown_2.11                MatrixGenerics_1.6.0
    ## [153] pROC_1.18.0                   base64enc_0.1-3
    ## [155] shiny_1.7.1                   lubridate_1.8.0
    ## [157] restfulr_0.0.13

# References

Chédin, Frédéric, Stella R Hartono, Lionel A Sanz, and Vincent Vanoosthuyse. 2021. “Best Practices for the Visualization, Mapping, and Manipulation of R-Loops.” The EMBO Journal 40 (4). https://doi.org/10.15252/embj.2020106394.

Favorov, Alexander, Loris Mularoni, Leslie M. Cope, Yulia Medvedeva, Andrey A. Mironov, Vsevolod J. Makeev, and Sarah J. Wheelan. 2012. “Exploring Massive, Genome Scale Datasets with the Genometricorr Package.” Edited by Hilmar Lapp. PLoS Computational Biology 8 (5): e1002529. https://doi.org/10.1371/journal.pcbi.1002529.

Gorthi, Aparna, July Carolina Romero, Eva Loranc, Lin Cao, Liesl A. Lawrence, Elicia Goodale, Amanda Balboni Iniguez, et al. 2018. “EWSFLI1 Increases Transcription to Cause R-Loops and Block Brca1 Repair in Ewing Sarcoma.” Nature 555 (7696): 387–91. https://doi.org/10.1038/nature25748.

Jenjaroenpun, Piroon, Thidathip Wongsurawat, Surya Pavan Yenamandra, and Vladimir A. Kuznetsov. 2015. “QmRLFS-Finder: A Model, Web Server and Stand-Alone Tool for Prediction and Analysis of R-Loop Forming Sequences.” Nucleic Acids Research 43 (W1): W527–W534. https://doi.org/10.1093/nar/gkv344.