knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file(), 
                     cache = FALSE)

Summary

The purpose of this analysis was to find the locations of R-loops in EUFA cells with and without BRCA2 and PAF1 complementation using the DRIP-Seq data generated from this study.

Questions to answer

From GitHub:

Questions:

  1. How does BRCA2 and PAF1 status impact R-loops? Are there regions of the genome which are particularly impacted?
  2. What genes are R-loops differentially found in? What pathways do they relate to?
  3. How do these results relate to BRCA2, PAF1, and XRN2 binding sites? How do they relate to the siXRN2 DRIP-Seq results?

Overall Results

From the analysis, we found:

BRCA2 and PAF1 complementation resolves R-loops throughout the EUFA genome. Interestingly, almost all the R-loops lost by BRCA2 complementation are also lost by PAF1 complementation, indicating an immense amount of overlap. The differential loss of R-loop abundance occurs at genes which are related to a variety of biology processes, such as DNA Repair and Ribosome. These genes are also bound by transcription factors such as E2F6, ZBTB7A, TAF1, and BRCA1.

From comparison with ChIP-Seq data, it was found that this happens are sites which are likely bound by both PAF1 and BRCA2. However, limitations in the approach taken here make it challenging to directly observe every instance where this happens. After overlapping these results with the known binding sites of XRN2, BRCA1, and PAF1, we found that a majority of differentially abundant R-loops also overlap with at least one of these factors.

Where multiple overlaps occur, there was strong enrichment for pathways related to MYC, E2F6, Ribosome biogenesis, breast cancer, estrogen treatment, UBTF, ATF, BRCA1 and others. These results suggest a strong relationship between XRN2, PAF1, and BRCA2 in the control of R-loop formation in promoter proximal regions of important gene sets for biological pathways relevant to the biology of interest.

Analysis

# Set the library path to the conda env that has DiffBind 2.16
library(tidyverse)
library(ChIPpeakAnno)
library(ChIPseeker)
library(clusterProfiler)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(EnsDb.Hsapiens.v86)
library(org.Hs.eg.db)
library(biomaRt)
library(annotatr)
library(enrichR)
setEnrichrSite("Enrichr") # Human genes
dbs <- listEnrichrDbs()
if (! file.exists("misc/annotatr_annotations.rda")) {
  annots = c('hg38_basicgenes', "hg38_genes_promoters", 'hg38_genes_intergenic', 
           "hg38_genes_intronexonboundaries", "hg38_genes_exonintronboundaries",
           'hg38_enhancers_fantom', "hg38_genes_firstexons" )
  annotations = build_annotations(genome = "hg38", annotations = annots)
  annotations$type[annotations$type == "hg38_genes_1to5kb"] <- "hg38_genes_upstream-1to5kb"
  annotations <- annotations[which(width(annotations) >= 1),]
  annotation_order <- c(
    "hg38_enhancers_fantom",
    "hg38_genes_upstream-1to5kb",
    "hg38_genes_promoters",
    "hg38_genes_5UTRs",
    "hg38_genes_firstexons",
    "hg38_genes_exons",
    "hg38_genes_exonintronboundaries",
    "hg38_genes_introns",
    "hg38_genes_intronexonboundaries",
    "hg38_genes_3UTRs",
    "hg38_genes_intergenic"
  )
  save(annotations, annotation_order, file = "misc/annotatr_annotations.rda")
} else {
  load("misc/annotatr_annotations.rda")
}

txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
annoData <- ChIPpeakAnno::toGRanges(EnsDb.Hsapiens.v86, feature="gene")

EUFA R-loop location analysis

Preliminary

Wrangling the peaks

files <- list(
  "EUFA_rep1" = "data/drip_seq/eufa_eufab2_drip/drip_seq/rseq_output/EUFA_1_S35_L004/peaks_macs_unstranded/EUFA_1_S35_L004_hg38_peaks.xls",
  "EUFA_rep2" = "data/drip_seq/eufa_eufab2_drip/drip_seq/rseq_output/EUFA_2_S36_L004/peaks_macs_unstranded/EUFA_2_S36_L004_hg38_peaks.xls"
)
EUFAgr <- lapply(files, function(file) {
  tmpdf <- read_tsv(file, skip = 30) %>%
    rename_with(~ gsub(pattern = "-log10\\((.+)\\)",
                       replacement = "\\1", .x)) %>%
    mutate(name = gsub(pattern = "called_peaks/", replacement = "", name)) %>%
    as.data.frame()
  gr <- toGRanges(tmpdf)
  gr <- gr[which(width(gr) < 10000),]
})
sapply(names(EUFAgr), function(peakName){
  peaks <- EUFAgr[[peakName]]
  paste0("number of peaks: ", length(peaks$length))
})
##                EUFA_rep1                EUFA_rep2 
## "number of peaks: 70319" "number of peaks: 62655"

Overlap of replicates

# Find the overlap of the replicates
olEUFA <- findOverlapsOfPeaks(EUFAgr)
olEUFA <- addMetadata(olEUFA, colNames="qvalue", FUN=mean) 

Venn diagram of replicate overlap

tmp = makeVennDiagram(olEUFA, fill = c("skyblue", "firebrick"),
                margin = .05)  # Decent sized overlap

Comparison of P Adjusted Value between overlapping/non-overlapping

macs2 assigns an adjusted p value to each peak, indicating the level of confidence in the peak calling. By calculating the overlap of peaks between biological replicates, we are (in theory) finding the genuine peaks. We can verify this by comparing the P adjusted values of the peaks that overlapped to the ones which didn’t (below).

olQV <- olEUFA$peaklist$`EUFA_rep1///EUFA_rep2`$qvalue
nonolQV <- c(olEUFA$peaklist$EUFA_rep2$qvalue,
             olEUFA$peaklist$EUFA_rep1$qvalue)
tibble(
  qval = c(olQV, nonolQV),
  group = c(rep("Overlapping", length(olQV)), 
            rep("Non-overlapping", length(nonolQV)))
) %>%
  ggplot(mapping = aes(y = qval, x = group, fill = group)) +
  geom_boxplot() +
  ylab("P Adjusted Value (-log10)") +
  xlab(NULL) +
  theme_bw(base_size = 15) +
  ggpubr::rremove("legend") +
  ggpubr::stat_compare_means(comparisons = list(c("Overlapping", "Non-overlapping")),
                             label = "p.signif", size = 5) +
  scale_y_continuous(limits = c(-5, 220)) +
  labs(title = "Overlapping/Non-overlapping R-loop sites P-Adj Value")

The results indicate that the overlapping peaks are only a slightly more robust measure of EUFA R-loop locations, indicated by the fact that they are more significant. However, there are other ways to assess this…

Comparing the R-loop profile of overlapping and non-overlapping peaks

Another way to assess the robustness of our overlapping peakset is to check the binding profile using a tool like ChIPseeker. In particular, we expect to see that most genuine EUFA R-loop locations are found closer to the TSS and enriched in genic regions.

Metaplot around TSS

plotAvgProf(tagMatrixListEUFA, xlim=c(-3000, 3000), facet = "row") +
  labs(title = "Overlapping vs Non-overlapping EUFA R-loop peaks around TSS")

