• Parameter values
  • Load required packages and helper functions
  • Read data
    • IP-MS
    • ChIP
  • Helper functions
  • Extract data for IP-MS heatmaps
  • Calculate peak distances to tRNA or rRNA genes
  • Facts and numbers - ChIP-seq heatmap
  • Figure 2
    • Settings
    • Summary heatmaps using results from the statistical tests (150 mM NaCl)
      • Overall heatmap
    • Summary heatmaps using results from the statistical tests (500 mM NaCl)
      • Overall heatmap
    • Supplementary plots - IP-MS data for HOT TFs
    • Peak vs. sample heatmap (binary)
    • Peak vs. sample heatmap (enrichments)
    • Supplementary plots - ChIP enrichment replicate-pair correlations
    • Supplementary plots - Bound promoters
    • Supplementary plots - PCA (IP-MS)
    • Put together
    • Supplementary figure
  • Session info

Parameter values

params
## $rds150
## [1] "data/ipms_150_sce.rds"
## 
## $rds500
## [1] "data/ipms_500_sce.rds"
## 
## $complexes
## [1] "data/complexes.json"
## 
## $idmap
## [1] "data/id_mapping_table.txt"
## 
## $baitclass
## [1] "data/ipms_bait_class.txt"
## 
## $dbdtxt
## [1] "data/TF_DBD_annotation.txt"
## 
## $peakcsv
## [1] "data/fused_peaks_filtered.csv.gz"
## 
## $peakenr
## [1] "data/fused_peaks_filtered_enrichments.csv.gz"
## 
## $peakatacfile
## [1] "data/fused_peaks_filtered_counts-external_atac.csv.gz"
## 
## $peakchipfile
## [1] "data/fused_peaks_filtered_counts-external_chip.csv.gz"
## 
## $rdsrnaseq_del
## [1] "data/rnaseq_del.rds"
## 
## $txcsv
## [1] "data/promoters-1kb_tes-1kb_annotated.csv.gz"
## 
## $gtf
## [1] "reference/Schizosaccharomyces_pombe.ASM294v2.55.gtf"
## 
## $ipms150untag
## [1] "data/ipms_150_untag_sce.rds"

Load required packages and helper functions

suppressPackageStartupMessages({
    library(SingleCellExperiment)
    library(ggplot2)
    library(cowplot)
    library(dplyr)
    library(forcats)
    library(jsonlite)
    library(ComplexHeatmap)
    library(colorspace)
    library(circlize)
    library(scattermore)
    library(tidyr)
    library(einprot)
})

## Source scripts with required helper functions and settings
source("params/plot_settings.R")
source("params/get_testres_function.R")
source("params/ipms_params.R")
source("params/mapping_functions.R")
source("params/ipms_heatmap_functions.R")

## Heatmap sizes
w_ipmsHm <-  1.2 * 84 # heatmap body width (mm)
h_ipmsHm <- 178 # total heatmap height (mm)
w_peaksHm <-  84 # heatmap body width (mm)
h_peaksHm <- 178 # total heatmap height (mm)

fs <- 8  # small font size
fm <- 9  # medium font size
fl <- 10 # large font size
ft <- 12 # title font size

Read data

IP-MS

We load two SingleCellExperiment objects, containing data and results from the low- and high-salt experiments, respectively.

sce150 <- readRDS(params$rds150)
sce500 <- readRDS(params$rds500)

We then extract the full list of baits, as well as a classification according to the experiment(s) where a certain bait was pulled down.

idmap <- read.delim(params$idmap)
baitclass <- read.delim(params$baitclass)

In the following set of IPs, we identified the bait as not pulled down.

(nobait_150 <- baitclass$Gene_name[baitclass$class == "nobait"])
## [1] "Cha4"        "Cuf2"        "Atf21"       "Cbf12"       "Sre2"       
## [6] "SPBC1348.12" "Untagged"

We also read a list of known complexes, which will be used for annotation purposes.

complexes <- jsonlite::read_json(params$complexes, simplifyVector = TRUE)
complexes$`MBF transcription complex`$color <- main_colors[5]
complexes$`CCAAT-binding factor complex`$color <- main_colors[1]
complexes$`Atf1-Pcr1`$color <- main_colors[3]
names(complexes)
##  [1] "NuA4"                            "SAGA"                           
##  [3] "19S proteasome - base"           "Clr6 HDAC - complex I'"         
##  [5] "SREBP-SCAP"                      "SPBC16G5.16-SPBC530.08"         
##  [7] "Atf1-Pcr1"                       "chaperonin-containing T-complex"
##  [9] "CCAAT-binding factor complex"    "Dal2"                           
## [11] "Cyc8(Ssn6)-Tup1"                 "Rad24-Rad25"                    
## [13] "Ste11-Cdc2-Cig2"                 "Pap1-Prr1"                      
## [15] "MBF transcription complex"

Finally we load a table classifying each TF to its DNA-binding domain family.

dbd <- read.delim(params$dbdtxt) |>
    dplyr::mutate(Gene_name = .capitalize(Gene_name))

ChIP

We start by reading all tables and objects needed for the ChIP part of the figure. Count tables are in addition normalized, and we also calculate averages over replicates or all samples where needed:

# `peakgr`: ChIP-seq peaks (as a GRanges object)
peakgr <- as(read.csv(params$peakcsv, row.names = 1), "GRanges")
peakgr$peaktype <- factor(peakgr$peaktype,
                          levels = c("common peaks (ubiquitous)",
                                     "common peaks (frequent)",
                                     "specific peaks"))

# `peakenr`: peak IP-enrichments log2(IP/input)
peakenr <- as.matrix(read.csv(params$peakenr, row.names = 1))


# ... average replicates (`peakenrAvg`)
grp <- sub("_rep[12]$", "", colnames(peakenr))
peakenrAvg <- do.call(cbind, lapply(split(colnames(peakenr), grp)[unique(grp)],
                                    function(j) {
                                        rowMeans(peakenr[, j, drop = FALSE])
                                    }))

# `rnaseq_del`: a `SummarizedExperiment` with RNA-seq quantifications
rnaseq_del <- readRDS(params$rdsrnaseq_del)

# `txannot`: data.frame with transcript start site (TSS) and
#            transcript end site (TES) annotation
txannot <- read.csv(params$txcsv, row.names = 1)

# `atac_raw`: public ATAC-seq read counts in ChIP-seq peaks from this study
atac_raw <- as.matrix(read.csv(params$peakatacfile, row.names = 1))

