---
title: "SingleMoleculeGenomicsIO"
output:
    rmarkdown::html_document:
        toc: true
        toc_float: true
vignette: >
  %\VignetteIndexEntry{SingleMoleculeGenomicsIO}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# Overview

`SingleMoleculeGenomicsIO` provides tools for reading and processing 
single-molecule genomics data in R. These include functions for reading base 
modifications from `modBam` files (standard bam files with modification data 
stored in `MM` and `ML` tags, see
[SAMtags.pdf](https://samtools.github.io/hts-specs/SAMtags.pdf)), 
e.g. generated from Oxford Nanopore or PacBio data using tools like 
[Dorado](https://github.com/nanoporetech/dorado) or equivalent, 
from text files generated by
[modkit](https://nanoporetech.github.io/modkit), or from `bam` files where 
modification information is encoded via sequence mismatches (e.g., resulting
from bisulfite sequencing or deamination-based protocols).
The package also provides an efficient representation of such data as R objects
(based on the widely used `SummarizedExperiment` class), and functions to 
manipulate and summarize such objects.
`SingleMoleculeGenomicsIO` does not provide functions for downstream analyses
(e.g. differential footprint identification or nucleosome analyses), or
visualization. Such functionality can for example be found in
[footprintR](https://fmicompbio.github.io/footprintR/) or
[nomeR](https://github.com/fmi-basel/compbio-nomeR).

In this vignette, we will focus on reading data from these different formats
and outline how to process the obtained data, e.g. to filter out 
low-quality reads, subset to specific sequence contexts, or regroup the reads
according to some annotation of interest.

Current contributors include:

- [Panagiotis Papasaikas](https://github.com/ppapasaikas)
- [Sebastien Smallwood](https://github.com/smalseba)
- [Charlotte Soneson](https://github.com/csoneson)
- [Michael Stadler](https://github.com/mbstadler)

# Installation

`SingleMoleculeGenomicsIO` can be installed from GitHub via the 
`BiocManager` package:

```
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("fmicompbio/SingleMoleculeGenomicsIO")
```

# Load required packages

```{r load-pkgs, message=FALSE}
library(SingleMoleculeGenomicsIO)
library(SummarizedExperiment)
library(Biostrings)
library(SparseArray)
library(BiocParallel)
```

# Data representation

The figure below illustrates the basic structure of the data container used 
within `SingleMoleculeGenomicsIO`. All reading functions described in the 
next section generate a container of this type.

![](SE-read-summary-level.png)

As we see, the data is stored in a `RangedSummarizedExperiment` object. This is 
a standard Bioconductor container that can contain one or more assays (storing 
the observed data), together with annotations for rows (genomic positions) and 
columns (samples or individual reads), as well as general metadata in a list 
format. 

The assay data in the object can be either on the individual read level, or 
represent data summarized across all reads for a given sample. As we will see 
below, an object can contain both read-level and summary-level data. 

# Load data from `modBam` files

In this section we will mainly focus on reading data from `modBam` files, and 
explore the generated object. For more information about reading other types of 
supported files (see above), we refer to section \@ref(other-formats). 
As mentioned, regardless of the reading function that is used, the resulting
`SummarizedExperiment` object will have the same format. 

## Load read-level data

A small example dataset with 6mA modifications representing
accessibility, in the form of two `modBam` files generated using the
[Dorado](https://github.com/nanoporetech/dorado) aligner, is provided with 
`SingleMoleculeGenomicsIO`. Here, we read this into R using the `readModBam()` 
function. Typically, only a subset of the full `modBam` file, corresponding 
to a specific genomic region of interest, is loaded. 

As these `modBam` files represent 6mA modifications, we set `modbase = "a"`
in the `readModBam()`. For a full list of
possible modifications, see the section about base modifications in the
[SAM tag specification](https://samtools.github.io/hts-specs/SAMtags.pdf)). 
We also provide a sequence reference (genome sequence) from which to extract 
the genomic sequence context for modified bases. This can be helpful during 
analysis, e.g. to filter out observations that are likely to represent 
sequencing error. 


```{r data-files}
# read-level 6mA data generated by Dorado
modbamfiles <- system.file("extdata",
                           c("6mA_1_10reads.bam", "6mA_2_10reads.bam"),
                            package = "SingleMoleculeGenomicsIO")
names(modbamfiles) <- c("sample1", "sample2")

# file with sequence of the reference genome in fasta format
reffile <- system.file("extdata", "reference.fa.gz",
                       package = "SingleMoleculeGenomicsIO")
```


```{r read-data}
se <- readModBam(bamfiles = modbamfiles,
                 regions = "chr1:6940000-7000000",
                 modbase = "a",
                 verbose = TRUE, 
                 sequenceContextWidth = 1, 
                 sequenceReference = reffile,
                 BPPARAM = SerialParam())
se
```

This will create a `RangedSummarizedExperiment` object with rows corresponding 
to individual genomic positions:

```{r row-data}
# rows are positions...
rowRanges(se)
```

The positions are stranded, as base modifications can be observed on either of 
the two strands:

```{r tabulate-strand}
table(strand(se))
```

We notice that (as expected) most of the positions for which we obtain a
modification value indeed correspond to an 'A' in the genome sequence. 
Discrepancies here may arise from, e.g., sequencing errors or differences 
between the genome sequence of the studied subjects and the reference sequence.

```{r tabulate-seqcontext}
table(rowRanges(se)$sequenceContext)
```

The columns of `se` correspond to samples (recall that in our case we loaded 
two `modBam` files):

```{r column-data}
# ... and columns are samples
colData(se)
```

The sample names in the `sample` column are obtained from the input files
(here `modbamfiles`), or if the files are not named will be automatically
assigned (using `"s1"`, `"s2"` and so on, assuming that each file corresponds
to a separate sample).

Read-level information is contained in a list under `readInfo`, with one
element per sample:

```{r coldata-read-info}
se$readInfo
```

In general, columns in `colData` that provide read-level information can be
obtained using `getReadLevelColDataNames()`:

```{r coldata-read-level-colnames}
getReadLevelColDataNames(se)
```

Each sample typically contains several reads (see `se$n_reads`),
which correspond to the rows of `se$readInfo`, for example for `"sample1"`:

```{r coldata-read-info-sample1}
se$readInfo$sample1
```

## Explore assay data

At this point, the `RangedSummarizedExperiment` object contains a single assay 
named `mod_prob`, which is a `DataFrame` with modification probabilities.

```{r mod-prob-assay}
assayNames(se)

m <- assay(se, "mod_prob")
m
```

In order to store read-level data for variable numbers of reads per sample,
each sample column is not just a simple vector, but a position-by-read
`NaMatrix` (a sparse matrix object defined in the
[`SparseArray`](https://bioconductor.org/packages/SparseArray) package).

```{r mod-prob-columns1}
m$sample1
```

If a 'flat' matrix is needed, in which columns correspond to
reads instead of samples, it can be created using `as.matrix()`.
Notice the number of columns in each of the objects:

```{r mod-prob-columns2}
ncol(se)           # 2 samples
ncol(m$sample1)    # 3 reads in "sample1"
ncol(m$sample2)    # 2 reads in "sample2"
ncol(as.matrix(m)) # 5 reads in total
```

In general, assays that provide read-level information can be
obtained using `getReadLevelAssayNames()`:

```{r se-read-level-assaynames}
getReadLevelAssayNames(se)
```

One advantage of the grouping of reads per sample is that we can easily
perform per-sample operations using `lapply` (returns a `list`) or `endoapply`
(returns a `DataFrame`):

```{r mod-prob-lapply}
lapply(m, ncol)
endoapply(m, ncol)
```

The `NaMatrix` objects do not store `NA` values explicitly and thus use much
less memory compared to a normal (dense) matrix. The `NA` values here occur
at positions (rows) that are not covered by the read of that column, and there
are typically a large fraction of NA values:

```{r nnacount}
prod(dim(m$sample1)) # total number of values
nnacount(m$sample1)  # number of non-NA values
```

**Important: ** These `NA` values have to be excluded for example when
calculating the average modification probability for a position
using `mean(..., na.rm = TRUE)`, otherwise the result will be `NA` for
all incompletely covered positions:

```{r mod-prob-zeros}
# find row corresponding to position "chr1:6928850:-"
(idx <- which(seqnames(se) == "chr1" & start(se) == 6928850 & strand(se) == "-"))

# modification probabilities at the given position
m[idx, ]

# WRONG: take the mean of all values (including NAs)
lapply(m[idx, ], mean)

# CORRECT: exclude the NA values (na.rm = TRUE)
lapply(m[idx, ], mean, na.rm = TRUE)
```

Note that for sample2, in which no read overlaps the selected position,
`mean(..., na.rm = TRUE)` returns `NaN`.

In most cases it is preferable to use the convenience functions in
`SingleMoleculeGenomicsIO`, such as `flattenReadLevelAssay()` (see next
section), to calculate summaries over reads instead of the `lapply(...)` used
above for illustration.

## Summarize read-level data

Summarized data can be obtained from the read-level data by calculating
a per-position summary over the reads in each sample, and are stored as 
additional assays in the resulting object. The `flattenReadLevelAssay()` 
function can calculate several relevant summaries:

```{r summarise-read-level}
se_summary <- flattenReadLevelAssay(
    se, statistics = c("Nmod", "Nvalid", "FracMod", "Pmod", "AvgConf"))
# new assays are added to the object
assayNames(se_summary)
```

As discussed above, this will correctly exclude non-observed values
and return zero values (for count assays `"Nmod"`, "`Nvalid`") or `NaN` for
assays with fractions (`"FracMod"`, `"Pmod"` and `"AvgConf"`), for positions
without any observed data (see for example `"Pmod"` in `"sample2"`):

```{r summarized-pmod}
assay(se_summary, "Pmod")[idx, ]
```

The summary statistics to calculate are selected using the `statistics` argument.
By default, `flattenReadLevelAssay()` will count the number of modified (`Nmod`)
and total (`Nvalid`) reads at each position and sample, and calculate the
fraction of modified bases from the two (`FracMod`).

```{r default-statistics}
assay(se_summary, "Nmod")[idx, ]
assay(se_summary, "Nvalid")[idx, ]
assay(se_summary, "FracMod")[idx, ]
```

In the above example, we additionally calculate the average modification
probability (`Pmod`) and the average confidence of the modification calls
per position (`AvgConf`). As we have not changed the default `keepReads = TRUE`,
we keep the read-level assay from the input object (`mod_prob`) in the returned
object. Because of the grouping of reads by sample in `mod_prob`, the dimensions
of read-level or summary-level objects remain identical. However, as the 
summary-level assays only contain a single value per position and sample, they
are regular matrices rather than the hierarchical `DataFrames` used for the 
read-level data.

```{r collapsed-assays}
# read-level data is retained in "mod_prob" assay
assayNames(se_summary)

# ... which groups the reads by sample
assay(se_summary, "mod_prob")

# the dimensions of  read-level `se` and summarized `se_summary` are identical
dim(se)
dim(se_summary)

# summary-level assays are regular matrices
class(assay(se_summary, "FracMod"))
```


## Load summary-level data

As illustrated above, summary-level data (e.g. average modification fraction at
each position) can be calculated from read-level data. Alternatively, reads can
also be directly summarized while reading from `modBam` files, which is faster
and reduces the size of the object in memory. In addition, there are other file
formats that only store summary-level data (e.g. `bed` files) that can also be
imported using `SingleMoleculeGenomicsIO` (see below).

To read summary-level data with `readModBam()`, we set the `level` argument
to `"summary"`:

```{r read-data-summary}
se_summary2 <- readModBam(bamfiles = modbamfiles,
                          regions = "chr1:6940000-7000000",
                          modbase = "a",
                          level = "summary",
                          verbose = TRUE, 
                          BPPARAM = SerialParam())
se_summary2
```

Notice that `se_summary2` contains only summary-level assays (`Nmod`, `Nvalid` 
and `FracMod`), and no read-level modifications anymore (assay `mod_prob`).
Otherwise, it resembles the read-level object `se`, for example it covers the
same genomic positions:

```{r row-data-summary}
# rows are positions...
rowRanges(se)
rowRanges(se_summary2)
```

The summary-level assays generated from `se` (see `se_summary` above) are also
identical to the ones directly read from the `modBam` file into `se_summary2`:

```{r identical-summaries}
identical(assay(se_summary, "Nvalid"),
          assay(se_summary2, "Nvalid"))
```

## Load data for a random subset of the reads

By default, `readModBam()` will load all reads overlapping the regions 
indicated by the `regions` argument. However, it is also possible to load a 
random subset of reads (this can be useful e.g. to calculate distributions 
of quality scores, using reads aligning to different parts of the genome but 
without the need to load the entire `modBam` file into R). This is 
achieved by setting the `nAlnsToSample` argument to a non-zero value 
(representing the desired number of reads per sample), and indicating which 
reference sequences to sample reads from. 

```{r sample-reads}
se_sample <- readModBam(bamfiles = modbamfiles, 
                        modbase = "a", 
                        nAlnsToSample = 5, 
                        seqnamesToSampleFrom = "chr1", 
                        verbose = TRUE, 
                        BPPARAM = SerialParam(RNGseed = 1327828L))
se_sample$n_reads
```

Note that the exact number of sampled reads may not be identical to the 
requested number. One reason for this is that the requested number is combined
with the total number of alignments to calculate the fraction of reads to be 
sampled (in this small case, 0.5), and each alignment will be retained with 
this probability. 

## Label reads by sequence variants

The `readModBam` and `readMismatchBam` functions allow to automatically label
reads using variant positions (`variantPositions` argument):

```{r read-varpositions}
# create GPos object with reference sequence variant positions
varpos <- GPos(seqnames = "chr1", pos = c(6937731, 6937788, 6937843, 6937857,
                                          6937873, 6937931, 6938932, 6939070))

# load and label reads at varpos
se_varpos <- readModBam(bamfiles = modbamfiles, 
                        modbase = "a", 
                        regions = "chr1:6925000-69500000",
                        variantPositions = varpos,
                        verbose = FALSE, 
                        BPPARAM = SerialParam(RNGseed = 1327828L))
```

Read bases at `variantPositions` will be extracted and concatenated to create
a variant label for each read, using `'-'` for non-overlapping positions.
The resulting labels are stored for each sample in the `variant_label` column
of `readInfo`:

```{r readlabels-sample1}
se_varpos$readInfo$sample1$variant_label
```

In a diploid experimental system with known heterozygous positions, this allows
for example to measure allele specific signals by grouping reads based on these
variant labels. A step-by-step example illustrating this procedure is
available in the [footprintR book](https://fmicompbio.github.io/footprintR_book/grouping-reads.html#sec-group-by-sequence).

# Load data from other file types {#other-formats}

In the previous sections we illustrated how to read data from `modBam` files
and extract both read-level and summary-level information. As mentioned in the 
beginning of this document, `SingleMoleculeGenomicsIO` can also be used to read
data from other file formats that are often used to store data from single 
molecule genomics experiments. 

In cases where read-level modification data has already been extracted from the
`modBam` files using [modkit](https://nanoporetech.github.io/modkit), these 
files can be read into a `RangedSummarizedExperiment` using the 
`readModkitExtract()` function.

```{r read-modkit-extract}
# define path to modkit extract output file
extrfile <- system.file("extdata", "modkit_extract_rc_5mC_1.tsv.gz",
                        package = "SingleMoleculeGenomicsIO")

# read complete file
readModkitExtract(extrfile, modbase = "m", filter = NULL,
                  sequenceContextWidth = 3,
                  sequenceReference = reffile, 
                  BPPARAM = SerialParam())
# only read entries that are not suggested to be filtered out
readModkitExtract(extrfile, modbase = "m", filter = "modkit",
                  sequenceContextWidth = 3,
                  sequenceReference = reffile, 
                  BPPARAM = SerialParam())
```

Summary-level data can also be provided as `bedMethyl` files (containing
total and methylated counts for individual positions). These can be read using
the `readBedMethyl()` function, which again returns a 
`RangedSummarizedExperiment` object (this time only containing summary-level 
assays).

```{r read-data-bedmethyl}
# collapsed 5mC data for a small genomic window in bedMethyl format
bedmethylfiles <- system.file("extdata",
                              c("modkit_pileup_5mC_1.bed.gz",
                                "modkit_pileup_5mC_2.bed.gz"),
                              package = "SingleMoleculeGenomicsIO")

names(bedmethylfiles) <- c("sample1", "sample2")
readBedMethyl(bedmethylfiles,
              modbase = "m",
              sequenceContextWidth = 3,
              sequenceReference = reffile, 
              BPPARAM = SerialParam())
```

In addition to the `modBam` files, where modification information is 
encoded in the `MM` and `ML` tags, `SingleMoleculeGenomicsIO` can also load 
data from `bam` files where modification information has to be inferred from 
the observed read sequence. This includes bisulfite sequencing data (where 
unmethylated and typically inaccessible Cs are converted to Ts) and data
generated using deamination protocols (where instead the accessible Cs are
converted to Ts). Because of  these conversions, such data is aligned to the
reference genome using dedicated alignment tools like the ones provided with
the `QuasR` R package or the `Bismark` aligner. `SingleMoleculeGenomicsIO`
supports reading `bam` files generated by either of these, using the
`readMismatchBam()` function. 
The genomic sequence context to focus on, as well as specifications about 
which bases represent modified and unmodified states, respectively, are 
determined based on the arguments `sequenceContext`, `readBaseUnmod` and 
`readBaseMod`.

```{r read-data-mismatchbam}
# bisulfite sequencing bam file aligned with QuasR::qAlign
quasrbisfile <- system.file("extdata", "BisSeq_quasr_paired.bam",
                            package = "SingleMoleculeGenomicsIO")

# read data
readMismatchBam(bamfiles = quasrbisfile, 
                bamFormat = "QuasR", 
                regions = "chr1:6950000-6955000", 
                sequenceContext = "C", 
                sequenceReference = reffile,
                readBaseUnmod = "T", 
                readBaseMod = "C", 
                BPPARAM = SerialParam())

# bisulfite sequencing bam file aligned with Bismark
bismarkbisfile <- system.file("extdata", "BisSeq_bismark_paired.bam",
                              package = "SingleMoleculeGenomicsIO")

# read data
readMismatchBam(bamfiles = bismarkbisfile, 
                bamFormat = "Bismark", 
                regions = "chr1:6950000-6955000", 
                sequenceContext = "C", 
                sequenceReference = reffile,
                readBaseUnmod = "T", 
                readBaseMod = "C", 
                BPPARAM = SerialParam())
```

The table below compares the capabilities of the provided reading functions. 

| Function | Description | Read-level output | Summary-level output | Region-based access | Sampling | Variant read labels |
| :------: | :---------: | :---------------: | :------------------: | :-----------------: | :------: | :-----------------: |
| `readModBam` | Read `modBam` file(s) | Yes | Yes | Yes | Yes | Yes |
| `readModkitExtract` | Read `modkit` extract file(s) | Yes | No | No | No | No |
| `readBedMethyl` | Read `bedMethyl` file(s) | No | Yes | No | No | No |
| `readMismatchBam` | Read `bam` file(s) with modifications encoded as as sequence mismatches | Yes | Yes | Yes | Yes | Yes |


# Expand the `SummarizedExperiment` object to single base resolution

The objects returned by the reading functions above contain data only for 
the positions with an observable modification state. For
cases where single nucleotide resolution is required (i.e., where the rows 
in the object must cover each position within a region), the 
`expandSEToBaseSpace` will return a suitable object. Note that any strand 
information in the original positions will be ignored.

```{r expand-se}
se_expand <- expandSEToBaseSpace(se)
dim(se_expand)
```

# Filter positions

`SingleMoleculeGenomicsIO` provides several ways to filter out individual 
genomic positions, based on e.g. the genomic sequence context or the 
observed read coverage. Such filtering can be important to reduce noise in 
the data for downstream analyses. 

To illustrate the position-level filtering, let us consider the `se` object 
obtained from the two `modBam` files above: 

```{r show-se-before-filtering}
se
```

Recall that we added the genomic sequence context of each position when reading 
this data (in `rowRanges(se)$sequenceContext`). We could also add this manually
to the `RangedSummarizedExperiment`, via the `addSeqContext()` function. 
We use the `filterPositions()` function to filter the rows of `se`, retaining 
only rows where the nucleotide in the corresponding genomic position is
an `A`, and where the read coverage is at least 2. Note that for the sequence 
context filter, any IUPAC code can be used - for example, `NCG` would match 
trinucleotide contexts with a central C followed by a G (and preceded by any 
base).


```{r filter-positions}
se_filt <- filterPositions(
    se, 
    filters = c("sequenceContext", "coverage"),
    sequenceContext = "A",
    assayNameCov = "mod_prob",
    minCov = 2
)
dim(se_filt)
```

In this case, the filtering reduced the number of unique positions from 
`r nrow(se)` to `r nrow(se_filt)` (a reduction by 
`r 100*(1 - signif(nrow(se_filt)/nrow(se), 3))`%). 
At the same time, the number of non-NA 
values in the `mod_prob` assay is only reduced from 
`r length(nnawhich(as.matrix(assay(se))))` to 
`r length(nnawhich(as.matrix(assay(se_filt))))` (a reduction by 
`r 100*(1 - signif(length(nnawhich(as.matrix(assay(se_filt))))/length(nnawhich(as.matrix(assay(se)))), 3))`%), indicative of the fact that 
the filtered out positions are generally covered by only a few reads. 

# Filter reads

In addition to the position filtering illustrated above, 
`SingleMoleculeGenomicsIO` provides functionality for calculating read-level 
quality scores and filtering out entire reads with low quality. The 
calculation of the quality scores is performed using the `addReadStats()` 
function, which will return a `RangedSummarizedExperiment` object with the 
quality scores added as a new column in the `colData` (by default named `QC`).

```{r add-read-stats}
# start from read-level data
# ... some read-level statistics were added already when reading the data
se$readInfo$sample1

# ... calculate additional read statistics
se <- addReadStats(se, name = "QC", 
                   BPPARAM = SerialParam())
se$QC$sample1
```

As a result, the `colData` now contains two columns with read-level information:

```{r readlevel-coldata-after-qc}
getReadLevelColDataNames(se)
```

The distribution of the quality statistics can also be plotted, which is often 
helpful in order to decide on cutoffs for downstream analysis. Note that in the 
current example data, each sample contains only a few reads and hence the plot
below is very sparse. 

```{r plot-read-stats}
plotReadStats(se, readInfoCol = "readInfo", qcCol = "QC")
```

The `filterReads()` function can be used to filter the set of reads based on 
various quality indicators. This step requires that the read statistics that 
will be used for the filtering have been calculated as indicated above.

```{r filter-reads}
se_filtreads <- filterReads(se, readInfoCol = "readInfo",
                            qcCol = "QC", minQscore = 10, 
                            minAlignedFraction = 0.8)
# number of retained reads for each sample
lapply(assay(se_filtreads, "mod_prob"), ncol)
# tables with the reason behind each filtered read
metadata(se_filtreads)$filteredOutReads
```

Both position and read filters above were applied after reading the data from 
the `modBam` file into a `RangedSummarizedExperiment` object. An alternative to 
this is to create a filtered `modBam` file once, and then use that for all 
downstream analyses. This can be done using the `filterReadsModBam()` file, 
which accepts the same filter criteria as `filterReads()` above. Here we 
illustrate the process using one of the example `modBam` files provided with 
the package. 

```{r filter-reads-modbam}
# input file
modbamfile <- system.file("extdata", "6mA_1_10reads.bam", 
                          package = "SingleMoleculeGenomicsIO")

# output file
filtmodbamfile <- tempfile(fileext = ".bam")

res <- filterReadsModBam(infiles = modbamfile, 
                         outfiles = filtmodbamfile, 
                         modbase = "a", 
                         indexOutfiles = FALSE, 
                         minQscore = 10, 
                         minAlignedFraction = 0.8,
                         BPPARAM = BiocParallel::SerialParam(),
                         verbose = TRUE)
res
```


# Subset and regroup reads

After importing read-level data with `SingleMoleculeGenomicsIO`, the set of 
reads can be subset using the `subsetReads()` function. The subsetting can be 
done either by specifying the read names or indices to retain, or by asking 
for a random subset of reads. 

```{r subset-reads}
# subset to only the first two reads per sample
(reads_to_keep <- lapply(assay(se, "mod_prob"), function(s) 1:2))
se_subset <- subsetReads(se, reads = reads_to_keep)
lapply(assay(se_subset, "mod_prob"), ncol)

# subset to specified reads
(reads_to_keep <- unlist(lapply(assay(se, "mod_prob"), function(s) colnames(s)[c(1, 2)])))
se_subset <- subsetReads(se, reads = reads_to_keep)
lapply(assay(se_subset, "mod_prob"), ncol)

# request a random set of two reads per sample
se_subset <- subsetReads(se, randomSubset = 2)
lapply(assay(se_subset, "mod_prob"), ncol)
```

See `?subsetReads` for more subsetting options. 

In some situations, a _regrouping_ of the reads is more suitable than a 
subsetting. By default, the reads will be grouped by the sample they originate 
from (i.e., each column in the `SummarizedExperiment` corresponds to a sample). 
The `regroupReads()` and `regroupReadsByColData()` functions allow us to 
reorganize the reads in such a way that each column of the object instead 
corresponds to another read property (e.g., the strand from which the 
read originates). The new groups can be defined manually as a named list 
of read identifiers (for `regroupReads()`), or be generated automatically from 
a categorical column in the read-level column annotations (for 
`regroupReadsByColData()`). Here we illustrate how to regroup the reads by 
the strand. 

```{r regroup-reads}
# read strand is stored in the column annotations (ref_strand)
se$readInfo$sample1

# regroup reads by strand, ignoring the sample information
se_regroup <- regroupReadsByColData(se, colNames = "readInfo:ref_strand", 
                                    withinSample = FALSE)
colnames(se_regroup)

# regroup reads by strand, within each sample
se_regroup <- regroupReadsByColData(se, colNames = "readInfo:ref_strand", 
                                    withinSample = TRUE)
colnames(se_regroup)
```

# Add read segment annotations

Read segment annotation allow us to store information about specific subregions
of reads. These could, for example, represent inferred locations of nucleosomes
or transcription factor binding sites. The segments need to be provided as a 
list of `IRangesList` objects (one per sample), with each element of an 
`IRangesList` providing the annotations for one read. 

```{r annotate-segments}
# define segments
segs <- list(sample1 = IRangesList(
    "sample1-d52a5f6a-a60a-4f85-913e-eada84bfbfb9" = IRanges(
        start = 6925834, width = 140),
    "sample1-233e48a7-f379-4dcf-9270-958231125563" = IRanges(
        start = c(6926000, 6926200), width = 140)), 
    sample2 = IRangesList())

# add segment annotation to se
se <- annotateReadSegments(se, irlList = segs, name = "nucl")
se$nucl
```

# Count state pairs

The import functions exemplified above all return a `SummarizedExperiment` 
object with read- and/or summary-level data. For certain applications, 
other data representations are more suitable. The `countStatePairs()` 
function allows us to parse a `modBam` file and tabulate the number of pairs 
of bases at a given distance and with a given combination of modification 
states. 

```{r count-state-pairs}
sp <- countStatePairs(bamfile = modbamfiles[1], 
                      modbase = "a", 
                      windowSize = 200, 
                      BPPARAM = BiocParallel::SerialParam(), 
                      verbose = FALSE)
sp
```


# Session info

```{r session-info}
sessioninfo::session_info()
```