Distance to TSS plot

plotDistToTSS(peakAnnoListEUFA) + 
  labs(title = "EUFA peak locations relative to TSS")

Feature distribution plot

plotAnnoBar(peakAnnoListEUFA) +
  labs(title = "Feature Overlap with EUFA Peaks")

However, we find that this isn’t completely clearcut. It seems that there are genuine R-loops occuring in the non-overlapping group. Therefore, the final peakset will be considered as R-loop peaks which either:

  1. Are found in both replicates (overlapping peaks)
  2. Have a P adjusted value < the 50% quantile of all p values
cutoff <- quantile(c(olEUFA$all.peaks$EUFA_rep1$qvalue, olEUFA$all.peaks$EUFA_rep2$qvalue))[3]
nonOlKeep <- olEUFA$uniquePeak[olEUFA$uniquePeaks$qvalue > cutoff, c("qvalue")]
EUFApeaks <- c(nonOlKeep, olEUFA$peaklist$`EUFA_rep1///EUFA_rep2`[,c("qvalue")])
names(EUFApeaks) <- NULL

Save EUFA R-loop peaks

EUFApeaks <- keepStandardChromosomes(EUFApeaks, pruning.mode = "coarse")
rtracklayer::export(EUFApeaks, con = "analysis/diff_drip_brca2_paf1/results/EUFApeaks.bed")

Results

What kinds of genomic features do these R-loops appear in?

EUFA_annot = annotate_regions(
  regions = EUFApeaks,
  annotations = annotations,
  ignore.strand = TRUE,
  quiet = F)
regioneR::getMask(genome = "hg38")
## GRanges object with 0 ranges and 0 metadata columns:
##    seqnames    ranges strand
##       <Rle> <IRanges>  <Rle>
##   -------
##   seqinfo: no sequences
EUFArandom <- regioneR::randomizeRegions(
  A = EUFApeaks,  
  genome = "hg38"
)
EUFArandom_annot = annotate_regions(
  regions = EUFArandom,
  annotations = annotations,
  ignore.strand = TRUE,
  quiet = FALSE)
plot_annotation(
  annotated_regions = EUFA_annot,  quiet = FALSE,
  annotated_random = EUFArandom_annot,
  annotation_order = annotation_order,
  plot_title = 'EUFA R-loop Site Feature Overlaps (vs Random)',
  x_label = 'Annotations',
  y_label = 'Count') 

What genes are they in? What pathways are these genes involved in?

EUFAgeneAnno <- annotatePeakInBatch(EUFApeaks, output = "overlap", maxgap = 1000,
                                  AnnotationData=annoData)
EUFAgeneAnno <- addGeneIDs(EUFAgeneAnno,
                         "org.Hs.eg.db",
                         IDs2Add = "SYMBOL")

EUFA R-loop genes saved to TSV file.

EUFAgenes <- as.data.frame(EUFAgeneAnno) %>%
  dplyr::filter(! is.na(SYMBOL)) %>%
  dplyr::select(SYMBOL) %>%
  dplyr::filter(! grepl(SYMBOL, pattern = ";|/")) %>%
  dplyr::distinct(SYMBOL) %>%
  write_tsv(file = "analysis/diff_drip_brca2_paf1/results/EUFA_bound_genes.tsv") %>%
  pull(SYMBOL) 

Pathway enrichment with Enrichr

The ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X and KEGG_2019_HUMAN databases were queried to enrich for relevant gene sets. The analysis with all possible genes sets is available permanently at this link

  1. KEGG Pathways
KEGG_EUFA <- EUFAenriched$KEGG_2019_Human
KEGG_EUFA %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top KEGG Pathways Enriched in EUFA R-loop genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

  1. ChEA Pathways (Transcription Factors)
ChEA_EUFA <- EUFAenriched$`ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X`
ChEA_EUFA %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top ChEA Pathways Enriched in EUFA R-loop genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

NOTE: There may be too many genes (> 11,000) in the list for meaningful enrichment results… One way to test the purity of the analysis is to also enrich a randomized set of R-loops. Random enricher link

EUFArandom <- keepStandardChromosomes(EUFArandom, pruning.mode = "coarse")
EUFArandgeneAnno <- annotatePeakInBatch(EUFArandom, output = "overlap", maxgap = 1000,
                                  AnnotationData=annoData)
EUFArandgeneAnno <- addGeneIDs(EUFArandgeneAnno,
                         "org.Hs.eg.db",
                         IDs2Add = "SYMBOL")
EUFArandgenes <- as.data.frame(EUFArandgeneAnno) %>%
  dplyr::filter(! is.na(SYMBOL)) %>%
  dplyr::select(SYMBOL) %>%
  dplyr::filter(! grepl(SYMBOL, pattern = ";|/")) %>%
  dplyr::distinct(SYMBOL) %>%
  write_tsv(file = "analysis/diff_drip_brca2_paf1/results/EUFA_bound_randomized_genes.tsv") %>%
  pull(SYMBOL) 
  1. Randomized R-loops KEGG Pathways
KEGG_EUFA <- EUFArandenriched$KEGG_2019_Human
KEGG_EUFA %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top KEGG Pathways Enriched in EUFA R-loop randomized genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

  1. Randomized R-loops ChEA Pathways (Transcription Factors)
ChEA_EUFA <- EUFArandenriched$`ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X`
ChEA_EUFA %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top ChEA Pathways Enriched in EUFA randomized R-loop genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

EUFAB2 R-loop location analysis

Preliminary

Wrangling the peaks

files <- list(
  "EUFAB2_rep1" = "data/drip_seq/eufa_eufab2_drip/drip_seq/rseq_output/EUFA_BRCA2_1_S37_L004/peaks_macs_unstranded/EUFA_BRCA2_1_S37_L004_hg38_peaks.xls",
  "EUFAB2_rep2" = "data/drip_seq/eufa_eufab2_drip/drip_seq/rseq_output/EUFA_BRCA2_2_S38_L004/peaks_macs_unstranded/EUFA_BRCA2_2_S38_L004_hg38_peaks.xls"
)
EUFAB2gr <- lapply(files, function(file) {
  tmpdf <- read_tsv(file, skip = 30) %>%
    rename_with(~ gsub(pattern = "-log10\\((.+)\\)",
                       replacement = "\\1", .x)) %>%
    mutate(name = gsub(pattern = "called_peaks/", replacement = "", name)) %>%
    as.data.frame()
  gr <- toGRanges(tmpdf)
  gr <- gr[which(width(gr) < 10000),]
})
sapply(names(EUFAB2gr), function(peakName){
  peaks <- EUFAB2gr[[peakName]]
  paste0("number of peaks: ", length(peaks$length))
})
##              EUFAB2_rep1              EUFAB2_rep2 
## "number of peaks: 55312" "number of peaks: 62471"

Overlap of replicates

# Find the overlap of the replicates
olEUFAB2 <- findOverlapsOfPeaks(EUFAB2gr)
olEUFAB2 <- addMetadata(olEUFAB2, colNames="qvalue", FUN=mean) 

Venn diagram of replicate overlap

tmp = makeVennDiagram(olEUFAB2, fill = c("skyblue", "firebrick"),
                margin = .05)  # Decent sized overlap

