library(mesa)
#> Loading required package: qsea
#> Loading required package: dplyr
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
#>
#> Attaching package: 'mesa'
#> The following object is masked from 'package:qsea':
#>
#> getPCA
library(stringr)
library(tibble)In this vignette we cover the generation of a qseaSet using the
functions provided in mesa, rather than those in
qsea. Please see the qsea package vignette
vignette("qsea_tutorial", package = "qsea") for more
detail, although we cover much of it here.
The first thing to prepare is the sampleTable: a data
frame containing metadata for the samples to be analysed. This has three
required columns:
sample_name : a name for each sample, which must be
unique. This is the primary identifier used to describe each sample, so
it is good to be descriptive.group : a column to specify a group for each sample,
for use inside certain functions for averaging over samples of the same
type. We suggest only using this for replicates of the same (or very
similar) samples, and adding extra columns for other comparisons that
may be of interest.file_name : a path to a bam file, containing the
aligned, deduplicated reads mapped to the relevant reference genome,
from material which has undergone the methylation enrichment/capture
process (e.g. MeDIP-seq or MBD-seq). Ideally
samtools fixmate has been run to supply mate information
onto the read tags, but filtering for quality is not required.There are also two optional fields that may be used during generation:
input_file : a path to an additional bam file
containing aligned, deduplicated reads for material from the same
sample, without performing enrichment. This is whole-genome sequencing,
but may be sequenced at a much lower depth (e.g. ultra-shallow WGS, ~1-2
million fragments is sufficient for calling copy number variation at a
window size of 1Mb). If present, this will be used on a sample-by-sample
basis to calculate copy number variation, at a chromosome scale, in
order to correct for this in the expected read depths. This is most
relevant when studying cancer, as particular cancers often show
recurrent patterns of copy number changes. Note that below ~0.5 million
fragments the calculated copy number patterns will become especially
noisy (for human samples).
gender or sex : if one of these columns
is included, then any chromosome called X, Y,
chrX or chrY will be set with the appropriate
zygosity for samples with M, m or
male in this column. See qsea.createSet.R in
the qsea source for details; within the National Biomarker
Centre we only consider the autosomes.
In addition, any further metadata columns may be included as required, and we recommend that anything you might want to query or compare the data on is included here. This includes experimental conditions (tissue, treatment, tumour status, etc.).
We also highly recommend the use of a tool such as NGSCheckMate for verifying that (human) samples come from the same individual or not. This is useful both for ensuring that the enriched and non-enriched samples do indeed come from the same individual, and to check that paired samples are actually paired (for instance when comparing pre- and post-treatment).
For the example here, we use the example data in the
MEDIPSData package. This contains pairs of tumour and
normal tissue from the same patients, with MeDIP-seq performed as the
enrichment step. These bam files are downsampled to only contain
chromosomes 20-22, and are aligned to hg19. Note, we highly recommend
using the hg38/GRCh38 reference genome rather than hg19 as used
here!
sampleTable <- tribble(~sample_name, ~group, ~tumour, ~patient, ~file_name,
"Normal1", "Normal1", "Normal", "Patient1", system.file(
"extdata", "NSCLC_MeDIP_1N_fst_chr_20_21_22.bam",
package = "MEDIPSData", mustWork = TRUE
),
"Normal2", "Normal1", "Normal", "Patient2", system.file(
"extdata", "NSCLC_MeDIP_2N_fst_chr_20_21_22.bam",
package = "MEDIPSData", mustWork = TRUE
),
"Tumour1", "Tumour1", "Tumour", "Patient1", system.file(
"extdata", "NSCLC_MeDIP_1T_fst_chr_20_21_22.bam",
package = "MEDIPSData", mustWork = TRUE
),
"Tumour2", "Tumour2", "Tumour", "Patient2", system.file(
"extdata", "NSCLC_MeDIP_2T_fst_chr_20_21_22.bam",
package = "MEDIPSData", mustWork = TRUE
)
)The qsea tutorial details the steps involved in
generating a qseaSet. We provide a function,
makeQset, which wraps these up into one function, although
you may still use the rest of the mesa functionality with
or without using this to create your qseaSet.
Before creating the qseaSet, you need to have:
sampleTable of metadata, including the file paths
to the bam files.BSgenome package, which must be
installed. You can check which BSgenome packages are available on
Bioconductor via BSgenome::available.genomes(), or make
your own if necessary
(vignette("BSgenomeForge", package = "BSgenome")).qseaSet, the
genome is split into fixed windows with equal size, and all analysis is
performed on these windows as independent units. We generally use 300bp
windows, as that is a little bigger than our average DNA fragment sizes,
which are in the 150-250bp range. If the DNA fragments are particularly
long, a larger window size may make sense. This should also be informed
by your depth of sequencing, as a smaller window size will be more prone
to noise.For our hg19 example data, we set
BSgenome = "BSgenome.Hsapiens.UCSC.hg19", specify the three
chromosomes included in the test dataset and give some approximate
fragment length and SD values.
makeQset(sampleTable,
BSgenome = "BSgenome.Hsapiens.UCSC.hg19",
chrSelect = paste0("chr",20:22),
windowSize = 300,
fragmentLength = 200,
fragmentSD = 50,
CNVmethod = "None")
#>
#> Attaching package: 'generics'
#> The following object is masked from 'package:dplyr':
#>
#> explain
#> The following objects are masked from 'package:base':
#>
#> as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
#> setequal, union
#>
#> Attaching package: 'BiocGenerics'
#> The following object is masked from 'package:dplyr':
#>
#> combine
#> The following objects are masked from 'package:stats':
#>
#> IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#>
#> anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#> colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#> get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
#> order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
#> rbind, Reduce, rownames, sapply, saveRDS, table, tapply, unique,
#> unsplit, which.max, which.min
#>
#> Attaching package: 'S4Vectors'
#> The following objects are masked from 'package:dplyr':
#>
#> first, rename
#> The following object is masked from 'package:utils':
#>
#> findMatches
#> The following objects are masked from 'package:base':
#>
#> expand.grid, I, unname
#>
#> Attaching package: 'IRanges'
#> The following objects are masked from 'package:dplyr':
#>
#> collapse, desc, slice
#>
#> Attaching package: 'Biostrings'
#> The following object is masked from 'package:base':
#>
#> strsplit
#> Considering 541532 regions with total size 162459600
#> ==== Creating qsea set ====
#> restricting analysis on chr20, chr21, chr22
#> Dividing provided regions in 300nt windows
#> Scanning one file at a time. Use e.g. setMesaParallel(nCores = 4) to parallelise this process.
#> Reading NSCLC_MeDIP_1N_fst_chr_20_21_22.bam, with minimum mapping quality 10.
#> No samtools fixmate MQ (mate quality) tags on the bam file, using R1 MAPQ only.
#> Dropping 20294 proper pairs that have neither read with mapq >= 10
#> Dropping 179 proper pairs that do not lie within [70, 1000]
#> Using high quality R1 unpaired reads as well as pairs
#> Estimating fragment length for 0 unpaired R1 reads as 186, based on 879638 proper pairs.
#> Reading NSCLC_MeDIP_2N_fst_chr_20_21_22.bam, with minimum mapping quality 10.
#> No samtools fixmate MQ (mate quality) tags on the bam file, using R1 MAPQ only.
#> Dropping 22848 proper pairs that have neither read with mapq >= 10
#> Dropping 312 proper pairs that do not lie within [70, 1000]
#> Using high quality R1 unpaired reads as well as pairs
#> Estimating fragment length for 0 unpaired R1 reads as 188, based on 1132561 proper pairs.
#> Reading NSCLC_MeDIP_1T_fst_chr_20_21_22.bam, with minimum mapping quality 10.
#> No samtools fixmate MQ (mate quality) tags on the bam file, using R1 MAPQ only.
#> Dropping 12436 proper pairs that have neither read with mapq >= 10
#> Dropping 146 proper pairs that do not lie within [70, 1000]
#> Using high quality R1 unpaired reads as well as pairs
#> Estimating fragment length for 0 unpaired R1 reads as 183, based on 637854 proper pairs.
#> Reading NSCLC_MeDIP_2T_fst_chr_20_21_22.bam, with minimum mapping quality 10.
#> No samtools fixmate MQ (mate quality) tags on the bam file, using R1 MAPQ only.
#> Dropping 11131 proper pairs that have neither read with mapq >= 10
#> Dropping 168 proper pairs that do not lie within [70, 1000]
#> Using high quality R1 unpaired reads as well as pairs
#> Estimating fragment length for 0 unpaired R1 reads as 189, based on 632731 proper pairs.
#> Generated objects
#> Calculated coverage
#> No CNV being calculated
#> Get genomic positions of "CG" ...
#> searching chr20 for "CG"...
#> found 717722 occurances of CG in chr20
#> estimating expected number of CGs per fragment for windows of chr20...
#> ...done
#> searching chr21 for "CG"...
#> found 380444 occurances of CG in chr21
#> estimating expected number of CGs per fragment for windows of chr21...
#> ...done
#> searching chr22 for "CG"...
#> found 578097 occurances of CG in chr22
#> estimating expected number of CGs per fragment for windows of chr22...
#> ...done
#> selecting windows with low CpG density for background read estimation
#> 4.396% of the windows have enrichment pattern density of at most 0.05 per fragment
#> and are used for background reads estimation
#> qseaSet object generated successfully
#> qseaSet
#> =======================================
#> 4 Samples:
#> Normal1, Normal2, Tumour1, Tumour2
#> 541532 Regions in 3 chromosomes of BSgenome.Hsapiens.UCSC.hg19A range of messages are produced as output during this process. Note
that by default the process is performed without parallelisation; this
can be controlled via the setMesaParallel function and will
vastly speed up the process when multiple cores are available (but use
more RAM).
This function combines several steps:
BSgenome, chrSelect, windowSize
and badRegions parameters.fragmentLength and fragmentSD parameters to
define the expected contribution of fragments to the likelihood that
they overlap the window (see the Calculating CpG Density section
below).The chosen chromosomes (via chrSelect) of the
BSgenome object will be tiled in windowSize
windows. These can be filtered for regions which have poor mappability
or are overrepresented in non-enriched sequencing, by providing a
GRanges object to the badRegions option. If provided, this
will remove any window which overlaps with the GRanges object; ensuring
that the windows remain fixed when this filtering is changed. The 2019
ENCODE exclusion regions are provided as a data object called
ENCODEbadRegions.
By default, the coverage is calculated using a custom method called
coverageMethod = "PairedAndR1s". This assumes paired-end
sequencing, and will first take all the (high quality) paired-end reads
within the bam file, then it will take (high quality) first reads (R1s)
whose mate (R2) did not map and extend them to the average fragment size
of the paired-end reads. Here, high quality means the MAPQ value is
above minMapQual, although for paired ends we accept pairs
where at least one end has a MAPQ above this value.
There are a number of options for controlling this process:
properPairsOnly - when set to TRUE, only the proper
pairs will be included, not the R1s without a mapped mate.minInsertSize and maxInsertSize control
the size distribution of DNA fragments to be included. This enables
filtering out of very short fragments, which may be noisy, extremely
long fragments, which are likely to be assigned to the wrong bin, and
also the prospect of performing in silico size selection, for instance
to enrich cell-free DNA for short fragments which are known to be more
likely to come from the tumour.minReferenceLength - this parameter controls the
minimum distance along the reference genome which a fragment must cover
to be considered. This helps filter out some R1s which have only a small
portion of the read mapped correctly, and may help with regions with
strange artifactual binding (e.g. large insertions), but it is unlikely
to make much difference when properPairsOnly = TRUE.As part of this process, a number of extra variables are added to the
libraries slot in the qseaSet; here we will
explain what these are. If properPairsOnly = TRUE, the
r1 columns will not be present.
The filtering that is undertaken is as follows, if
properPairsOnly = TRUE:
Rsamtools. This is given by
qsea_initial_reads.qsea_initial_pairs; all other
reads are discarded.minMapQual
threshold. This therefore allows properly paired reads where one end is
multi-mapping. This step is reliant on the presence of the
MQ tags on the bam file, which are added by
samtools fixmate during bam file generation. This will be
stored as qsea_mapq_filtered_pairs.minInsertSize and maxInsertSize; this number
is then stored as qsea_size_filtered_pairs. This is also
stored as total_fragments (used in qsea).chrSelect with
the regions in badRegions removed, to give the
valid_fragments. This is the final number of fragments of
DNA that lie within the windows under consideration.When properPairsOnly = FALSE, the first read will be
retained from any read pair that is not considered “properly paired”,
alongside any R1s whose mate did not map or has been otherwise filtered
out. The rescuing process for proper pairs will not be performed (and
they will be kept as singleton R1s), and these R1s will be filtered for
MAPQ and presence in the windows only. There will be additional columns
qsea_initial_r1s and qsea_filtered_r1s, while
total_fragments and valid_fragments will
include both types of fragment.
Note that the coverageMethod = "PairedAndR1s" option
does use more RAM than the qsea method, which can be
selected via coverageMethod = "qseaPaired"; single-end
reads are not currently supported.
“Proper pairs”, as determined by the bwa mem aligner at
any rate, are paired-end sequencing reads where both ends have mapped
sufficiently near to each other, and they point towards each other on
opposite strands of the chromosome. However, it should be noted that how
near is near enough to be called as proper pairs is calculated within
bwa mem on an individual sample basis. This allowable
distance is calculated by determining the mean and standard deviation of
the fragment length amongst the mapped reads in that sample, with a
cutoff of 4 standard deviations away from the mean by default. In
general, for sheared genomic DNA this is not a big issue, because the
distribution is sufficiently broad that the standard deviation is large
enough. However, for circulating cell-free DNA, the distribution of
fragment sizes has a very sharp peak at ~167bp, and hence the global
fragment standard deviation is quite low. This can lead to read pairs of
length ~400bp which otherwise look like they should be “proper pairs”
not being called as such by the aligner. We therefore perform a rescuing
process, to ensure that we keep these read pairs by finding any read
pairs that meet the orientation conditions for being proper pairs and
are within minInsertSize to maxInsertSize in
length. This process is computationally expensive, so we only perform it
when properPairsOnly = TRUE. Note that bwa mem
has an option, -I, which allows you to specify the insert
size distribution yourself. This is highly recommended when working with
cell-free DNA.
The copy number variation is calculated over a set of larger genomic
windows as set by CNVwindowSize; this is typically much
larger than the windowSize due to coverage limitations.
This will be performed using the same options as the enriched samples.
Once the coverage is determined, the R package HMMcopy will
be used, with default parameters, to calculate the copy number
variation. The non-enriched (also known as “Input”) bam files should be
specified by adding a column called input_file to the
sampleTable. To use this you need to provide
HMMcopy with a pair of data frames containing the GC
content and mappability per region you wish to calculate over, as
detailed in vignette("HMMcopy", package = "HMMcopy"). There
is an exception if you are using the human genome hg38 (or
GRCh38), in which case files are provided for window sizes
of 50kb, 500kb or 1000kb. These
are stored in objects called e.g. gc_hg38_1000kb and
map_hg38_1000kb; similar objects need to be generated for
other genomes or window sizes using the tools gcCounter or
mapCounter from the HMMcopy
suite.
Alternatively, copy number variation may be calculated using an
external method, for instance using the ichorCNA
package. This may then be added to the qseaSet object using the
qsea function addCNV.
As a further option, qsea provides a method for
determining the copy number from the methylation enriched samples, using
only reads which do not contain a CG; this may be selected using
CNVmethod = "MeCap". We do note that for our MBD-seq
experiments we found there was reasonable overall agreement, but this
method did give spurious CNV compared to the non-enriched reads,
although that may be due to the selectivity of our MBD-seq enrichment
step being relatively high.
No copy number variation will be calculated if
CNVmethod = "None".
One essential part of the qseaSet is the
CpG_density for each region. This is calculated using the
genome sequence, alongside an expected fragment size distribution via
the qsea function
qsea:::estimatePatternDensity().
This process works by first locating all CG
dinucleotides across the genome, before modelling a distribution around
each, representing the likelihood for a fragment including that
dinucleotide to reach adjacent positions. This distribution is assumed
to be normally distributed, with a mean and standard deviation. Thus
each CG contributes to the density of the base pairs around
it. Finally, these values are averaged across the windows to give a
CpG_density column on the regions of the
qseaSet. This method is designed to capture the fact that
some captured fragments of DNA may be heavily methylated at one end, but
if the fragment is long it may be counted in an adjacent window.
This calculation therefore requires the fragmentLength
and fragmentSD parameters to be specified, to define these
normal distributions, and is associated with the whole
qseaSet. This is therefore not sample dependent, and so if
samples have vastly different fragment distributions then there will be
some expected variation not captured in this process. This does have a
significant effect on the calculated beta values for each sample, so it
is important to set the fragmentLength and
fragmentSD values in an appropriate way. We also provide a
fragmentType option, which hardcodes these two parameters
to values derived from our initial sheared or cfDNA datasets. When
fragmentType = "Sheared" is used this sets
fragmentLength = 213, fragmentSD = 60, while
fragmentType = "cfDNA" sets
fragmentLength = 167, fragmentSD = 38, but we do emphasise
that you should check your own data. Ideally the
CpG_density would be a property of the fragment
distribution of each individual sample, but mesa does not
support this at this time.
Please note that these fragmentLength and
fragmentSD parameters are only used for calculating this CG
density metric; they are not used to control what reads are included in
the coverage calculations.
In any ME-seq experiment, there will be off-target fragments
sequenced due to the imperfect nature of the enrichment process. These
are used to calculate the offset; a measure of the
background level of fragments in windows that do not contain CGs. The
maximum CpG_density value for these windows is controlled
by maxPatternDensity. If this is set too low, then
qsea will give a warning saying that the offset cannot be
calculated. In particular, this could be an issue if using a targeted
enrichment step in your library preparation method.
Finally, the qsea beta value calculation is applied.
This requires the CpG_density and the coverage, as well as
some assumption about which regions are methylated, or at least how many
of them are. The qsea
manuscript details three methods for this assumption, two of which
require external data in the form of methylation arrays or
bisulphite/enzymatic sequencing. The third method, known as “blind
calibration” in the qsea manuscript, is the only one implemented here.
This works on a global assumption about the typical methylation levels
across the human genome: that the majority of isolated CGs are
methylated, while the majority of CpG islands are unmethylated. This
blind calibration method therefore assumes an average methylation level
of 76% in windows with a CpG_density of 1, falling along a
straight line to 25% in windows with a CpG_density of 15.
This method is the only one implemented here, as it does not require
additional data.
We have found this blind calibration sufficient for the human genome, including being able to build classifiers using methylation arrays as our data source.
If you wish to use array data for normalisation, then this Bioconductor support post may be a good source of information, but we have only used the blind calibration method.
One thing to note is that this function does not set the
library_factor for the qseaSet, which
qsea suggests doing via the Trimmed Mean of M-values (TMM)
method. The way qsea does this is to choose one sample (by
default the first sample), and use this as the reference sample to
calculate a scaling factor to adjust the total fragments by when
calculating the NRPM (normalised reads per million) values.
TMM normalisation is a method used in RNA-seq to compensate for genes which are excessively expressed in a particular condition. In RNA-seq, the distribution of read counts varies over ~4 orders of magnitude, with some highly expressed genes having tens of thousands of reads while others have tens or hundreds. If in a particular tissue there is a group of genes which are extremely highly expressed that are not expressed in another tissue, then the total pool of available reads will be dominated by those genes in the first tissue, comparatively depressing the other genes. So just using the total size ends up biased towards these highly expressed genes, making it look like the other genes are less expressed, and giving many more differentially expressed genes than might actually be the case.
For example, if there are 100 reads in a sample, and 3 genes, with values:
then we want to be careful not to call the last two genes as differentially expressed, when they are really just affected by the lack of the first gene. One option could be to normalise with the median of the read distribution, although, as many genes are not expressed, this is often implemented as an upper quartile.
TMM normalisation however takes some subset of the “genes” (windows in our case), having trimmed off the very high/very low expressed “genes”, to calculate this scaling factor for total number of fragments.
However, we find that for methylation, after filtering for excluded regions, extremely few windows in any sample have more than 100 reads. So the distribution of read counts only spans roughly 0-100, removing the main reason for TMM normalisation.
There are a couple of reasons why we therefore choose not to apply TMM:
calculateDMRs calculation), there is no effect
of the library_factor, as it is normalised out in the beta
value fitting process.qsea).
This therefore means that the normalisation of all the samples depends
on which other samples they were normalised with, meaning that comparing
between batches would need renormalising every time. This would also
mean that filtering samples pre- or post-normalisation would give
different results.However, we note that these considerations might not apply for
everyone; the qsea::addLibraryFactors() function may be
used to apply this normalisation if required.
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 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=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: UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] BSgenome.Hsapiens.UCSC.hg19_1.4.3 BSgenome_1.80.0
#> [3] rtracklayer_1.72.0 BiocIO_1.22.0
#> [5] Biostrings_2.80.2 XVector_0.52.0
#> [7] GenomicRanges_1.64.0 Seqinfo_1.2.0
#> [9] IRanges_2.46.0 S4Vectors_0.50.3
#> [11] BiocGenerics_0.58.1 generics_0.1.4
#> [13] tibble_3.3.1 stringr_1.6.0
#> [15] mesa_0.99.8.9000 dplyr_1.2.1
#> [17] qsea_1.38.0
#>
#> loaded via a namespace (and not attached):
#> [1] SummarizedExperiment_1.42.0 rjson_0.2.23
#> [3] xfun_0.61 bslib_0.12.0
#> [5] plyranges_1.32.0 htmlwidgets_1.6.4
#> [7] Biobase_2.72.0 lattice_0.23-1
#> [9] vctrs_0.7.3 tools_4.6.1
#> [11] bitops_1.1-0 curl_8.1.0
#> [13] parallel_4.6.1 pkgconfig_2.0.3
#> [15] Matrix_1.7-6 data.table_1.18.6.1
#> [17] desc_1.4.3 cigarillo_1.2.1
#> [19] lifecycle_1.0.5 compiler_4.6.1
#> [21] Rsamtools_2.28.0 textshaping_1.0.5
#> [23] statmod_1.5.2 codetools_0.2-20
#> [25] MEDIPSData_1.48.0 GenomeInfoDb_1.48.0
#> [27] htmltools_0.5.9 sass_0.4.10
#> [29] RCurl_1.98-1.20 yaml_2.3.12
#> [31] pillar_1.11.1 pkgdown_2.2.1
#> [33] crayon_1.5.3 jquerylib_0.1.4
#> [35] BiocParallel_1.46.0 limma_3.68.5
#> [37] cachem_1.1.0 DelayedArray_0.38.2
#> [39] abind_1.4-8 gtools_3.9.5
#> [41] tidyselect_1.2.1 digest_0.6.39
#> [43] stringi_1.8.9 purrr_1.2.2
#> [45] restfulr_0.0.17 fastmap_1.2.0
#> [47] grid_4.6.1 cli_3.6.6
#> [49] SparseArray_1.12.3 magrittr_2.0.5
#> [51] S4Arrays_1.12.1 XML_3.99-0.25
#> [53] withr_3.0.3 UCSC.utils_1.8.0
#> [55] rmarkdown_2.32 httr_1.4.9
#> [57] matrixStats_1.5.0 otel_0.2.0
#> [59] ragg_1.5.2 zoo_1.9-1
#> [61] memoise_2.0.1 evaluate_1.0.5
#> [63] knitr_1.52 rlang_1.3.0
#> [65] glue_1.8.1 HMMcopy_1.54.0
#> [67] jsonlite_2.0.0 R6_2.6.1
#> [69] MatrixGenerics_1.24.0 GenomicAlignments_1.48.0
#> [71] systemfonts_1.3.2 fs_2.1.0