---
title: "Introduction to Mesa"
output: rmarkdown::html_vignette
vignette: >
    %\VignetteIndexEntry{Introduction to Mesa}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8}
---

```{r chunk, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.width = 7,
    fig.height = 5,
    out.width = "100%",
    dev = "ragg_png"
)
```

```{r setup}
library(mesa)
```


# Introduction
This vignette introduces the `mesa` package, which is designed to help analyse
methylation enrichment sequencing data. Methylation enrichment sequencing data
is generated by performing next-generation sequencing (NGS) on DNA which has
been enriched for those fragments of DNA which contain methylation on the
cytosine of CG dinucleotides. This may be, for example, via
immunoprecipitation (MeDIP-seq), or via a pulldown performed with a protein in
the methyl-binding-domain (MBD) family (MBD-seq).

Once this step has been performed, the enriched set of DNA fragments is then
sequenced to generate FASTQ files. These should be aligned to a reference genome
and deduplicated (this deduplication MUST be done, to prevent PCR duplicates
from looking like a peak) to produce bam files using standard alignment
techniques. For details on how to generate the qseaSet from bam files, see
`vignette("generation")`. If the enrichment step has worked correctly, there
should be peaks of reads present at regions of the genome where the DNA
contained methylated CGs, with minimal background reads elsewhere.

# Example Dataset
For the examples here, we will use a dataset consisting of 10 samples: 5 tumours
and 5 matched adjacent normal samples. For size reasons, this only includes a
small part of chromosome 7. As noted above, `vignette("generation")` describes
how to generate the qseaSet object.

First we will shorten the name of the example dataset to just be `qseaSet`, and
see what it contains. Throughout this vignette we assume familiarity with the
tidyverse `vignette("paper", package = "tidyverse")`, and with the pipe operator
`|>` to chain commands together.

```{r example-dataset}
library(qsea)
qseaSet <- exampleTumourNormal
```

As the underlying structure to hold our methylation data we use the `qsea`
package. The tutorial vignette for `qsea` is available via
`vignette("qsea_tutorial", package = "qsea")`, but we will cover some of
these concepts here.

The qseaSet object contains a number of different pieces of information
regarding the experiment, in one single object. If we evaluate it, we see some
summary information. This includes the number of samples included, their names
and how many windows are present.

```{r chunk-2}
qseaSet
```

# Parts of a qseaSet
The qseaSet contains multiple parts to it. The main ones to be aware of are:

* Sample table
* Counts matrix
* Windows
* Copy number variation (CNV) matrix
* Library information

We will go into some detail about each of these components now.

## Sample Table
The sample table is the sample metadata (data about the data). This is a data
frame with one row for each sample in your qseaSet. It *MUST* contain a column
called `sample_name`, as this is the primary identifier used in `qsea`, as well
as a `group` column which indicates replicate information to collapse across (if
present) for certain functions. Any other columns can then be added onto this
table, and it should contain any metadata to describe the different samples
contained within the qseaSet.

We access the sample table using the function `getSampleTable`:

```{r sample-table}
qseaSet |> 
    getSampleTable() |>
    as_tibble()
```

We can also get a vector of the samples present inside the qseaSet by using
`getSampleNames`:

```{r chunk-3}
qseaSet |> 
    getSampleNames()
```

As the sample table is so important, mesa extends a number of `dplyr`
(tidyverse) functions that let you modify it. For instance, if we want to add a
new column to the sample table, we can use `mutate`:

```{r chunk-4}
qseaSet |>
    mutate(new_column = "New") |>
    getSampleTable()
```

Sometimes you have further information about some or all samples that you want
to include, at which point you can use `left_join` to combine with a new data
frame. Missing data is fine, and you can join on a different column by using
`by`.

```{r chunk-5}
tumourVAFdata <- tribble(~sample_name, ~VAF, #VAF = Variant Allele Frequency
                            "Colon1_T",23,
                            "Colon2_T",40,
                            "Lung1_T",15,
                            "Lung3_T",5)

qseaSet |>
    left_join(tumourVAFdata, by = join_by(sample_name)) |>
    getSampleTable() |>
    dplyr::select(sample_name, VAF)

```