Comparison of P Adjusted Value between overlapping/non-overlapping

macs2 assigns an adjusted p value to each peak, indicating the level of confidence in the peak calling. By calculating the overlap of peaks between biological replicates, we are (in theory) finding the genuine peaks. We can verify this by comparing the P adjusted values of the peaks that overlapped to the ones which didn’t (below).

olQV <- olEUFAB2$peaklist$`EUFAB2_rep1///EUFAB2_rep2`$qvalue
nonolQV <- c(olEUFAB2$peaklist$EUFAB2_rep2$qvalue,
             olEUFAB2$peaklist$EUFAB2_rep1$qvalue)
tibble(
  qval = c(olQV, nonolQV),
  group = c(rep("Overlapping", length(olQV)), 
            rep("Non-overlapping", length(nonolQV)))
) %>%
  ggplot(mapping = aes(y = qval, x = group, fill = group)) +
  geom_boxplot() +
  ylab("P Adjusted Value (-log10)") +
  xlab(NULL) +
  theme_bw(base_size = 15) +
  ggpubr::rremove("legend") +
  ggpubr::stat_compare_means(comparisons = list(c("Overlapping", "Non-overlapping")),
                             label = "p.signif", size = 5) +
  scale_y_continuous(limits = c(-5, 135)) +
  labs(title = "Overlapping/Non-overlapping R-loop sites P-Adj Value")

The results indicate that the overlapping peaks are only a bit more robust measure of EUFAB2 R-loop locations, indicated by the fact that they are more significant. However, there are other ways to assess this…

Comparing the R-loop profile of overlapping and non-overlapping peaks

Another way to assess the robustness of our overlapping peakset is to check the binding profile using a tool like ChIPseeker. In particular, we expect to see that most genuine EUFAB2 R-loop locations are found closer to the TSS and enriched in genic regions.

Metaplot around TSS

plotAvgProf(tagMatrixListEUFAB2, xlim=c(-3000, 3000), facet = "row") +
  labs(title = "Overlapping vs Non-overlapping EUFAB2 R-loop peaks around TSS")

Distance to TSS plot

plotDistToTSS(peakAnnoListEUFAB2) + 
  labs(title = "EUFAB2 peak locations relative to TSS")

Feature distribution plot

plotAnnoBar(peakAnnoListEUFAB2) +
  labs(title = "Feature Overlap with EUFAB2 Peaks")

However, we find that this isn’t completely clearcut. It seems that there are genuine R-loops occuring in the non-overlapping group. Therefore, the final peakset will be considered as R-loop peaks which either:

  1. Are found in both replicates (overlapping peaks)
  2. Have a P adjusted value < the 50% quantile of all p values
cutoff <- quantile(c(olEUFAB2$all.peaks$EUFAB2_rep1$qvalue, olEUFAB2$all.peaks$EUFAB2_rep2$qvalue))[3]
nonOlKeep <- olEUFAB2$uniquePeak[olEUFAB2$uniquePeaks$qvalue > cutoff, c("qvalue")]
EUFAB2peaks <- c(nonOlKeep, olEUFAB2$peaklist$`EUFAB2_rep1///EUFAB2_rep2`[,c("qvalue")])
names(EUFAB2peaks) <- NULL

Save EUFAB2 R-loop peaks

EUFAB2peaks <- keepStandardChromosomes(EUFAB2peaks, pruning.mode = "coarse")
rtracklayer::export(EUFAB2peaks, con = "analysis/diff_drip_brca2_paf1/results/EUFAB2peaks.bed")

Results

What kinds of genomic features do EUFAB2 R-loops appear in?

EUFAB2_annot = annotate_regions(
  regions = EUFAB2peaks,
  annotations = annotations,
  ignore.strand = TRUE,
  quiet = F)
regioneR::getMask(genome = "hg38")
## GRanges object with 0 ranges and 0 metadata columns:
##    seqnames    ranges strand
##       <Rle> <IRanges>  <Rle>
##   -------
##   seqinfo: no sequences
EUFAB2random <- regioneR::randomizeRegions(
  A = EUFAB2peaks,  
  genome = "hg38"
)
EUFAB2random_annot = annotate_regions(
  regions = EUFAB2random,
  annotations = annotations,
  ignore.strand = TRUE,
  quiet = FALSE)
plot_annotation(
  annotated_regions = EUFAB2_annot,  quiet = FALSE,
  annotated_random = EUFAB2random_annot,
  annotation_order = annotation_order,
  plot_title = 'EUFAB2 R-loop Site Feature Overlaps (vs Random)',
  x_label = 'Annotations',
  y_label = 'Count') 

What genes are they in? What pathways are these genes involved in?

EUFAB2geneAnno <- annotatePeakInBatch(EUFAB2peaks, output = "overlap", maxgap = 1000,
                                  AnnotationData=annoData)
EUFAB2geneAnno <- addGeneIDs(EUFAB2geneAnno,
                         "org.Hs.eg.db",
                         IDs2Add = "SYMBOL")

EUFAB2 R-loop genes saved to TSV file.

EUFAB2genes <- as.data.frame(EUFAB2geneAnno) %>%
  dplyr::filter(! is.na(SYMBOL)) %>%
  dplyr::select(SYMBOL) %>%
  dplyr::filter(! grepl(SYMBOL, pattern = ";|/")) %>%
  dplyr::distinct(SYMBOL) %>%
  write_tsv(file = "analysis/diff_drip_brca2_paf1/results/EUFAB2_bound_genes.tsv") %>%
  pull(SYMBOL) 

Pathway enrichment with Enrichr

The ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X and KEGG_2019_HUMAN databases were queried to enrich for relevant gene sets. The analysis with all possible genes sets is available permanently at this link

  1. KEGG Pathways
KEGG_EUFAB2 <- EUFAB2enriched$KEGG_2019_Human
KEGG_EUFAB2 %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top KEGG Pathways Enriched in EUFAB2 R-loop genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

  1. ChEA Pathways (Transcription Factors)
ChEA_EUFAB2 <- EUFAB2enriched$`ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X`
ChEA_EUFAB2 %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top ChEA Pathways Enriched in EUFAB2 R-loop genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

EUFAPAF1 R-loop location analysis

Preliminary

Wrangling the peaks

files <- list(
  "EUFAPAF1_rep1" = "data/drip_seq/eufapaf1_drip/EUFA_PAF_1_S30/peaks_macs_unstranded/EUFA_PAF_1_S30_hg38.unstranded_peaks.xls",
  "EUFAPAF1_rep2" = "data/drip_seq/eufapaf1_drip/EUFA_PAF_2_S31/peaks_macs_unstranded/EUFA_PAF_2_S31_hg38.unstranded_peaks.xls",
  "EUFAPAF1_rep3" = "data/drip_seq/eufapaf1_drip/EUFA_PAF_3_S32/peaks_macs_unstranded/EUFA_PAF_3_S32_hg38.unstranded_peaks.xls"
)
EUFAPAF1gr <- lapply(files, function(file) {
  tmpdf <- read_tsv(file, skip = 30) %>%
    rename_with(~ gsub(pattern = "-log10\\((.+)\\)",
                       replacement = "\\1", .x)) %>%
    as.data.frame()
  gr <- toGRanges(tmpdf)
  gr <- gr[which(width(gr) < 10000),]
})
sapply(names(EUFAPAF1gr), function(peakName){
  peaks <- EUFAPAF1gr[[peakName]]
  paste0("number of peaks: ", length(peaks$length))
})
##            EUFAPAF1_rep1            EUFAPAF1_rep2            EUFAPAF1_rep3 
## "number of peaks: 48512" "number of peaks: 67392" "number of peaks: 42934"