# ... calculate reads per million (`atac_cpm`) and
#     per kilobase and million (`atac_rpkm`)
#     and average experiments (`atac_cpm_avg`, `atac_rpkm_avg`)
atac_cpm <- sweep(x = atac_raw[, -1], MARGIN = 2,
                  STATS = colSums(atac_raw[, -1]), FUN = "/") * 1e6
atac_rpkm <- atac_cpm / atac_raw[, "width"] * 1e3
atac_cpm_avg <- rowMeans(atac_cpm)
atac_rpkm_avg <- rowMeans(atac_rpkm)

# `chip_raw`: public ChIP-seq read counts in ChIP-seq peaks from this study
chip_raw <- as.matrix(read.csv(params$peakchipfile, row.names = 1))

# ... calculate reads per million (`chip_cpm`) and
#     per kilobase and million (`chip_rpkm`)
#     and calculate enrichments (`chip_enr`, log2 IP/input) averaged over replicates
chip_cpm <- sweep(x = chip_raw[, -1], MARGIN = 2,
                  STATS = colSums(chip_raw[, -1]), FUN = "/") * 1e6
chip_rpkm <- chip_cpm / chip_raw[, "width"] * 1e3
chip_series <- sub("^ChIP_([^_]+)_.+$", "\\1", colnames(chip_cpm))
chip_enr <- do.call(
    cbind,
    lapply(split(colnames(chip_cpm), chip_series),
           function(nms) {
               grp <- sub("ChIP_(GSE[0-9]+)_GSM[0-9]+_(.+)_rep[0-9]$",
                          "\\1_\\2", nms)
               stopifnot(any(is_input <- grep("_Input", grp)))
               enr <- do.call(
                   cbind,
                   lapply(unique(grp[-is_input]), function(grp1) {
                       rowMeans(log2(chip_cpm[, nms[grp == grp1], drop = FALSE] + 1)) -
                           rowMeans(log2(chip_cpm[, nms[is_input], drop = FALSE] + 1))
                   }))
               colnames(enr) <- unique(grp[-is_input])
               enr
           }))


# `genegr`: chromosomal ranges of genes (as GRanges object)
genegr <- rtracklayer::import(params$gtf, feature.type = "gene")

Helper functions

We define a number of additional helper functions for data extraction and summarization, which will be used below.

# for each region in `from`, calculate the distance to its nearest element in `to`
# if there are no element in `to` on the sequence of `from`, set distance to NA
calcDistanceToNearest <- function(from, to) {
    stopifnot(exprs = {
        is(from, "GRanges")
        is(to, "GRanges")
    })
    tmp <- distanceToNearest(x = from, subject = to)
    d <- rep(NA, length(from))
    d[queryHits(tmp)] <- mcols(tmp)$distance
    return(d)
}

Extract data for IP-MS heatmaps

hmdata_150 <- makeHeatmapData(sce = sce150, adjpthr = adjpThreshold, 
                              log2fcthr = log2fcThreshold, conc = 150,
                              idmap = idmap, baitclass = baitclass)
hmdata_500 <- makeHeatmapData(sce = sce500, adjpthr = adjpThreshold, 
                              log2fcthr = log2fcThreshold, conc = 500,
                              idmap = idmap, baitclass = baitclass)

Calculate peak distances to tRNA or rRNA genes

We want to identify ChIP-seq peaks that are near a tRNA or rRNA gene. We can obtain the coordinates of these genes from genegr and the coordinates of peaks from peakgr, and calculate the distance between nearest peak-gene pairs using distanceToNearest(). Any pair with a distance below maxdist will be classified as “near” in the plots below.

Remark: Nearest distances can be NA in cases where for example a peak resides on a sequence (chromosome) that does not contain any tRNA or rRNA gene, or vice versa, such as the sequence AB325691 (which contains gap-filling sequence between SPBPB21E7.09 and SPBPB10D8.01 in chromosome 2) or the mitochondrial sequence MT.

# distance less than `maxdist` are defined as "near"
maxdist <- 100

# tRNA
# ... PomBase and Ensembl_Fungi annotate 196 and 183 tRNAs, respectively
#     we combine the two and obtain 198 annotated tRNAs
table(genegr$gene_biotype == "tRNA", genegr$source)
##        
##         PomBase Ensembl_Fungi
##   FALSE    6818            71
##   TRUE      196           183
is_tRNA_pombase <- genegr$source == "PomBase" & genegr$gene_biotype == "tRNA"
is_tRNA_ensembl <- genegr$source == "Ensembl_Fungi" & genegr$gene_biotype == "tRNA"
is_tRNA <- is_tRNA_pombase | (is_tRNA_ensembl & !genegr %in% genegr[is_tRNA_pombase])
sum(is_tRNA)
## [1] 198
# ... now we measure distances between nearest pairs
#     either from peak to nearest tRNA
dist.peak2tRNA <- calcDistanceToNearest(from = peakgr, to = genegr[is_tRNA])
#     or from tRNA to nearest peak
dist.tRNA2peak <- calcDistanceToNearest(from = genegr[is_tRNA], to = peakgr)

# rRNA
# ... all 36 5S_rRNA genes in our annotation stem from Ensembl_Fungi
table(genegr$gene_name == "5S_rRNA", genegr$source)
##        
##         PomBase Ensembl_Fungi
##   FALSE    4524           218
##   TRUE        0            36
is_rRNA <- !is.na(genegr$gene_name) & genegr$gene_name == "5S_rRNA"
sum(is_rRNA)
## [1] 36
# ... now we measure distances between nearest pairs
#     either from peak to nearest rRNA
dist.peak2rRNA <- calcDistanceToNearest(from = peakgr, to = genegr[is_rRNA])
#     or from rRNA to nearest peak
dist.rRNA2peak <- calcDistanceToNearest(from = genegr[is_rRNA], to = peakgr)

Facts and numbers - ChIP-seq heatmap

# percent of peaks in low H3 regions (log2 IP/input < 1.0)
usePeaks <- peakgr$peaktype != "common peaks (ubiquitous)"
cat(sprintf("Percent of binding sites in regions of low H3 ChIP-seq enrichment: %.1f %%\n",
            mean(chip_enr[usePeaks, "GSE108668_H3"] < 1.0) * 100))
## Percent of binding sites in regions of low H3 ChIP-seq enrichment: 97.1 %
cat(sprintf("Percent of binding sites in regions of high H3K14ac ChIP-seq enrichment: %.1f %%\n",
            mean(chip_enr[usePeaks, "GSE108668_H3K14ac"] >= 1.0) * 100))