We can filter the samples present in the qseaSet by using the `filter` function.
This takes the sample table, filters using the logical criteria given, and
returns a qseaSet with (potentially) fewer samples. All the subcomponents of the
qseaSet are reduced appropriately behind the scenes.

```{r chunk-6}
qseaSet |>
    filter(tumour == "Tumour")
```

There are also functions for sorting the samples of the qseaSet: `sort` and
`arrange` are implemented:

```{r chunk-7}
qseaSet |>
    arrange(tumour)
```

We can also alter columns of the sample table via the `select` function from
`dplyr`. This may be used in both positive and negative selection modes. Note
that select will always retain the columns `sample_name` and `group` in the
output, as they are required for a valid qseaSet.

```{r chunk-8}
qseaSet |>
    dplyr::select(tumour) |>
    getSampleTable()
```

```{r chunk-9}
qseaSet |>
    dplyr::select(-tumour) |>
    getSampleTable()
```

Note that, as usual in R, none of these functions change the original object;
the result needs to be assigned to a new variable (or the same variable).

## Counts Matrix
The counts matrix holds how many fragments of DNA (which may be a single read or
a paired read) were assigned to each window when the qsea object was built. Note
that qsea uses the *centre* of the fragment to assign the window. This counts
matrix should not generally be used directly (we will cover how to make tables
of data later), but you can access it via:

```{r counts-matrix}
qseaSet |>
    getCounts() |>
    head()
```

Note that this object does not have any information as to where each window is
located. That information is all kept inside the regions slot (see below). The
separate vignette `vignette("data-and-qc")` details how to access tables of
data.

## Windows
The regions object is a GenomicRanges (GRanges) object, representing the genomic
windows that the qseaSet covers. This is similar to a data frame, but with
additional information attached to it. In this, the chromosomes are represented
by the `seqnames` column, and the start-stop positions are part of the `ranges`
column. We would highly recommend the use of the `plyranges` package to work
with GRanges objects, see `vignette("an-introduction", package = "plyranges")`.

For DNA methylation this also contains a column called `CpG_density`, which is a
normalised value representing how many CG dinucleotides were present in the
reference genome for that window. For methylation enrichment sequencing, regions
with several methylated CG dinucleotides are more likely to be successfully
pulled down than isolated methylated CGs, so there is a bias towards CG-dense
regions. This bias is stronger when using an MBD protein for the enrichment
step, rather than immunoprecipitation.

```{r windows}
qseaSet |>
    getWindows()
```

We can filter the regions of the qseaSet by using the `filterWindows` function.
Again, this uses `dplyr::filter` to do the logical subsetting. This function can
also use the `seqnames` (i.e. the chromosome), `start` or `end` columns to
filter on.

```{r chunk-10}
qseaSet |>
    filterWindows(CpG_density >= 10)

qseaSet |>
    filterWindows(start >= 26000000)
```

We can also use an already existing GRanges object (or a data frame that has the
columns `seqnames`, `start` and `end`) to perform the intersection with by using
`filterByOverlaps` (this calls the `plyranges` function `filter_by_overlaps`):

```{r chunk-11}
regionOfInterest <- data.frame(seqnames = 7, start = 25000000, end = 26000000)

qseaSet |>
    filterByOverlaps(regionOfInterest)
```

There also exists a `filterByNonOverlaps` function, which will keep the windows
that do not overlap the given regions (using the plyranges function
`filter_by_non_overlaps`):

```{r chunk-12}
qseaSet |>
    filterByNonOverlaps(regionOfInterest)
```

## Library Information
The qseaSet also keeps information about samples in a slot called libraries.
This includes, for instance, the number of fragments in each sample. When we
generate qseaSets through the `makeQset` function (see
`vignette("generation")`), this also includes a set of enrichment-based metrics.
We can use the `addLibraryInformation` function to move these onto the
sampleTable, in order to be able to see or filter on them.

```{r library-information}
qseaSet |>
    addLibraryInformation() |>
    getSampleTable() |>
    as_tibble() |>
    dplyr::select(sample_name, total_fragments, valid_fragments, relH)
```