Overlap of replicates

# Find the overlap of the replicates
olEUFAPAF1 <- findOverlapsOfPeaks(EUFAPAF1gr)
olEUFAPAF1 <- addMetadata(olEUFAPAF1, colNames="qvalue", FUN=mean) 

Venn diagram of replicate overlap

tmp = makeVennDiagram(olEUFAPAF1, fill = c("skyblue", "firebrick", "goldenrod"),
                margin = .05)  # Decent sized overlap

Comparison of P Adjusted Value between overlapping/non-overlapping

macs2 assigns an adjusted p value to each peak, indicating the level of confidence in the peak calling. By calculating the overlap of peaks between biological replicates, we are (in theory) finding the genuine peaks. We can verify this by comparing the P adjusted values of the peaks that overlapped to the ones which didn’t (below).

olQV <- c(olEUFAPAF1$peaklist$`EUFAPAF1_rep1///EUFAPAF1_rep2`$qvalue,
          olEUFAPAF1$peaklist$`EUFAPAF1_rep2///EUFAPAF1_rep3`$qvalue,
          olEUFAPAF1$peaklist$`EUFAPAF1_rep1///EUFAPAF1_rep3`$qvalue,
          olEUFAPAF1$peaklist$`EUFAPAF1_rep1///EUFAPAF1_rep2///EUFAPAF1_rep3`$qvalue)
nonolQV <- c(olEUFAPAF1$peaklist$EUFAPAF1_rep2$qvalue,
             olEUFAPAF1$peaklist$EUFAPAF1_rep3$qvalue,
             olEUFAPAF1$peaklist$EUFAPAF1_rep1$qvalue)
tibble(
  qval = c(olQV, nonolQV),
  group = c(rep("Overlapping", length(olQV)), 
            rep("Non-overlapping", length(nonolQV)))
) %>%
  ggplot(mapping = aes(y = qval, x = group, fill = group)) +
  geom_boxplot() +
  ylab("P Adjusted Value (-log10)") +
  xlab(NULL) +
  theme_bw(base_size = 15) +
  ggpubr::rremove("legend") +
  ggpubr::stat_compare_means(comparisons = list(c("Overlapping", "Non-overlapping")),
                             label = "p.signif", size = 5) +
  scale_y_continuous(limits = c(-5, 70)) +
  labs(title = "Overlapping/Non-overlapping R-loop sites P-Adj Value")

The results indicate that the overlapping peaks are only a bit more robust measure of EUFAPAF1 R-loop locations, indicated by the fact that they are more significant. However, there are other ways to assess this…

Comparing the R-loop profile of overlapping and non-overlapping peaks

Another way to assess the robustness of our overlapping peakset is to check the binding profile using a tool like ChIPseeker. In particular, we expect to see that most genuine EUFAPAF1 R-loop locations are found closer to the TSS and enriched in genic regions.

Metaplot around TSS

plotAvgProf(tagMatrixListEUFAPAF1, xlim=c(-3000, 3000), facet = "row") +
  labs(title = "Overlapping vs Non-overlapping EUFAPAF1 R-loop peaks around TSS")

Distance to TSS plot

plotDistToTSS(peakAnnoListEUFAPAF1) + 
  labs(title = "EUFAPAF1 peak locations relative to TSS")

Feature distribution plot

plotAnnoBar(peakAnnoListEUFAPAF1) +
  labs(title = "Feature Overlap with EUFAPAF1 Peaks")

However, we find that this isn’t completely clearcut. It seems that there are genuine R-loops occuring in the non-overlapping group. Therefore, the final peakset will be considered as R-loop peaks which either:

  1. Are found in both replicates (overlapping peaks)
  2. Have a P adjusted value < the 50% quantile of all p values
cutoff <- quantile(c(olEUFAPAF1$all.peaks$EUFAPAF1_rep1$qvalue, olEUFAPAF1$all.peaks$EUFAPAF1_rep3$qvalue,
                     olEUFAPAF1$all.peaks$EUFAPAF1_rep2$qvalue))[3]
nonOlKeep <- olEUFAPAF1$uniquePeak[olEUFAPAF1$uniquePeaks$qvalue > cutoff, c("qvalue")]
EUFAPAF1peaks <- c(nonOlKeep, olEUFAPAF1$peaklist$`EUFAPAF1_rep1///EUFAPAF1_rep2`[,c("qvalue")],
                   olEUFAPAF1$peaklist$`EUFAPAF1_rep2///EUFAPAF1_rep3`[,c("qvalue")],
                   olEUFAPAF1$peaklist$`EUFAPAF1_rep1///EUFAPAF1_rep3`[,c("qvalue")],
                   olEUFAPAF1$peaklist$`EUFAPAF1_rep1///EUFAPAF1_rep2///EUFAPAF1_rep3`[,c("qvalue")])
names(EUFAPAF1peaks) <- NULL

Save EUFAPAF1 R-loop peaks

EUFAPAF1peaks <- keepStandardChromosomes(EUFAPAF1peaks, pruning.mode = "coarse")
rtracklayer::export(EUFAPAF1peaks, con = "analysis/diff_drip_brca2_paf1/results/EUFAPAF1peaks.bed")

Results

What kinds of genomic features do EUFAPAF1 R-loops appear in?

EUFAPAF1_annot = annotate_regions(
  regions = EUFAPAF1peaks,
  annotations = annotations,
  ignore.strand = TRUE,
  quiet = F)
regioneR::getMask(genome = "hg38")
## GRanges object with 0 ranges and 0 metadata columns:
##    seqnames    ranges strand
##       <Rle> <IRanges>  <Rle>
##   -------
##   seqinfo: no sequences
EUFAPAF1random <- regioneR::randomizeRegions(
  A = EUFAPAF1peaks,  
  genome = "hg38"
)
EUFAPAF1random_annot = annotate_regions(
  regions = EUFAPAF1random,
  annotations = annotations,
  ignore.strand = TRUE,
  quiet = FALSE)
plot_annotation(
  annotated_regions = EUFAPAF1_annot,  quiet = FALSE,
  annotated_random = EUFAPAF1random_annot,
  annotation_order = annotation_order,
  plot_title = 'EUFAPAF1 R-loop Site Feature Overlaps (vs Random)',
  x_label = 'Annotations',
  y_label = 'Count') 

What genes are they in? What pathways are these genes involved in?

EUFAPAF1geneAnno <- annotatePeakInBatch(EUFAPAF1peaks, output = "overlap", maxgap = 1000,
                                  AnnotationData=annoData)
EUFAPAF1geneAnno <- addGeneIDs(EUFAPAF1geneAnno,
                         "org.Hs.eg.db",
                         IDs2Add = "SYMBOL")