## Percent of binding sites in regions of high H3K14ac ChIP-seq enrichment: 11.4 %

Figure 2

Settings

colLabel_150 <- expand.grid(c("Untagged", "Ace2", "Rst2", "Pap1"),
                            c("tube", "plate")) |>
    tidyr::unite(col = "col", Var1, Var2) |>
    dplyr::pull(col)
colLabel_500 <- expand.grid(c("Untagged", "Atf1", "Moc3", "Ntu1",
                              "Pcr1", "Ntu2"),
                            c("tube", "plate")) |>
    tidyr::unite(col = "col", Var1, Var2) |>
    dplyr::pull(col)
rowLabel <- NULL

complexLabel <- c("NuA4", "SAGA")

Summary heatmaps using results from the statistical tests (150 mM NaCl)

Overall heatmap

hm1a <- Heatmap(hmdata_150$tstat, 
                cluster_rows = FALSE, 
                cluster_columns = FALSE, 
                col = makeHeatmapCol(stringency = "high"), 
                border = TRUE, 
                border_gp = gpar(lwd = 0.5),
                column_split = hmdata_150$colsplit, 
                name = "Mod. t-stat.", 
                na_col = binary_heatmap_colors["FALSE"],
                column_names_gp = gpar(fontsize = fs),
                row_names_gp = gpar(fontsize = fs),
                row_title = "Proteins significantly enriched in at least one experiment",
                row_title_gp = gpar(fontsize = fl),
                show_row_names = FALSE, show_column_names = FALSE,
                use_raster = FALSE, show_heatmap_legend = TRUE,
                right_annotation = makeComplexAnnotation(
                    hmdata_150$tstat, complexes[complexLabel], 
                    idmap = idmap, show_legend = TRUE, fontsize = fl, bg_color = bg_color),
                bottom_annotation = makeColLabels(hmdata_150$tstat, colLabel_150, fl),
                width = unit(w_ipmsHm, "mm"), # heatmap body width
                heatmap_height = unit(h_ipmsHm, "mm"), # whole heatmap height
                column_title_gp = gpar(fontsize = ft),
                heatmap_legend_param = list(title_gp = gpar(fontsize = fl), 
                                            legend_direction = "horizontal",
                                            title_position = "topcenter",
                                            border = "gray10",
                                            border_gp = gpar(lwd = 0.5)))

set.seed(42L)
hm_150 <- grid.grabExpr(
    hm1a <- draw(hm1a, merge_legend = TRUE,
                 heatmap_legend_side = "bottom",
                 annotation_legend_side = "bottom"),
    width = w_ipmsHm, height = h_ipmsHm
)

plot_grid(hm_150)

Summary heatmaps using results from the statistical tests (500 mM NaCl)

Overall heatmap

hm2a <- Heatmap(hmdata_500$tstat, 
                cluster_rows = FALSE, 
                cluster_columns = FALSE, 
                col = makeHeatmapCol(stringency = "high"), 
                border = TRUE, 
                border_gp = gpar(lwd = 0.5),
                column_split = hmdata_500$colsplit, 
                name = "Mod. t-stat.", 
                na_col = binary_heatmap_colors["FALSE"],
                column_names_gp = gpar(fontsize = fs),
                row_names_gp = gpar(fontsize = fs),
                row_title = "Proteins significantly enriched in at least one experiment",
                row_title_gp = gpar(fontsize = fl),
                show_row_names = FALSE, show_column_names = FALSE,
                use_raster = FALSE, show_heatmap_legend = TRUE,
                right_annotation = makeComplexAnnotation(
                    hmdata_500$tstat, complexes[complexLabel],
                    idmap = idmap, show_legend = FALSE, fontsize = fl, bg_color = bg_color),
                bottom_annotation = makeColLabels(hmdata_500$tstat, colLabel_500, fl),
                width = unit(w_ipmsHm, "mm"), # heatmap body width
                heatmap_height = unit(h_ipmsHm, "mm"), # whole heatmap height
                column_title_gp = gpar(fontsize = ft),
                heatmap_legend_param = list(title_gp = gpar(fontsize = fl), 
                                            legend_direction = "horizontal",
                                            title_position = "topcenter",
                                            border = "gray10",
                                            border_gp = gpar(lwd = 0.5)))

hm_500 <- grid.grabExpr(
    hm2a <- draw(hm2a, merge_legend = TRUE,
                 heatmap_legend_side = "bottom",
                 annotation_legend_side = "bottom"),
    width = w_ipmsHm, height = h_ipmsHm
)

plot_grid(hm_500)

Supplementary plots - IP-MS data for HOT TFs

sce150untag <- readRDS(params$ipms150untag)

res150untag <- .getTestCols(sce150untag, adjp_cutoff = 0.01, logfc_cutoff = 1)
tstat <- res150untag$tstat
colnames(tstat) <- .getProteinNameFromComparison(colnames(tstat))
tstat <- tstat[rownames(tstat) %in% colnames(tstat), ]
tstat[is.na(tstat)] <- 0
dim(tstat)
## [1] 85 89
hot_tfs <- c("Php3", "Sak1", "Pcr1", "Prr1", "Atf1", "Rst2", "Adn2", 
             "Adn3", "Hsr1", "Phx1", "Pho7")
stopifnot(all(hot_tfs %in% rownames(tstat)))
stopifnot(all(hot_tfs %in% colnames(tstat)))
stopifnot(all(rownames(tstat) %in% colnames(tstat)))
tstat <- tstat[c(hot_tfs, setdiff(rownames(tstat), hot_tfs)), ]
tstat <- tstat[, c(rownames(tstat), setdiff(colnames(tstat), rownames(tstat)))]

hmipmshot <- Heatmap(
    tstat, 
    col = makeHeatmapCol(stringency = "high"), 
    border = TRUE,
    row_split = ifelse(rownames(tstat) %in% hot_tfs, "HOT TFs", "Other TFs"), 
    column_split = ifelse(colnames(tstat) %in% hot_tfs, "HOT TFs", "Other TFs"),
    cluster_rows = FALSE, cluster_columns = FALSE,
    row_title_gp = gpar(fontsize = fl), show_row_names = FALSE, 
    row_names_gp = gpar(fontsize = fs),
    column_title_gp = gpar(fontsize = fl), show_column_names = FALSE,
    heatmap_legend_param = list(labels_gp = gpar(fontsize = fs),
                                border = "gray10",
                                legend_direction = "horizontal",
                                legend_width = unit(24, "mm"),
                                title_position = "lefttop",
                                title_gp = gpar(fontsize = fl)),
    na_col = bg_color,
    name = "Mod. t-stat.")

