dinoR-vignette
Introduction
The objective of dinoR is to detect potential
differences in nucleosome binding/accessibility patterns between two
groups of NOMe-seq experiments. We recently developed a method, called
guided NOMe-seq, which allows for obtaining high coverage NOMe-seq data
from more than 1000 regions of interest (ROIs) in the genome at low
cost, thereby allowing us to profile multiple samples at once
(manuscript submitted). For example, as shown below, we generated
NOMe-seq data targeting ROIs spanning CTCF or REST motifs from wild-type
(WT) and ADNP knock-out mouse ES cells, in duplicates. To facilitate
processing of these samples, with a focus on detecting differences in
nucleosome positioning between WT and ADNP KO cells, we developed
dinoR. dinoR focuses on comparing ROIs which
are grouped together based on similar features (for example, which
transcription factor binds the motif in the ROI), and statistical
analysis of NOMe-seq pattern differences between samples. For quality
control and plotting methylation across a single ROI, please refer to
the Bioconductor package SingleMoleculeFootprinting.
Installation
dinoR is available on Bioconductor and can be installed via:
if (!require("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("dinoR")then we can load dinoR and other necessary packages:
suppressPackageStartupMessages({
library(dinoR)
library(ggplot2)
library(dplyr)
})
#> Warning: multiple methods tables found for 'setequal'
#> Warning: replacing previous import 'BiocGenerics::setequal' by
#> 'S4Vectors::setequal' when loading 'IRanges'
#> Warning: replacing previous import 'BiocGenerics::setequal' by
#> 'S4Vectors::setequal' when loading 'GenomeInfoDb'
#> Warning: replacing previous import 'BiocGenerics::setequal' by
#> 'S4Vectors::setequal' when loading 'GenomicRanges'
#> Warning: replacing previous import 'BiocGenerics::setequal' by
#> 'S4Vectors::setequal' when loading 'XVector'
#> Warning: replacing previous import 'BiocGenerics::setequal' by
#> 'S4Vectors::setequal' when loading 'SummarizedExperiment'
#> Warning: replacing previous import 'BiocGenerics::setequal' by
#> 'S4Vectors::setequal' when loading 'S4Arrays'
#> Warning: replacing previous import 'BiocGenerics::setequal' by
#> 'S4Vectors::setequal' when loading 'DelayedArray'
#> Warning: replacing previous import 'BiocGenerics::setequal' by
#> 'S4Vectors::setequal' when loading 'SparseArray'
#> Warning: replacing previous import 'S4Arrays::read_block' by
#> 'DelayedArray::read_block' when loading 'SummarizedExperiment'Load the NOMe-seq data for ADNP Knock-Out and
WT mouse ES cells (two replicates each)
To obtain the example data shown below, we used biscuit to map 300bp paired-end reads to the genome, UMI-tools to remove duplicated UMIs, and the fetch-NOMe package to get the protection from GCH methylation calls for each read pair (fragment) overlapping a region of interest (ROI). We then use the R package NOMeConverteR to convert the resulting tibble into the RangedSummarizedExperiment object that is loaded below. This represents an efficient way of sharing NOMe-seq data.
Alternatively, mapping and obtaining the GCH methylation calls can be
done using QuasR,
or any other bisulfite aware mapping/quantification tool. The results
then need to be converted into a RangedSummarizedExperiment. For an
example NOMe-seq RangedSummarizedExperiment see the example data loaded
below, or the data generated by using
dinoR::createExampleData.
Importantly, the ROIs should all be the same length and centered around a transcription factor motif, with the strand of the motif indicated. That will ensure that the genomic positions around the motif are sorted according to motif strand, which will allow the user to observe potential asymmetries in protection from methylation relative to the TF motif.
data(NomeData)
NomeData
#> class: RangedSummarizedExperiment
#> dim: 219 4
#> metadata(0):
#> assays(5): nFragsFetched nFragsNonUnique nFragsBisFailed nFragsAnalyzed
#> reads
#> rownames(219): Adnp_chr8_47978653_47979275
#> Adnp_chr6_119394879_119395501 ... Rest_chr4_140283342_140283964
#> Rest_chr7_64704080_64704702
#> rowData names(1): motif
#> colnames(4): AdnpKO_1 AdnpKO_2 WT_1 WT_2
#> colData names(2): samples groupThe reads assay contains GPos objects with the GCH methylation data in two sparse logical matrices, one for protection from methylation, and one for methylation.
assays(NomeData)[["reads"]][1, 1]
#> [[1]]
#> UnstitchedGPos object with 623 positions and 2 metadata columns:
#> seqnames pos strand | protection methylation
#> <Rle> <integer> <Rle> | <lgCMatrix> <lgCMatrix>
#> [1] chr8 47978653 + | FALSE:FALSE:FALSE:... FALSE:FALSE:FALSE:...
#> [2] chr8 47978654 + | FALSE:FALSE:FALSE:... FALSE:FALSE:FALSE:...
#> [3] chr8 47978655 + | FALSE:FALSE:FALSE:... FALSE:FALSE:FALSE:...
#> [4] chr8 47978656 + | FALSE:FALSE:FALSE:... FALSE:FALSE:FALSE:...
#> [5] chr8 47978657 + | FALSE:FALSE:FALSE:... FALSE:FALSE:FALSE:...
#> ... ... ... ... . ... ...
#> [619] chr8 47979271 + | FALSE:FALSE:FALSE:... FALSE:FALSE:FALSE:...
#> [620] chr8 47979272 + | FALSE:FALSE:FALSE:... FALSE:FALSE:FALSE:...
#> [621] chr8 47979273 + | FALSE:FALSE:FALSE:... FALSE:FALSE:FALSE:...
#> [622] chr8 47979274 + | FALSE:FALSE:FALSE:... FALSE:FALSE:FALSE:...
#> [623] chr8 47979275 + | FALSE:FALSE:FALSE:... FALSE:FALSE:FALSE:...
#> -------
#> seqinfo: 53 sequences from an unspecified genome; no seqlengthsMeta plots across ROIs with common TF motifs in the center
We generate metaplots, grouping our ROIs into those that have REST, CTCF, or ADNP bound to the motifs in their center. We use 2 samples from WT mouse ES cells, and two samples from ADNP KO mouse ES cells. We exclude any ROI - sample combinations which contain less than 10 reads (nr=10).
avePlotData <- metaPlots(NomeData = NomeData, nr = 10, ROIgroup = "motif")
# plot average plots
ggplot(avePlotData, aes(x = position, y = protection)) +
geom_point(alpha = 0.5) +
geom_line(aes(x = position, y = loess), col = "darkblue", lwd = 2) +
theme_classic() +
facet_grid(rows = vars(type), cols = vars(sample), scales = "free") +
ylim(c(0, 100)) +
geom_hline(
yintercept = c(10, 20, 30, 40, 50, 60, 70, 80, 90),
alpha = 0.5, color = "grey", linetype = "dashed"
)
We can already see that while the NOMe footprints around REST and CTCF bound motifs don’t change, there are clear differences between WT and ADNP KO cells around the ADNP bound motifs.
Determine fragment counts for five chromatin patterns:
TF, open, upNuc, downNuc, Nuc
To quantify the differences visible in above meta plots, we adopted and slightly modified the approach of Sönmezer et al., 2021. We classify each fragment according to five types of footprints: transcription factor bound (TF), open chromatin, and nucleosome (we distinguish also upstream positioned nucleosome (upNuc), downstream positioned nucleosome (downNuc), and all other nucleosome (Nuc) footprints). To do this we use three windows (-50:-25, -8:8, 25:50) around the motif center (which should correspond to the ROI center of the provided ROIs).
NomeData <- footprintCalc(NomeData)
NomeData
#> class: RangedSummarizedExperiment
#> dim: 219 4
#> metadata(0):
#> assays(6): nFragsFetched nFragsNonUnique ... reads footprints
#> rownames(219): Adnp_chr8_47978653_47979275
#> Adnp_chr6_119394879_119395501 ... Rest_chr4_140283342_140283964
#> Rest_chr7_64704080_64704702
#> rowData names(1): motif
#> colnames(4): AdnpKO_1 AdnpKO_2 WT_1 WT_2
#> colData names(2): samples groupThis results in the addition of an assay called footprints, which assigns each fragment one of the five footprint classes described above.
assays(NomeData)[["footprints"]][[1, 1]][1:10]
#> f1:M04764:177:000000000-JPNT5:1:2119:13500:10508_TGAGCTAG_CATCGACTC
#> "tf"
#> f1:M04764:177:000000000-JPNT5:1:2115:9271:12803_TGAGCTAG_AGGCACCTG
#> "tf"
#> f1:M04764:177:000000000-JPNT5:1:2114:19141:6070_TGAGCTAA_ATCTACCCA
#> "tf"
#> f1:M04764:177:000000000-JPNT5:1:2106:11802:15963_TGAGCTAG_CCACTGTTA
#> "tf"
#> f1:M04764:177:000000000-JPNT5:1:2106:9228:15151_TGAGCTAG_TCAATCAAG
#> "tf"
#> f1:M04764:177:000000000-JPNT5:1:2104:15018:10388_TGAGCTAG_CCAGACGCC
#> "tf"
#> f1:M04764:177:000000000-JPNT5:1:2104:13159:6640_TGAGCTAG_TGATGGACA
#> "downNuc"
#> f1:M04764:177:000000000-JPNT5:1:2113:21587:8114_TGAGCTAG_TCTATCAGC
#> "tf"
#> f1:M04764:177:000000000-JPNT5:1:2110:18921:3356_TGAGCTAG_TTAAAGTCT
#> "tf"
#> f1:M04764:177:000000000-JPNT5:1:2104:27164:6456_TGAGCTAG_TATTTCTAC
#> "tf"Then we count the number of fragments in each sample-ROI combination supporting each footprint class. This results in the addition of an assay for each footprint category, and one for the combined counts across all footprint categories (“all”).
NomeData <- footprintQuant(NomeData)
assays(NomeData)[7:12]
#> List of length 6
#> names(6): tf open upNuc Nuc downNuc allNote that if a fragment does not have methylation protection data in all three windows needed for classification, the fragment will not be used.
Next we can test for differential abundance of footprints between ADNP KO and WT samples.
Calculate differential NOMe-seq footprint abundance between ADNP KO and WT
We use edgeR to check for differences in abundance between wild type and ADNP KO samples for each footprint type fragment count compared to the total fragment counts. Library sizes for TMM normalization are calculated on the total fragment counts.
res <- diNOMeTest(NomeData,
WTsamples = c("WT_1", "WT_2"),
KOsamples = c("AdnpKO_1", "AdnpKO_2")
)
res
#> # A tibble: 1,040 × 10
#> logFC logCPM F PValue FDR contrasts ROI motif logadjPval
#> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <chr> <dbl>
#> 1 1.69 10.9 17.3 0.0000332 0.00691 open_vs_all Adnp_chr4… Adnp 2.16
#> 2 1.53 10.4 8.84 0.00299 0.225 open_vs_all Adnp_chr6… Adnp 0.647
#> 3 1.22 9.34 8.50 0.00358 0.225 open_vs_all Adnp_chr8… Adnp 0.647
#> 4 -1.46 10.5 7.83 0.00518 0.225 open_vs_all Ctcf_chr2… Adnp 0.647
#> 5 1.05 10.6 7.75 0.00542 0.225 open_vs_all Adnp_chr1… Adnp 0.647
#> 6 1.13 10.2 6.46 0.0111 0.385 open_vs_all Adnp_chr1… Adnp 0.414
#> 7 1.00 10.8 4.88 0.0273 0.651 open_vs_all Adnp_chr4… Adnp 0.186
#> 8 -1.12 10.9 4.78 0.0288 0.651 open_vs_all Ctcf_chr7… Adnp 0.186
#> 9 0.566 9.86 4.61 0.0318 0.651 open_vs_all Adnp_chr2… Adnp 0.186
#> 10 -1.00 9.62 4.57 0.0326 0.651 open_vs_all Rest_chr7… Adnp 0.186
#> # ℹ 1,030 more rows
#> # ℹ 1 more variable: regulated <chr>We can then simply plot the number of regulated ROIs within each ROI type…
res %>%
group_by(contrasts, motif, regulated) %>%
summarize(n = n()) %>%
ggplot(aes(x = motif, y = n, fill = regulated)) +
geom_bar(stat = "identity") +
facet_grid(~contrasts) +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
scale_fill_manual(values = c("orange", "grey", "blue3"))
#> `summarise()` has grouped output by 'contrasts', 'motif'. You can override
#> using the `.groups` argument.
…or display the results in MA plots.
ggplot(res, aes(y = logFC, x = logCPM, col = regulated)) +
geom_point() +
facet_grid(~contrasts) +
theme_bw() +
scale_color_manual(values = c("orange", "grey", "blue3"))
Calculate the percentage of fragments in each footprint type and
plot a (clustered) heatmap comparing percentages in WT and ADNP KO
Compare the footprint percentages and significance testing
results for ADNP KO and WT
compareFootprints(footprint_percentages, res,
WTsamples = c("WT_1", "WT_2"),
KOsamples = c("AdnpKO_1", "AdnpKO_2"), plotcols =
c("#f03b20", "#a8ddb5", "#bdbdbd")
)
We can see that in ADNP KO samples, transcription factor footprints significantly increase around ADNP motifs, while nucleosome footprints decrease.
Combining the nucleosome patterns
In case we are not interested in the upstream and downstream nucleosome patterns, but would rather keep all nucleosome pattern fragments within the nucleosome group, we can do that using the option combineNucCounts=TRUE.
res <- diNOMeTest(NomeData,
WTsamples = c("WT_1", "WT_2"),
KOsamples = c("AdnpKO_1", "AdnpKO_2"), combineNucCounts = TRUE
)
footprint_percentages <- footprintPerc(NomeData, combineNucCounts = TRUE)
compareFootprints(footprint_percentages, res,
WTsamples = c("WT_1", "WT_2"),
KOsamples = c("AdnpKO_1", "AdnpKO_2"), plotcols =
c("#f03b20", "#a8ddb5", "#bdbdbd")
)
Session information
sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
#> [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: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] dplyr_1.1.4 ggplot2_3.5.1
#> [3] dinoR_1.3.0 SummarizedExperiment_1.37.0
#> [5] Biobase_2.67.0 GenomicRanges_1.59.0
#> [7] GenomeInfoDb_1.43.0 IRanges_2.41.0
#> [9] S4Vectors_0.45.0 BiocGenerics_0.53.1
#> [11] generics_0.1.3 MatrixGenerics_1.19.0
#> [13] matrixStats_1.4.1 rmarkdown_2.29
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 circlize_0.4.16 shape_1.4.6.1
#> [4] rjson_0.2.23 xfun_0.49 bslib_0.8.0
#> [7] GlobalOptions_0.1.2 lattice_0.22-6 vctrs_0.6.5
#> [10] tools_4.4.2 parallel_4.4.2 fansi_1.0.6
#> [13] tibble_3.2.1 highr_0.11 cluster_2.1.6
#> [16] pkgconfig_2.0.3 Matrix_1.7-1 RColorBrewer_1.1-3
#> [19] lifecycle_1.0.4 GenomeInfoDbData_1.2.13 farver_2.1.2
#> [22] stringr_1.5.1 compiler_4.4.2 statmod_1.5.0
#> [25] munsell_0.5.1 codetools_0.2-20 ComplexHeatmap_2.23.0
#> [28] clue_0.3-65 htmltools_0.5.8.1 sys_3.4.3
#> [31] buildtools_1.0.0 sass_0.4.9 yaml_2.3.10
#> [34] tidyr_1.3.1 pillar_1.9.0 crayon_1.5.3
#> [37] jquerylib_0.1.4 limma_3.63.1 DelayedArray_0.33.1
#> [40] cachem_1.1.0 iterators_1.0.14 abind_1.4-8
#> [43] foreach_1.5.2 locfit_1.5-9.10 tidyselect_1.2.1
#> [46] digest_0.6.37 stringi_1.8.4 purrr_1.0.2
#> [49] labeling_0.4.3 maketools_1.3.1 cowplot_1.1.3
#> [52] fastmap_1.2.0 grid_4.4.2 colorspace_2.1-1
#> [55] cli_3.6.3 SparseArray_1.7.1 magrittr_2.0.3
#> [58] S4Arrays_1.7.1 utf8_1.2.4 withr_3.0.2
#> [61] edgeR_4.5.0 UCSC.utils_1.3.0 scales_1.3.0
#> [64] XVector_0.47.0 httr_1.4.7 png_0.1-8
#> [67] GetoptLong_1.0.5 evaluate_1.0.1 knitr_1.48
#> [70] doParallel_1.0.17 rlang_1.1.4 glue_1.8.0
#> [73] jsonlite_1.8.9 R6_2.5.1 zlibbioc_1.52.0- Introduction
- Installation
- Load the NOMe-seq data for ADNP Knock-Out and
- WT mouse ES cells (two replicates each)
- Meta plots across ROIs with common TF motifs in the center
- Determine fragment counts for five chromatin patterns:
- TF, open, upNuc, downNuc, Nuc
- Calculate differential NOMe-seq footprint abundance between ADNP KO and WT
- Calculate the percentage of fragments in each footprint type and
- plot a (clustered) heatmap comparing percentages in WT and ADNP KO
- Compare the footprint percentages and significance testing
- results for ADNP KO and WT
- Combining the nucleosome patterns
- Session information