EUFAPAF1 R-loop genes saved to TSV file.

EUFAPAF1genes <- as.data.frame(EUFAPAF1geneAnno) %>%
  dplyr::filter(! is.na(SYMBOL)) %>%
  dplyr::select(SYMBOL) %>%
  dplyr::filter(! grepl(SYMBOL, pattern = ";|/")) %>%
  dplyr::distinct(SYMBOL) %>%
  write_tsv(file = "analysis/diff_drip_brca2_paf1/results/EUFAPAF1_bound_genes.tsv") %>%
  pull(SYMBOL) 

Pathway enrichment with Enrichr

The ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X and KEGG_2019_HUMAN databases were queried to enrich for relevant gene sets. The analysis with all possible genes sets is available permanently at this link

  1. KEGG Pathways
KEGG_EUFAPAF1 <- EUFAPAF1enriched$KEGG_2019_Human
KEGG_EUFAPAF1 %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top KEGG Pathways Enriched in EUFAPAF1 R-loop genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

  1. ChEA Pathways (Transcription Factors)
ChEA_EUFAPAF1 <- EUFAPAF1enriched$`ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X`
ChEA_EUFAPAF1 %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top ChEA Pathways Enriched in EUFAPAF1 R-loop genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

EUFA, EUFAB2, and EUFAPAF1 R-loop Overlap Analysis

# Find the overlap of the replicates
ol <- findOverlapsOfPeaks(peakList)
ol <- addMetadata(ol, colNames="qvalue", FUN=mean) 
tmp = makeVennDiagram(ol, fill = c("skyblue", "firebrick", "goldenrod"),
                margin = .05)  # Decent sized overlap

Note: These results are difficult to interpret. With DRIP-Seq it is necessary to use a differential stats model as the pure peak locations is influenced by a large number of technical factors. Considering that we have two batchs here, it isn’t possible from overlap alone to determine the impact of BRCA2 and PAF1 on R-loops.

Differential method to find how BRCA2 and PAF1 impacts R-loop formation

This approach uses the DESeq2 differential expression stats model provided by the DiffBind package. First, the consensus peaks are derived from the union of all R-loop peaks. Then, the alignment files (.bam files) are quantified across these consensus peaks to build a “peak count matrix”. Then, DESeq2 is used to find the differential R-loops between conditions. In this case, we are searching for the R-loops which decrease with EUFAB2 or EUFAPAF1 compared to EUFA alone

PCA showing the difference between EUFAB2, EUFAPAF1, and EUFA

load("analysis/diff_drip_brca2_paf1/dbadata.rda")
pcaplot

EUFAB2 vs EUFA Differential R-loop Loss (i.e., “What R-loops does BRCA2 degrade?”)

MA Plot showing the effect of BRCA2 on R-loop abundance

maplot1

Volcano Plot showing the effect of BRCA2 on R-loop abundance

volcano1

dbgr <- dbgr1
dbgrsigup <- dbgr[dbgr$FDR < .05 & dbgr$Fold > 0,]
dbgrsigdown <- dbgr[dbgr$FDR < .05 & dbgr$Fold < 0,]
rtracklayer::export(dbgrsigdown, con = "analysis/diff_drip_brca2_paf1/results/diffDRIP_EUFAB2_peaks.bed")

Where are the sites that show decreased R-loops with BRCA2 (B2vsEUFA DA R-loops)? What features do they overlap with? How does this compare with typical R-loop peaks and the sites of increased R-loops?

Metaplot around TSS

plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), facet = "row") +
  labs(title = "B2vsEUFA DA R-loops around TSS")

Distance to TSS plot

plotDistToTSS(peakAnnoList) + 
  labs(title = "B2vsEUFA DA R-loop locations relative to TSS")

Feature distribution plot

plotAnnoBar(peakAnnoList) +
  labs(title = "Feature Overlap with B2vsEUFA DA R-loops")

peaks <- dbgrsigdown
peaks$FDR <- -log10(peaks$FDR)
names(peaks) <- NULL

Save differentially decreased (BRCA2-degraded) R-loop peaks

peaks <- keepStandardChromosomes(peaks, pruning.mode = "coarse")
rtracklayer::export(peaks, con = "analysis/diff_drip_brca2_paf1/results/EUFAB2vsEUFA_degraded_RLoops.bed")

What kinds of genomic features do these BRCA2-degraded R-loops bind to?

plot_annotation(
  annotated_regions = annot, quiet = FALSE,
  annotated_random = random_annot,
  annotation_order = annotation_order,
  plot_title = 'EUFAB2vsEUFA Degraded R-loops Feature Overlaps (vs Random)',
  x_label = 'Annotations',
  y_label = 'Count') 

What genes do these EUFAB2vsEUFA Degraded R-loops occur in/near? What pathways do they relate to?

geneAnno <- annotatePeakInBatch(peaks, output = "overlap", maxgap = 1000, 
                                  AnnotationData=annoData)
geneAnno <- addGeneIDs(geneAnno,
                         "org.Hs.eg.db",
                         IDs2Add = "SYMBOL")

BRCA2-degraded R-loop genes saved to TSV file.

genes <- as.data.frame(geneAnno) %>%
  dplyr::filter(! is.na(SYMBOL)) %>%
  dplyr::select(SYMBOL) %>%
  dplyr::distinct(SYMBOL) %>%
  dplyr::filter(! grepl(SYMBOL, pattern = ";")) %>%
  write_tsv(file = "analysis/diff_drip_brca2_paf1/results/diffDRIP_B2vsEUFA_genes.tsv") %>%
  pull(SYMBOL) 
as.data.frame(geneAnno) %>%
  dplyr::filter(! is.na(SYMBOL)) %>%
  arrange(desc(FDR)) %>%
  dplyr::select(c(12, 1, 2, 3, 9, 10, 11, 21)) %>%
  dplyr::mutate(FDR = 10^(-1*FDR)) %>%
  DT::datatable(rownames = FALSE)

Pathway enrichment with Enrichr

The ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X and KEGG_2019_HUMAN databases were queried to enrich for relevant gene sets. The analysis with all possible genes sets is available permanently at this link

  1. KEGG Pathways
KEGG_ <- enriched$KEGG_2019_Human
KEGG_ %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top KEGG Pathways Enriched in B2vsEUFA DA R-loops (down)",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

  1. ChEA Pathways (Transcription Factors)
ChEA_ <- enriched$`ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X`
ChEA_ %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top ChEA Pathways Enriched in B2vsEUFA DA R-loops (down)",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

EUFAPAF1 vs EUFA

MA Plot showing the effect of PAF1 on R-loop abundance

maplot2

Volcano Plot showing the effect of PAF1 on R-loop abundance

volcano2

dbgr <- dbgr2
dbgrsigup <- dbgr[dbgr$FDR < .05 & dbgr$Fold > 0,]
dbgrsigdown <- dbgr[dbgr$FDR < .05 & dbgr$Fold < 0,]
rtracklayer::export(dbgrsigdown, con = "analysis/diff_drip_brca2_paf1/results/diffDRIP_EUFAPAF1_peaks.bed")