These columns are as follows:

* *total_fragments*: The total number of fragments remaining after the read
    filtering steps described in `vignette("generation")`
* *valid_fragments*: The total number of fragments remaining after filtering out
    the unreliable regions. A key metric.
* *library_factor*: A qsea-calculated value related to the distribution of the
    reads over the genome
* *offset*: A qsea-calculated value related to the enrichment of the data, with
    how many of the reads are in regions with low CpG density
* *fragment_median*: The median of the fragment size distribution
* *fragment_mean*: The mean of the fragment size distribution
* *fragment_sd*: The standard deviation of the fragment size distribution
* *relH*: The relative enrichment of the pieces of genome covered by the
    fragments, based on (#CGs in fragments)/(#CGs in genome)
* *GoGe*: Another enrichment measurement, based on (normalised) #CGs/(#Cs * #Gs)
    inside the fragments.
* *fragments_without_pattern*: Number of fragments without a CG in the reference
    genome
* *prop_without_pattern*: Proportion of fragments without a CG in the reference
    genome

When copy number variation information is added through separate WGS, there will
also be equivalent columns starting with `input_`, corresponding to the
non-enriched (known as "Input") sample used for estimation of copy number
variation.

As there are so many columns here, we don't add them to the sampleTable by
default, for conciseness, hence the need to use the `addLibraryInformation`
function.

## Copy Number Variation
The qseaSet object can also contain information on the copy number variation
(CNV) that has been calculated for each sample. This is where the number of
copies of each chromosome, or parts of a chromosome, varies, for instance
during carcinogenesis. For known variation of entire chromosomes, e.g. sex
chromosomes, the `zygosity` option can be set, as detailed in the `qsea`
documentation.

This copy number variation is often characterised via (shallow) sequencing of
the DNA without performing the enrichment step, and is stored as a GRanges
object, accessed via the function `getCNV`. `mesa` contains routines to add this
from bam files in the generation step using the package `HMMcopy`, or it can be
added manually if calculated elsewhere (see `vignette("generation")`). The
`qsea` package does include the ability to use off-target reads to calculate
copy number variation, but we find that this gives more noise than parallel
shallow WGS on the same samples, at least in our MBD-seq.

We can access the copy number variation through the `qsea` function `getCNV`:

```{r copy-number-variation}
qseaSet |>
    getCNV() |> 
    head()
```

Here you can see that many of the values are 0 (or close to it), which
represents no significant copy number variation from neutral.

The copy number variation can be plotted using the function `plotCNVheatmap`
(note that this is different to the qsea function `plotCNV`, which uses bright
shades of red and blue):

```{r chunk-13}
qseaSet |>
    plotCNVheatmap()
```

## Parameters, Enrichment, Zygosity
There are also three other slots of the qseaSet, where it stores internal
information: `parameters`, a list of the parameters that were used to generate
the qseaSet; `enrichment`, the fitted enrichment parameters used to calculate
beta values; and `zygosity`, the expected number of copies of each chromosome
in each sample.

These should not typically need to be modified or accessed.

# Parallelisation

Finally, we will mention how parallelisation within the `mesa` functions is set
up. The functions inside the package use the `BiocParallel` package, and will
automatically check if parallelisation has been turned on. This is controlled
via the global option value of `mesa_parallel`:

```{r parallelisation}
options("mesa_parallel")
```

To set this, we use the `setMesaParallel` function:

```{r chunk-14}
setMesaParallel(useParallel = TRUE, nCores = 2)
```

The use of `nCores` here will set up a `MulticoreParam` architecture with that
many cores (a reasonable default for Linux or macOS). Specify a different
architecture with `BiocParallel::register` if you wish to vary this; see
the BiocParallel vignette for how to specify one,
`vignette("Introduction_To_BiocParallel", package = "BiocParallel")`. The most
recently registered instance will be used.

Setting this parameter globally via `options` means that we don't have to
explicitly pass `parallel = TRUE` to functions in order for them to
parallelise across the available cores.

## Session info
```{r session-info, echo=FALSE}
sessionInfo()
```