hm_ipms_hot <- grid.grabExpr(
    hmipmshot <- draw(hmipmshot, merge_legend = TRUE, 
                      heatmap_legend_side = "bottom")
)
plot_grid(hm_ipms_hot)

Peak vs. sample heatmap (binary)

The following creates heatmap of peaks (rows) versus ChIP-seq experiments (columns). The values are binary (TRUE or FALSE) and indicate if a peak was enriched in ChIP-seq experiment (had an average log2 IP/input enrichment values greater than 1.0).

# extract binary enrichment matrix from `peakgr`
m <- as.matrix(mcols(peakgr)[, grep("^is_enr_in[.]", colnames(mcols(peakgr)))])
mode(m) <- "numeric"
colnames(m) <- sub("is_enr_in.", "", colnames(m))

# specify ChIP-seq experiments for which to show labels
selLabels <- c(hot_tfs, "Untagged")
stopifnot(all(selLabels %in% colnames(m)))

# only show peaks that are enriched in at least one ChIP-seq experiment
sel_peaksBin <- rowSums(m) > 0

# prepare annotation data
#   `fractTFs`: fraction of ChIP-seq experiments that a peak is enriched in
#   `numpeaks`: the number of enriched peaks for each ChIP-seq experiment
#   `peaktype`: the type of each peak (specific, frequent or ubiquitous)
#   `near_ncRNA`: specifies if the peak is nearer than `maxdist` from
#                 a %S rRNA or tRNA gene
fracTFs <- rowMeans(m > 0)
numpeaks <- colSums(m > 0)
levels(peakgr$peaktype) <- c("Common (ubiquitous)", # use shorter names
                             "Common (frequent)",
                             "Specific peaks")
peaktype <- factor(gsub("[()]", "", sub(" ", "\n", peakgr$peaktype)),
                   levels = c("Common\nubiquitous",
                              "Common\nfrequent",
                              "Specific\npeaks"))
near_ncRNA <- factor(ifelse(!is.na(dist.peak2tRNA) & dist.peak2tRNA < maxdist,
                            "tRNA",
                            ifelse(!is.na(dist.peak2rRNA) & dist.peak2rRNA < maxdist,
                                   "5S rRNA", NA)),
                         levels = c("tRNA", "5S rRNA"))

# most tRNA and 5S rRNA genes are near "common (ubiquitous)" peaks:
table(peaktype[sel_peaksBin], near_ncRNA[sel_peaksBin])
##                     
##                      tRNA 5S rRNA
##   Common\nubiquitous  107      27
##   Common\nfrequent      1       0
##   Specific\npeaks       5       1
# parameters for annotation legends and colors
annotLegendParams <- list(`%GC` = list(at = c(20, 40, 60),
                                       labels_gp = gpar(fontsize = fs),
                                       border = "gray10",
                                       legend_direction = "horizontal",
                                       legend_width = unit(18, "mm"),
                                       title_position = "topcenter",
                                       title_gp = gpar(fontsize = fl)),
                          `TSS dist.` = list(at = c(0, 4000, 8000),
                                             labels_gp = gpar(fontsize = fs),
                                             border = "gray10",
                                             legend_direction = "horizontal",
                                             legend_width = unit(18, "mm"),
                                             title_position = "topcenter",
                                             title_gp = gpar(fontsize = fl)),
                          `tRNA or 5S rRNA` = list(labels_gp = gpar(fontsize = fs),
                                                   title = ""))
cols_ncRNA <- c("tRNA" = main_colors[5], "5S rRNA" = main_colors[3],
                "none" = bg_color)
annotCols <- list(
    `%GC` = colorRamp2(breaks = seq(min(peakgr$fracGC[sel_peaksBin]),
                                    max(peakgr$fracGC[sel_peaksBin]),
                                    length.out = 64) * 100,
                       colors = rev(hcl.colors(64, "Greens"))),
    ` width` = colorRamp2(breaks = seq(min(width(peakgr)[sel_peaksBin]),
                                       max(width(peakgr)[sel_peaksBin]),
                                       length.out = 64),
                          colors = rev(hcl.colors(64, "YlOrBr"))),
    `TSS dist.` = colorRamp2(breaks = c(0, 10^seq(2, log10(max(peakgr$nrst_TSS_dist)),
                                                  length.out = 63)),
                             colors = hcl.colors(64)),
    `tRNA or 5S rRNA` = cols_ncRNA)

# prepare heatmap annotations
#   `peakAnnot` for peaks (rows), left side
#   `fracTFsAnnot` for peaks (rows), right side
#   `sampleAnnot` for ChIP-seq experiments (columns), top
#   `sampleLabels` for ChIP-seq experiments (columns), bottom
peakAnnot <- HeatmapAnnotation(
                  which = "row",
                  `%GC` = peakgr$fracGC[sel_peaksBin] * 100,
                  `TSS dist.` = peakgr$nrst_TSS_dist[sel_peaksBin],
                  `tRNA or 5S rRNA` = near_ncRNA[sel_peaksBin],
                  col = annotCols, na_col = bg_color,
                  width = unit(21, "mm"), show_legend = TRUE,
                  annotation_name_gp = gpar(fontsize = fl),
                  annotation_legend_param = annotLegendParams)

fracTFsAnnot <- HeatmapAnnotation(
    which = "row",
    `Fraction\nof TFs` = anno_barplot(
        x = fracTFs[sel_peaksBin], gp = gpar(col = na_color),
        border = FALSE, bar_width = 1.0,
        axis_param = list(at = c(0, 0.4, 0.8),
                          gp = gpar(fontsize = fs)),
        width = unit(10, "mm")),
    annotation_name_gp = gpar(fontsize = fl))
sampleAnnot <- HeatmapAnnotation(
    which = "column",
    `Number of\nenriched\npeaks` = anno_barplot(
        x = numpeaks,
        gp = gpar(col = 0, fill = na_color),
        border = FALSE, bar_width = 1.0,
        axis_param = list(at = c(0, 150, 300),
                          gp = gpar(fontsize = fs),
                          facing = "outside",
                          direction = "normal"),
        width = unit(10, "mm")),
    annotation_name_gp = gpar(fontsize = fl),
    annotation_name_side = "left", annotation_name_rot = 0,
    show_legend = TRUE)