Where are the sites that show decreased R-loops with PAF1 (PAF1vsEUFA degraded R-loops)? What features do they overlap with?

Metaplot around TSS

plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), facet = "row") +
  labs(title = "Comparison of R-loops around TSS")

Distance to TSS plot

plotDistToTSS(peakAnnoList) + 
  labs(title = "Comparison of R-loops locations relative to TSS")

Feature distribution plot

plotAnnoBar(peakAnnoList) +
  labs(title = "Feature Overlap")

peaks <- dbgrsigdown
peaks$FDR <- -log10(peaks$FDR)
names(peaks) <- NULL

Save peaks

peaks <- keepStandardChromosomes(peaks, pruning.mode = "coarse")

What kinds of genomic features does it bind to?

plot_annotation(
  annotated_regions = annot, quiet = FALSE,
  annotated_random = random_annot,
  annotation_order = annotation_order,
  plot_title = 'PAF1vsEUFA degraded R-loops Feature Overlaps (vs Random)',
  x_label = 'Annotations',
  y_label = 'Count') 

What genes do these PAF1vsEUFA degraded R-loops occur in/near? What pathways do they relate to?

geneAnno <- annotatePeakInBatch(peaks, output = "overlap", maxgap = 1000, 
                                  AnnotationData=annoData)
geneAnno <- addGeneIDs(geneAnno,
                         "org.Hs.eg.db",
                         IDs2Add = "SYMBOL")

PAF1-degraded R-loop genes saved to TSV file.

genes <- as.data.frame(geneAnno) %>%
  dplyr::filter(! is.na(SYMBOL)) %>%
  dplyr::select(SYMBOL) %>%
  dplyr::distinct(SYMBOL) %>%
  dplyr::filter(! grepl(SYMBOL, pattern = ";")) %>%
  write_tsv(file = "analysis/diff_drip_brca2_paf1/results/diffDRIP_PAF1vsEUFA_genes.tsv") %>%
  pull(SYMBOL) 
as.data.frame(geneAnno) %>%
  dplyr::filter(! is.na(SYMBOL)) %>%
  arrange(desc(FDR)) %>%
  dplyr::select(c(12, 1, 2, 3, 9, 10, 11, 21)) %>%
  dplyr::mutate(FDR = 10^(-1*FDR)) %>%
  DT::datatable(rownames = FALSE)
## Warning in instance$preRenderHook(instance): It seems your data is too big
## for client-side DataTables. You may consider server-side processing: https://
## rstudio.github.io/DT/server.html

Pathway enrichment with Enrichr

The ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X and KEGG_2019_HUMAN databases were queried to enrich for relevant gene sets. The analysis with all possible genes sets is available permanently at this link

  1. KEGG Pathways
KEGG_ <- enriched$KEGG_2019_Human
KEGG_ %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top KEGG Pathways Enriched in PAF1vsEUFA degraded R-loops",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

  1. ChEA Pathways (Transcription Factors)
ChEA_ <- enriched$`ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X`
ChEA_ %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top ChEA Pathways Enriched in PAF1vsEUFA degraded R-loops",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

Compare the B2-degraded R-loops with the PAF1-degraded R-loops

# Find the overlap of the replicates
diffgrlist <- list(
  "BRCA2-degraded R-loops" =dbgr1[dbgr1$FDR < .05 & dbgr1$Fold < 0,],
  "PAF1-degraded R-loops" = dbgr2[dbgr2$FDR < .05 & dbgr2$Fold < 0,]
)
diffgrlist <- lapply(diffgrlist, function(x) {
  x$FDR <- -log10(x$FDR)
  x
})
oldiffgrlist <- findOverlapsOfPeaks(diffgrlist)
oldiffgrlist <- addMetadata(oldiffgrlist, colNames="FDR", FUN=mean) 

Venn diagram of replicate overlap

tmp = makeVennDiagram(oldiffgrlist, fill = c("skyblue", "firebrick"),
                margin = .1)  # Decent sized overlap

Comparison of P Adjusted Value between overlapping/non-overlapping

macs2 assigns an adjusted p value to each peak, indicating the level of confidence in the peak calling. By calculating the overlap of peaks between biological replicates, we are (in theory) finding the genuine peaks. We can verify this by comparing the P adjusted values of the peaks that overlapped to the ones which didn’t (below).

olQV <- c(oldiffgrlist$peaklist$`BRCA2-degraded R-loops///PAF1-degraded R-loops`$FDR)
nonolQV <- c(oldiffgrlist$peaklist$`PAF1-degraded R-loops`$FDR,
             oldiffgrlist$peaklist$`BRCA2-degraded R-loops`$FDR)
tibble(
  qval = c(olQV, nonolQV),
  group = c(rep("Overlapping", length(olQV)), 
            rep("Non-overlapping", length(nonolQV)))
) %>%
  ggplot(mapping = aes(y = qval, x = group, fill = group)) +
  geom_boxplot() +
  ylab("P Adjusted Value (-log10)") +
  xlab(NULL) +
  theme_bw(base_size = 15) +
  ggpubr::rremove("legend") +
  ggpubr::stat_compare_means(comparisons = list(c("Overlapping", "Non-overlapping")),
                             label = "p.signif", size = 5) +
  scale_y_continuous(limits = c(-5, 140)) +
  labs(title = "Overlapping/Non-overlapping R-loop sites P-Adj Value")

Comparing the R-loop profile of overlapping and non-overlapping peaks

Another way to assess the robustness of our overlapping peakset is to check the binding profile using a tool like ChIPseeker. In particular, we expect to see that most genuine diffbind R-loop locations are found closer to the TSS and enriched in genic regions.

Metaplot around TSS

plotAvgProf(tagMatrixListdiffbind, xlim=c(-3000, 3000), facet = "row") +
  labs(title = "Overlapping vs Non-overlapping differential R-loop peaks around TSS")

Distance to TSS plot

plotDistToTSS(peakAnnoListdiffbind) + 
  labs(title = "differential peak locations relative to TSS")

Feature distribution plot

plotAnnoBar(peakAnnoListdiffbind) +
  labs(title = "Feature Overlap with differential Peaks")

Save diffbind R-loop peaks

diffbindpeaks <- oldiffgrlist$peaklist$`BRCA2-degraded R-loops///PAF1-degraded R-loops`
paf1diffpeaks <- oldiffgrlist$peaklist$`PAF1-degraded R-loops`
diffbindpeaks <- keepStandardChromosomes(diffbindpeaks, pruning.mode = "coarse")
rtracklayer::export(diffbindpeaks, con = "analysis/diff_drip_brca2_paf1/results/BRCA2_PAF1_shared_RLoops_lost.bed")

Results

What kinds of genomic features do they appear in?

diffbind_annot = annotate_regions(
  regions = diffbindpeaks,
  annotations = annotations,
  ignore.strand = TRUE,
  quiet = F)
regioneR::getMask(genome = "hg38")
## GRanges object with 0 ranges and 0 metadata columns:
##    seqnames    ranges strand
##       <Rle> <IRanges>  <Rle>
##   -------
##   seqinfo: no sequences
diffbindrandom <- regioneR::randomizeRegions(
  A = diffbindpeaks,  
  genome = "hg38"
)
diffbindrandom_annot = annotate_regions(
  regions = diffbindrandom,
  annotations = annotations,
  ignore.strand = TRUE,
  quiet = FALSE)
