knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file(),
cache = FALSE)
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.
From GitHub:
Questions:
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.
# 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")
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"
# 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
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…
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:
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")
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')
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
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()
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)
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()
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()
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"
# 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
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…
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:
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")
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')
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
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()
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()
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"
# 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
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…
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:
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")
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')
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
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()
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()
# 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.
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
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")
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")
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')
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
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()
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()
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")
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")
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')
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
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()
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()
# 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
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")
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")
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')
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
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()
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()
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
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()
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()
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.
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)
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
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()
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()
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
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()
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()
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.
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.