sampleLabels <- columnAnnotation(
    TFnames = anno_mark(which = "column", side = "bottom",
                        at = match(selLabels, colnames(m)), 
                        labels = sub("^[^_]+_", "", selLabels),
                        labels_gp = gpar(fontsize = fl)))

# create main heatmap with annotations
hm3 <- Heatmap(m[sel_peaksBin, ],
               column_title = "Binarized ChIP enrichments",
               column_title_side = "top",
               column_title_gp = gpar(fontsize = ft),
               col = circlize::colorRamp2(
                   breaks = seq(0, max(m), length.out = 64),
                   colors = colorRampPalette(binary_heatmap_colors[c("FALSE", "TRUE")])(64)),
               border = TRUE, border_gp = gpar(lwd = 0.5),
               cluster_rows = TRUE, row_dend_width = unit(25, "mm"), show_row_dend = FALSE,
               cluster_columns = TRUE, column_dend_height = unit(10, "mm"),
               row_title_gp = gpar(fontsize = fl),
               row_split = peaktype[sel_peaksBin], cluster_row_slices = FALSE,
               left_annotation = peakAnnot,
               right_annotation = fracTFsAnnot,
               top_annotation = sampleAnnot,
               bottom_annotation = sampleLabels,
               show_row_names = FALSE, show_column_names = FALSE,
               use_raster = TRUE, show_heatmap_legend = FALSE,
               width = unit(w_peaksHm, "mm"), # heatmap body width
               heatmap_height = unit(h_peaksHm, "mm")) # whole heatmap height

# draw the heatmap and grab as a grid graphics `grob`
# (needed for combining figure panels)
# use defined random number seed to make clustering deterministic
set.seed(42L)
hm_peaksBin <- grid.grabExpr(
    hm3 <- draw(hm3, merge_legend = TRUE,
                heatmap_legend_side = "bottom",
                annotation_legend_side = "bottom"),
    width = w_peaksHm, height = h_peaksHm
)

plot_grid(hm_peaksBin)

Peak vs. sample heatmap (enrichments)

The following creates a heatmap similar to the peak-versus-experiment heatmap above, but using the ChIP-seq enrichment values (log2 IP/input) directly instead of binarized (TRUE, FALSE) values.

# use ChIP-seq enrichment matrix from this study (replicates averaged)
m <- peakenrAvg
stopifnot(all(selLabels %in% colnames(m)))

# to keep the same row order as in the binary heatmap `hm3`,
# we have to re-calculate the annotations and switch off the row clustering
sel_peaksEnr <- which(sel_peaksBin)[unlist(row_order(hm3), use.names = FALSE)]

# prepare annotation data
#   `chip_enr_sel`: selected public ChIP-seq experiments
#                   (log2 IP/input enrichments)
#   `chip_enr_sel_range`: value range for color scale (1% to 99%)
#   `chip_enr_sel_range_full`: value range for color scale (0% to 100%)
chip_enr_sel <- chip_enr[sel_peaksEnr, ]
chip_enr_sel_range <- quantile(abs(chip_enr_sel), probs = c(0.95)) * c(-1, 1)
chip_enr_sel_range_full <- range(chip_enr_sel)

# parameters for annotation legends and colors
annotCols2 <- list(
    `tRNA or 5S rRNA` = cols_ncRNA,
    H3 = colorRamp2(breaks = c(chip_enr_sel_range_full[1],
                               seq(chip_enr_sel_range[1],
                                   chip_enr_sel_range[2],
                                   length.out = 62),
                               chip_enr_sel_range_full[2]),
                    colors = viridisLite::magma(64)),
    H3K14ac = colorRamp2(breaks = c(chip_enr_sel_range_full[1],
                                    seq(chip_enr_sel_range[1],
                                        chip_enr_sel_range[2],
                                        length.out = 62),
                                    chip_enr_sel_range_full[2]),
                         colors = viridisLite::magma(64)),
    H3K9me2 = colorRamp2(breaks = c(chip_enr_sel_range_full[1],
                                    seq(chip_enr_sel_range[1],
                                        chip_enr_sel_range[2],
                                        length.out = 62),
                                    chip_enr_sel_range_full[2]),
                         colors = viridisLite::magma(64)),
    H3K9me3 = colorRamp2(breaks = c(chip_enr_sel_range_full[1],
                                    seq(chip_enr_sel_range[1],
                                        chip_enr_sel_range[2],
                                        length.out = 62),
                                    chip_enr_sel_range_full[2]),
                         colors = viridisLite::magma(64)))
annotLegendParams2 <- list(
    `tRNA or 5S rRNA` = list(labels_gp = gpar(fontsize = fs),
                             title = ""),
    H3 = list(at = c(-6, 0, 6),
              title = "Histone log2 IP/input",
              labels_gp = gpar(fontsize = fs),
              border = "gray10",
              legend_direction = "horizontal",
              legend_width = unit(30, "mm"),
              title_position = "topcenter",
              title_gp = gpar(fontsize = fl)))

# prepare heatmap annotations
#   `peakAnnot` for peaks (rows), left side
#   `peakAnnot2` for peaks (rows), right side
peakAnnot <- HeatmapAnnotation(
                  which = "row",
                  `tRNA or 5S rRNA` = near_ncRNA[sel_peaksEnr],
                  H3 = chip_enr_sel[, "GSE108668_H3"],
                  H3K14ac = chip_enr_sel[, "GSE108668_H3K14ac"],
                  H3K9me2 = chip_enr_sel[, "GSE182250_H3K9me2"],
                  H3K9me3 = chip_enr_sel[, "GSE182250_H3K9me3"],
                  col = annotCols2, na_col = bg_color,
                  width = unit(28, "mm"),
                  annotation_name_gp = gpar(fontsize = fl),
                  annotation_legend_param = annotLegendParams2,
                  show_legend = c(TRUE, TRUE, FALSE, FALSE, FALSE))
atac_trunc <- pmin(atac_rpkm_avg[sel_peaksEnr], 2000) # cap ATAC signal at 2000
peakAnnot2 <- HeatmapAnnotation(
    which = "row",
    `Accessibility\n1e3 RPKM\nATAC-seq` =
        anno_barplot(x = atac_trunc / 1e3,
                     gp = gpar(col = na_color),
                     border = FALSE, bar_width = 1.0,
                     baseline = 0,
                     axis_param = list(at = c(0, 1, 2),
                                       gp = gpar(fontsize = fs)),
                     width = unit(10, "mm")),
    annotation_name_gp = gpar(fontsize = fl))