plot_annotation(
  annotated_regions = diffbind_annot,  quiet = FALSE,
  annotated_random = diffbindrandom_annot,
  annotation_order = annotation_order,
  plot_title = 'diffbind R-loop Site Feature Overlaps (vs Random)',
  x_label = 'Annotations',
  y_label = 'Count') 

What genes are they in? What pathways are these genes involved in?

diffbindgeneAnno <- annotatePeakInBatch(diffbindpeaks, output = "overlap", maxgap = 1000,
                                  AnnotationData=annoData)
diffbindgeneAnno <- addGeneIDs(diffbindgeneAnno,
                         "org.Hs.eg.db",
                         IDs2Add = "SYMBOL")

diffbind R-loop genes saved to TSV file.

diffbindgenes <- as.data.frame(diffbindgeneAnno) %>%
  dplyr::filter(! is.na(SYMBOL)) %>%
  dplyr::select(SYMBOL) %>%
  dplyr::filter(! grepl(SYMBOL, pattern = ";|/")) %>%
  dplyr::distinct(SYMBOL) %>%
  write_tsv(file = "analysis/diff_drip_brca2_paf1/results/diffbind_bound_genes.tsv") %>%
  pull(SYMBOL) 

Pathway enrichment with Enrichr

The ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X and KEGG_2019_HUMAN databases were queried to enrich for relevant gene sets. The analysis with all possible genes sets is available permanently at this link

  1. KEGG Pathways
KEGG_diffbind <- diffbindenriched$KEGG_2019_Human
KEGG_diffbind %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top KEGG Pathways Enriched in diffbind R-loop genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

  1. ChEA Pathways (Transcription Factors)
ChEA_diffbind <- diffbindenriched$`ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X`
ChEA_diffbind %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top ChEA Pathways Enriched in diffbind R-loop genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

Check the PAF1-only (non-overlapping) degraded R-loops

paf1diffpeaks <- keepStandardChromosomes(paf1diffpeaks, pruning.mode = "coarse")
paf1diffgeneAnno <- annotatePeakInBatch(paf1diffpeaks, output = "overlap", maxgap = 1000,
                                  AnnotationData=annoData)
paf1diffgeneAnno <- addGeneIDs(paf1diffgeneAnno,
                         "org.Hs.eg.db",
                         IDs2Add = "SYMBOL")

diffbind R-loop genes saved to TSV file.

paf1diffgenes <- as.data.frame(paf1diffgeneAnno) %>%
  dplyr::filter(! is.na(SYMBOL)) %>%
  dplyr::select(SYMBOL) %>%
  dplyr::filter(! grepl(SYMBOL, pattern = ";|/")) %>%
  dplyr::distinct(SYMBOL) %>%
  write_tsv(file = "analysis/diff_drip_brca2_paf1/results/diffbind_paf1_only_bound_genes.tsv") %>%
  pull(SYMBOL) 

Pathway enrichment with Enrichr

The ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X and KEGG_2019_HUMAN databases were queried to enrich for relevant gene sets. The analysis with all possible genes sets is available permanently at this link

  1. KEGG Pathways
KEGG_diffbind <- paf1enriched$KEGG_2019_Human
KEGG_diffbind %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top KEGG Pathways Enriched in diffbind R-loop genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

  1. ChEA Pathways (Transcription Factors)
ChEA_diffbind <- paf1enriched$`ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X`
ChEA_diffbind %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top ChEA Pathways Enriched in diffbind R-loop genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

Compare with BRCA2, PAF1, and XRN2 binding sites (along with siXRN2 DRIP-Seq)

EUFAB2/PAF1 vs EUFA lost R-loops vs BRCA2/PAF1 ChIP

BRCA2peaks <- import("analysis/brca2_paf1_xrn2_binding/results/B2peaks.bed")
PAF1peaks <- import("analysis/brca2_paf1_xrn2_binding/results/PAF1peaks.bed")
XRN2peaks <- import("analysis/brca2_paf1_xrn2_binding/results/XRN2peaks.bed")
listOL <- list(
  "BRCA2-degraded R-loops" = diffgrlist$`BRCA2-degraded R-loops`,
  "PAF1-degraded R-loops" = diffgrlist$`PAF1-degraded R-loops`,
  "BRCA2 ChIP" = BRCA2peaks,
  "PAF1 ChIP" = PAF1peaks,
  "XRN2 ChIP" = XRN2peaks
)
olbpr <- findOverlapsOfPeaks(listOL)
tmp = makeVennDiagram(olbpr, fill = c("firebrick", "skyblue", "forestgreen", "goldenrod", "pink"), margin = .1)

This result likely indicates the uncertainty in our measurements and suggests that the overlap between BRCA2 and PAF1 is probably larger than we are able to measure currently.

EUFAB2/PAF1 vs EUFA lost R-loops vs BRCA2/PAF1/XRN2 (all overlapping)

listOL2 <- list(
  "BRCA2 ChIP" = BRCA2peaks,
  "PAF1 ChIP" = PAF1peaks,
  "XRN2 ChIP" = XRN2peaks
)
olbpr2 <- findOverlapsOfPeaks(listOL2)
XRN2degraded <- import("analysis/diff_drip_xrn2/results/diffDRIP_siXRN2_peaks.bed")
listOL3 <- list(
  "BRCA2-degraded R-loops" = diffgrlist$`BRCA2-degraded R-loops`,
  "PAF1-degraded R-loops" = diffgrlist$`PAF1-degraded R-loops`,
  "XRN2-degraded R-loops" = XRN2degraded,
  "BRCA2-PAF1-XRN2 sites" = c(olbpr2$peaklist$`BRCA2 ChIP///XRN2 ChIP`, olbpr2$peaklist$`PAF1 ChIP///XRN2 ChIP`,
                              olbpr2$peaklist$`BRCA2 ChIP///PAF1 ChIP`, olbpr2$peaklist$`BRCA2 ChIP///PAF1 ChIP///XRN2 ChIP`)
)
olbpr3 <- findOverlapsOfPeaks(listOL3)
tmp = makeVennDiagram(olbpr3, fill = c("firebrick", "skyblue", "goldenrod", "forestgreen"), margin = .05)

Enrichr for everything in the overlap

fullOL <- olbpr3$peaklist$`BRCA2-degraded R-loops///PAF1-degraded R-loops///XRN2-degraded R-loops///BRCA2-PAF1-XRN2 sites`
fullOL <- keepStandardChromosomes(fullOL, pruning.mode = "coarse")

geneAnno <- annotatePeakInBatch(myPeakList = fullOL, output = "overlap", maxgap = 1000,
                                  AnnotationData=annoData)
geneAnno <- addGeneIDs(geneAnno,
                         "org.Hs.eg.db",
                         IDs2Add = "SYMBOL")

Overlap genes saved to TSV file.

genes <- as.data.frame(geneAnno) %>%
  dplyr::filter(! is.na(SYMBOL)) %>%
  dplyr::select(SYMBOL) %>%
  dplyr::distinct(SYMBOL) %>%
  dplyr::filter(! grepl(SYMBOL, pattern = ";|/")) %>%
  write_tsv(file = "analysis/diff_drip_brca2_paf1/results/genes_in_all_overlap.tsv") %>%
  pull(SYMBOL) 

Pathway enrichment with Enrichr

The ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X and KEGG_2019_HUMAN databases were queried to enrich for relevant gene sets. The analysis with all possible genes sets is available permanently at this link

  1. KEGG Pathways
KEGG_ <- enriched$KEGG_2019_Human
KEGG_ %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top KEGG Pathways Enriched in DA R-loops and XRN2/BRCA2/PAF1-bound genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

  1. ChEA Pathways (Transcription Factors)
ChEA_ <- enriched$`ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X`
ChEA_ %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top ChEA Pathways Enriched in DA R-loops and XRN2/BRCA2/PAF1-bound genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

Enrichr for everything that has at least three overlaps

fullOL2 <- unlist(GRangesList(olbpr3$peaklist[grep(names(olbpr3$peaklist), pattern = ".+///.+///.+")]))
fullOL2 <- keepStandardChromosomes(fullOL2, pruning.mode = "coarse")

geneAnno <- annotatePeakInBatch(myPeakList = fullOL2, output = "overlap", maxgap = 1000,
                                  AnnotationData=annoData)
## Warning in annotatePeakInBatch(myPeakList = fullOL2, output = "overlap", : Found duplicated names in myPeakList. 
##                     Changing the peak names ...
geneAnno <- addGeneIDs(geneAnno,
                         "org.Hs.eg.db",
                         IDs2Add = "SYMBOL")

Overlap genes saved to TSV file.

genes2 <- as.data.frame(geneAnno) %>%
  dplyr::filter(! is.na(SYMBOL)) %>%
  dplyr::select(SYMBOL) %>%
  dplyr::distinct(SYMBOL) %>%
  dplyr::filter(! grepl(SYMBOL, pattern = ";|/")) %>%
  write_tsv(file = "analysis/diff_drip_brca2_paf1/results/genes_in_atleast_3_overlap.tsv") %>%
  pull(SYMBOL) 

Pathway enrichment with Enrichr

The ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X and KEGG_2019_HUMAN databases were queried to enrich for relevant gene sets. The analysis with all possible genes sets is available permanently at this link

  1. KEGG Pathways
KEGG_ <- enriched$KEGG_2019_Human
KEGG_ %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top KEGG Pathways Enriched in DA R-loops and XRN2/BRCA2/PAF1-bound genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

  1. ChEA Pathways (Transcription Factors)
ChEA_ <- enriched$`ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X`
ChEA_ %>%
  top_n(10, Combined.Score) %>% 
  arrange(desc(Combined.Score)) %>%
  mutate(Term = factor(Term, levels = rev(Term))) %>%
  ggplot(aes(x = Term, y = Combined.Score, fill = Adjusted.P.value)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 14) +
  xlab(NULL) +
  labs(title = "Top ChEA Pathways Enriched in DA R-loops and XRN2/BRCA2/PAF1-bound genes",
       fill = "Adjusted P Value") +
  ylab("Combined Score") +
  ggpubr::rotate() 

Discussion

Questions to answer

  1. How does BRCA2 and PAF1 status impact R-loops? Are there regions of the genome which are particularly impacted?
  2. What genes are R-loops differentially found in? What pathways do they relate to?
  3. How do these results relate to BRCA2, PAF1, and XRN2 binding sites? How do they relate to the siXRN2 DRIP-Seq results?

How does BRCA2 and PAF1 status impact R-loops? Are there regions of the genome which are particularly impacted?

Comparing the binding sites directly seemed to produce inconclusive results. Therefore, we implemented a differential binding model from the DiffBind package. This revealed 3,501 R-loop sites that are differential between EUFAB2 and EUFA. Of these, the vast majority (2,839 R-loop sites) were lost with BRCA2 complementation. It also revealed 19534 R-loops that are differentially abundance between EUFAPAF1 and EUFA. Of these, the vast majority (16,485) are lost with PAF1 complementation. Interestingly the majority of R-loops lost with BRCA2 complementation were also lost with PAF1 complementation. Considering the difference in quality between the BRCA2 and PAF1 datasets, this indicates that BRCA2 and PAF1 probably help to degrade R-loops at the same places in a much greater number of sites than we are capable of measuring, but that PAF1 retains some independence as well.

What genes are R-loops differentially found in? What pathways do they relate to?

From analysis of the DA R-loops, it was found that which showed evidence of being degraded by both BRCA2 and PAF1 are involved in a variety of processes. From the analysis with enrichr, we see the PRDM5, HIF1A, UBTF, GATA6, TAF1, BRCA1 and other TFs we had previously noticed. We also noticed pathways related to mRNA processing, translation, estrogen treatment, NF-kB, EMT, and cell cycle. This indicates that the R-loops degraded by PAF1 and BRCA2 may be involved in processes related to maintanence of mesodermal lineage, ribosome biogenesis, transcriptional regulation, and control of cell cycle. This fits in with the theory that R-loops can play a role in the determination of cellular state in mesodermal tissues like fibroblasts.

How do these results relate to BRCA2, PAF1, and XRN2 binding sites? These results were compared with the ChIP-Seq for BRCA2, PAF1, and XRN2. The overlap was relatively weak (46 sites) and the gene enrichment with enrichr showed weak enrichment for ribosome and TFs we had previously identified. However, the amount of noise in the data is basically intractible by the time we begin comparing so many different studies and modalities, so it was reasoned that a more accurate representation of the true overlap in EUFA cells would be achieved by finding the places where at least 3 of the 4 groups (B2, PAF1, XRN2 DRIP and the intersected set of ChIP sites) overlapped. With that constraint relaxed, we found 1,288 genes that show the overlap of XRN2, BRCA2, and PAF1 ChIP and DRIP. By using enrichr we found that these genes are involved in pathways related to XRN2, BRCA2, MYC, ATF2, TAF1, ZBTB7A, E2F6 and other transcription factors we had previously identified. We also found pathways related to ribosome biogenesis, mRNA processing, cell cycle, mesodermal differentiation, breast and cervical cancer, and estrogen treatment.

Conclusions and future directions

From this analysis, we learned BRCA2 and PAF1 complementation have a largely overlapping effect in reducing R-loops in EUFA cells. We also learning that this happens a genes related to biology we have already uncovered from previous analyses, such as mesodermal differentiation, cell cycle, ribosome biogenesis, mRNA processing, NF-KB, ATFs, UBTF, XRN2, BRCA1, E2F6, and others. By comparing with the results of the siXRN2 study and with ChIP-Seq for BRCA2, XRN2, and PAF1, we found further evidence that BRCA2, PAF1, and XRN2 are co-localized at proximal promoter regions and degrading R-loops in a coordinated fashion. However, due to an overabundance of noise in the data, we’ve reached the limit of our ability to gain further insight into where they are overlapping or about the dynamics of their interactions with R-loops. It will remain for future experiments, such as the ATAC-Seq and (eventually) GRO-Seq to add additional resolution to our understanding of the interaction of these factors.