# main heatmap value range for color scale
mx <- max(abs(m[sel_peaksEnr, ]))
qs <- quantile(abs(m[sel_peaksEnr, ]), .99)

# create main heatmap with annotations
hm4 <- Heatmap(m[sel_peaksEnr, ], name = "log2 IP/input",
               column_title = "ChIP enrichments",
               column_title_side = "top",
               column_title_gp = gpar(fontsize = ft),
               col = circlize::colorRamp2(
                   breaks = c(-mx, seq(-qs, qs, length.out = 62), mx),
                   colors = colorRampPalette(enrichment_heatmap_colors)(64)),
               cluster_rows = FALSE,
               cluster_columns = column_dend(hm3), column_dend_height = unit(10, "mm"),
               row_title_gp = gpar(fontsize = fl),
               row_split = peaktype[sel_peaksEnr], cluster_row_slices = FALSE,
               border = TRUE, border_gp = gpar(lwd = 0.5),
               left_annotation = peakAnnot,
               right_annotation = peakAnnot2,
               top_annotation = sampleAnnot,
               bottom_annotation = sampleLabels,
               show_row_names = FALSE, show_column_names = FALSE,
               use_raster = TRUE, show_heatmap_legend = TRUE,
               heatmap_legend_param = list(at = round(c(-mx, 0, mx), 1),
                                           labels_gp = gpar(fontsize = fs),
                                           border = "gray10",
                                           legend_direction = "horizontal",
                                           legend_width = unit(30, "mm"),
                                           title = "TF log2 IP/input",
                                           title_position = "topcenter",
                                           title_gp = gpar(fontsize = fl)),
               width = unit(w_peaksHm, "mm"), # heatmap body width
               heatmap_height = unit(h_peaksHm, "mm"))

# draw the heatmap and grab as a grid graphics `grob`
# (needed for combining figure panels)
# use defined random number seed to make clustering deterministic
set.seed(42L)
hm_peaksEnr <- grid.grabExpr(
    hm4 <- draw(hm4, merge_legend = TRUE,
                heatmap_legend_side = "bottom",
                annotation_legend_side = "bottom"),
    width = w_peaksHm, height = h_peaksHm
)

plot_grid(hm_peaksEnr)

Supplementary plots - ChIP enrichment replicate-pair correlations

pd <- peakenr |> as.data.frame() |>
    tibble::rownames_to_column("peakid") |>
    pivot_longer(col = !matches("peakid")) |>
    mutate(group = sub("_rep[12]$", "", name),
           name = factor(name, levels = unique(name)))

ylims <- c(-2.0, 4.0)
yticks <- c(-2, 0, 2, 4)

pL <- lapply(levels(fct_relevel(sort(unique(pd$group)), "Untagged", after = Inf)),
             function(grp1) {
    pd1 <- pd[pd$group == grp1, ] |> pivot_wider()
    pd1$cols <- densCols(x = pd1[[paste0(grp1, "_rep1")]],
                         y = pd1[[paste0(grp1, "_rep2")]],
                         nbin = 64, colramp = colorRampPalette(hcl.colors(32)))
    ggplot(pd1, aes(.data[[paste0(grp1, "_rep1")]],
                    .data[[paste0(grp1, "_rep2")]])) +
        geom_abline(intercept = 0, slope = 1, linetype = 1, color = "gray") +
        geom_hline(yintercept = 0, linetype = 2, color = "gray") +
        geom_vline(xintercept = 0, linetype = 2, color = "gray") +
        geom_scattermost(xy = as.data.frame(pd1[, paste0(grp1, c("_rep1", "_rep2"))]),
                         pointsize = 2, color = pd1$cols, pixels = c(256, 256)) +
        coord_fixed(xlim = ylims, ylim = ylims, expand = FALSE, clip = "off") +
        scale_x_continuous(breaks = yticks) +
        scale_y_continuous(breaks = yticks) +
        theme_cowplot(7) +
        labs(x = element_blank(), y = element_blank()) +
        annotate("text", x = -Inf, y = Inf, hjust = -0.05, vjust = 1.05,
                 label = grp1, color = "black", size = 3) +
        theme(legend.position = "none")
})

(gg_enrpairs <- plot_grid(plotlist = pL, ncol = 6))

Supplementary plots - Bound promoters

In txannot$prom.number.overlapping.highconf.peaks, we record how many TF peaks are found in the promoter regions of a gene (defined as a 1kb window centered on the annotated transcript start site).

Here we compare the expression levels (obtained from rnaseq_del) of genes twith different numbers of TFs bound at their promoters.

# calculate average expression of each gene
avg_tpm <- rowMeans(assay(rnaseq_del, "abundance"))

# add average expression to gene annotation and filter out
# genes that are not protein coding genes
pd <- txannot |>
    mutate(avg_tpm = avg_tpm[unique_einprot_id]) |>
    filter(!is.na(avg_tpm),
           gene_biotype == "protein_coding") |>
    mutate(prom.number.overlapping.highconf.peaks.categories = 
               cut(prom.number.overlapping.highconf.peaks,
                   breaks = c(-Inf, 0, 1, Inf),
                   labels = c("0", "1", ">1"),
                   right = TRUE))

# visualize
(gg_rnaviolin <- ggplot(pd, aes(prom.number.overlapping.highconf.peaks.categories, avg_tpm + 0.01)) +
        geom_violin(draw_quantiles = 0.5, fill = bg_color) +
        scale_y_log10(labels = scales::label_number(big.mark = "")) +
        labs(x = "Number of promoter TF peaks",
             y = "RNA level\n(TPM + 0.01, log-scale)") +
        theme_cowplot())

Supplementary plots - PCA (IP-MS)

## Low salt
sce150_pca <- einprot::doPCA(
    sce150[rowSums(!assay(sce150, 
                          metadata(sce150)$aNames$assayImputIndic)) >= 2, ],
    metadata(sce150)$aNames$assayImputed)$sce
df <- scuttle::makePerCellDF(sce150_pca, features = NULL, use.coldata = TRUE, 
                             use.dimred = TRUE) |>
    dplyr::mutate(Protocol = .capitalize(sub("_digest|_freeze", "", sub("SOP_", "", Protocol))))
pcvar150 <- attr(reducedDim(sce150_pca, paste0("PCA_", metadata(sce150)$aNames$assayImputed)), 
                 "percentVar")
(ggpca1 <- ggplot(
    df, 
    aes(x = .data[[paste0("PCA_", metadata(sce150)$aNames$assayImputed, ".1")]], 
        y = .data[[paste0("PCA_", metadata(sce150)$aNames$assayImputed, ".2")]], 
        color = Protocol)) + 
    geom_point(size = 3, alpha = 0.5) +
    scale_color_manual(values = c(Plate = main_colors[1], 
                                  Tube = main_colors[5]), name = "") + 
    coord_fixed() + 
    labs(x = paste0("PC1 (", round(pcvar150[1], 1), "%)"), 
         y = paste0("PC2 (", round(pcvar150[2], 1), "%)"), 
         title = "150 mM NaCl IP-MS") + 
    theme_cowplot() + 
    theme(legend.position = "bottom") + 
    guides(color = guide_legend(nrow = 1, byrow = TRUE)))

(ggpca2 <- ggplot(
    df |>
        dplyr::mutate(bait = ifelse(Gene_name == "atf1", "Atf1", "Other baits")) |>
        dplyr::mutate(prot = .capitalize(prot)), 
    aes(x = .data[[paste0("PCA_", metadata(sce150)$aNames$assayImputed, ".1")]], 
        y = .data[[paste0("PCA_", metadata(sce150)$aNames$assayImputed, ".2")]], 
        color = bait, shape = prot)) + 
    geom_point(size = 3, alpha = 0.5) +
    scale_color_manual(values = c(Atf1 = main_colors[3], 
                                  `Other baits` = na_color), name = "") + 
    scale_shape_manual(values = c(Plate = 16, Tube = 17), name = "") + 
    coord_fixed() + 
    labs(x = paste0("PC1 (", round(pcvar150[1], 1), "%)"), 
         y = paste0("PC2 (", round(pcvar150[2], 1), "%)"),
         title = "150 mM NaCl IP-MS") + 
    theme_cowplot() + 
    theme(legend.position = "bottom", legend.box = "vertical", 
          legend.box.just = "left", legend.margin = margin()))

sce500_pca <- einprot::doPCA(
    sce500[rowSums(!assay(sce500, metadata(sce500)$aNames$assayImputIndic)) >= 2, ],
    metadata(sce500)$aNames$assayImputed)$sce
df <- scuttle::makePerCellDF(sce500_pca, features = NULL, use.coldata = TRUE, 
                             use.dimred = TRUE)
pcvar500 <- attr(reducedDim(sce500_pca, paste0("PCA_", metadata(sce150)$aNames$assayImputed)), 
                 "percentVar")
(ggpca3 <- ggplot(
    df |>
        dplyr::mutate(bait = ifelse(Gene_name == "atf1", "Atf1", "Other baits")), 
    aes(x = .data[[paste0("PCA_", metadata(sce500)$aNames$assayImputed, ".1")]], 
        y = .data[[paste0("PCA_", metadata(sce500)$aNames$assayImputed, ".2")]], 
        color = bait)) + 
    geom_point(size = 3, alpha = 0.5) +
    scale_color_manual(values = c(Atf1 = main_colors[3], 
                                  `Other baits` = na_color), name = "") + 
    coord_fixed() + 
    labs(x = paste0("PC1 (", round(pcvar500[1], 1), "%)"), 
         y = paste0("PC2 (", round(pcvar500[2], 1), "%)"), 
         title = "500 mM NaCl IP-MS") + 
    theme_cowplot() + 
    theme(legend.position = "bottom", legend.box.just = "left"))

Put together

cowplot::plot_grid(
    hm_150,
    hm_500,
    hm_peaksEnr,
    hm_peaksBin,
    nrow = 2,
    labels = c("A", "B", "C", "D"),
    align = "v",
    axis = "b"
)

Supplementary figure

Assemble the panels into Supplementary Figure 2:

cowplot::plot_grid(
    cowplot::plot_grid(
        ggpca2,
        ggpca3,
        cowplot::plot_grid(NULL, gg_rnaviolin, ncol = 1,
                           rel_heights = c(0.1, 0.9)),
        cowplot::plot_grid(NULL, hm_ipms_hot, ncol = 1,
                           rel_heights = c(0.1, 0.9)),
        ncol = 1,
        rel_heights = c(1.1, 0.9, 0.8, 1.2),
        labels = c("", "", "C", "D")
    ),
    gg_enrpairs,
    rel_widths = c(0.65, 1), 
    scale = 0.9, nrow = 1,
    labels = c("A", "B"),
    vjust = 5.2, hjust = -2.2
)

Session info

Session info
sessionInfo()
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
## 
## Matrix products: default
## BLAS/LAPACK: FlexiBLAS OPENBLAS;  LAPACK version 3.10.1
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Zurich
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] einprot_0.9.4               tidyr_1.3.1                
##  [3] scattermore_1.2             circlize_0.4.16            
##  [5] colorspace_2.1-0            ComplexHeatmap_2.18.0      
##  [7] jsonlite_1.8.8              forcats_1.0.0              
##  [9] dplyr_1.1.4                 cowplot_1.1.3              
## [11] ggplot2_3.5.0               SingleCellExperiment_1.24.0
## [13] SummarizedExperiment_1.32.0 Biobase_2.62.0             
## [15] GenomicRanges_1.54.1        GenomeInfoDb_1.38.8        
## [17] IRanges_2.36.0              S4Vectors_0.40.2           
## [19] BiocGenerics_0.48.1         MatrixGenerics_1.14.0      
## [21] matrixStats_1.3.0          
## 
## loaded via a namespace (and not attached):
##   [1] R.methodsS3_1.8.2           ggalt_0.4.0                
##   [3] poweRlaw_0.80.0             DT_0.33                    
##   [5] Biostrings_2.70.3           vctrs_0.6.5                
##   [7] digest_0.6.35               png_0.1-8                  
##   [9] shape_1.4.6.1               pcaPP_2.0-4                
##  [11] ggrepel_0.9.5               iSEE_2.14.0                
##  [13] magick_2.8.3                MASS_7.3-60                
##  [15] reshape2_1.4.4              httpuv_1.6.15              
##  [17] foreach_1.5.2               withr_3.0.0                
##  [19] xfun_0.43                   survival_3.5-8             
##  [21] memoise_2.0.1               hexbin_1.28.3              
##  [23] ggbeeswarm_0.7.2            gmm_1.8                    
##  [25] systemfonts_1.0.6           zoo_1.8-12                 
##  [27] GlobalOptions_0.1.2         gtools_3.9.5               
##  [29] R.oo_1.26.0                 DEoptimR_1.1-3             
##  [31] GGally_2.2.1                KEGGREST_1.42.0            
##  [33] promises_1.3.0              httr_1.4.7                 
##  [35] restfulr_0.0.15             hash_2.2.6.3               
##  [37] rstudioapi_0.16.0           shinyAce_0.4.2             
##  [39] rrcovNA_0.5-1               miniUI_0.1.1.1             
##  [41] generics_0.1.3              babelgene_22.9             
##  [43] zlibbioc_1.48.2             ScaledMatrix_1.10.0        
##  [45] GenomeInfoDbData_1.2.11     SparseArray_1.2.4          
##  [47] xtable_1.8-4                stringr_1.5.1              
##  [49] ade4_1.7-22                 pracma_2.4.4               
##  [51] doParallel_1.0.17           evaluate_0.23              
##  [53] ExploreModelMatrix_1.14.0   S4Arrays_1.2.1             
##  [55] proDA_1.16.0                hms_1.1.3                  
##  [57] irlba_2.3.5.1               iSEEhex_1.4.0              
##  [59] shinyWidgets_0.8.5          magrittr_2.0.3             
##  [61] readr_2.1.5                 later_1.3.2                
##  [63] viridis_0.6.5               lattice_0.22-6             
##  [65] MsCoreUtils_1.14.1          genefilter_1.84.0          
##  [67] robustbase_0.99-2           XML_3.99-0.16.1            
##  [69] scuttle_1.12.0              pillar_1.9.0               
##  [71] nlme_3.1-164                iterators_1.0.14           
##  [73] caTools_1.18.2              compiler_4.3.2             
##  [75] beachmat_2.18.1             stringi_1.8.3              
##  [77] GenomicAlignments_1.38.2    tmvtnorm_1.6               
##  [79] plyr_1.8.9                  msigdbr_7.5.1              
##  [81] crayon_1.5.2                abind_1.4-5                
##  [83] BiocIO_1.12.0               scater_1.30.1              
##  [85] chron_2.3-61                bit_4.0.5                  
##  [87] sandwich_3.1-0              pcaMethods_1.94.0          
##  [89] codetools_0.2-19            BiocSingular_1.18.0        
##  [91] bslib_0.7.0                 GetoptLong_1.0.5           
##  [93] plotly_4.10.4               mime_0.12                  
##  [95] MultiAssayExperiment_1.28.0 splines_4.3.2              
##  [97] Rcpp_1.0.12                 sparseMatrixStats_1.14.0   
##  [99] Rttf2pt1_1.3.12             knitr_1.46                 
## [101] blob_1.2.4                  utf8_1.2.4                 
## [103] clue_0.3-65                 seqLogo_1.68.0             
## [105] AnnotationFilter_1.26.0     QFeatures_1.12.0           
## [107] DelayedMatrixStats_1.24.0   tibble_3.2.1               
## [109] sqldf_0.4-11                iSEEu_1.14.0               
## [111] Matrix_1.6-5                statmod_1.5.0              
## [113] tzdb_0.4.0                  svglite_2.1.3              
## [115] pkgconfig_2.0.3             tools_4.3.2                
## [117] cachem_1.0.8                RSQLite_2.3.6              
## [119] viridisLite_0.4.2           DBI_1.2.2                  
## [121] impute_1.76.0               fastmap_1.1.1              
## [123] rmarkdown_2.26              scales_1.3.0               
## [125] shinydashboard_0.7.2        Rsamtools_2.18.0           
## [127] sass_0.4.9                  patchwork_1.2.0            
## [129] ggstats_0.6.0               farver_2.1.1               
## [131] mgcv_1.9-1                  gsubfn_0.7                 
## [133] yaml_2.3.8                  rtracklayer_1.62.0         
## [135] cli_3.6.2                   purrr_1.0.2                
## [137] writexl_1.5.0               lifecycle_1.0.4            
## [139] mvtnorm_1.2-4               rintrojs_0.3.4             
## [141] BiocParallel_1.36.0         annotate_1.80.0            
## [143] gtable_0.3.5                rjson_0.2.21               
## [145] parallel_4.3.2              limma_3.58.1               
## [147] colourpicker_1.3.0          TFBSTools_1.40.0           
## [149] bitops_1.0-7                kableExtra_1.4.0           
## [151] bit64_4.0.5                 BiocNeighbors_1.20.2       
## [153] proto_1.0.0                 CNEr_1.38.0                
## [155] jquerylib_0.1.4             highr_0.10                 
## [157] shinyjs_2.1.0               imputeLCMD_2.1             
## [159] R.utils_2.12.3              rrcov_1.7-5                
## [161] lazyeval_0.2.2              shiny_1.8.1.1              
## [163] htmltools_0.5.8.1           proj4_1.0-14               
## [165] GO.db_3.18.0                glue_1.7.0                 
## [167] STRINGdb_2.14.3             TFMPvalue_0.0.9            
## [169] XVector_0.42.0              RCurl_1.98-1.14            
## [171] ComplexUpset_1.3.3          mclust_6.1                 
## [173] BSgenome_1.70.2             motifStack_1.46.0          
## [175] gridExtra_2.3               igraph_2.0.3               
## [177] extrafontdb_1.0             R6_2.5.1                   
## [179] ggiraph_0.8.9               gplots_3.1.3.1             
## [181] labeling_0.4.3              cluster_2.1.4              
## [183] stringdist_0.9.12           DirichletMultinomial_1.44.0
## [185] DelayedArray_0.28.0         tidyselect_1.2.1           
## [187] vipor_0.4.7                 plotrix_3.8-4              
## [189] ProtGenerics_1.34.0         maps_3.4.2                 
## [191] xml2_1.3.6                  ash_1.0-15                 
## [193] AnnotationDbi_1.64.1        rsvd_1.0.5                 
## [195] munsell_0.5.1               KernSmooth_2.23-22         
## [197] data.table_1.15.4           htmlwidgets_1.6.4          
## [199] RColorBrewer_1.1-3          rlang_1.1.3                
## [201] extrafont_0.19              uuid_1.2-0                 
## [203] norm_1.0-11.1               fansi_1.0.6                
## [205] Cairo_1.6-2                 beeswarm_0.4.0