---
title: "Using TENET"
output:
    BiocStyle::html_document:
        toc: true
        toc_depth: 2
package: TENET
author: Rhie Lab at the University of Southern California
date: "`r Sys.Date()`"
abstract: >
    TENET identifies key transcription factors (TFs) and regulatory elements
    (REs) linked to a specific cell type by finding significantly correlated
    differences in gene expression and RE methylation between case and control
    input datasets, and identifying the top genes by number of significant RE
    DNA methylation site links. It also includes many additional tools to aid
    in visualization and analysis of the results, including plots displaying
    and comparing methylation and expression data and RE DNA methylation site
    link counts, survival analysis, TF motif searching in the vicinity of
    linked RE DNA methylation sites, custom TAD and peak overlap analysis, and
    UCSC Genome Browser track file generation. A utility function is also
    provided to download methylation, expression, and patient survival data
    from The Cancer Genome Atlas (TCGA) for use in TENET or other analyses.
vignette: >
    %\VignetteIndexEntry{Using TENET}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteDepends{MotifDb}
    \usepackage[utf8]{inputenc}
---

\RaggedRight

```{r echo = FALSE, message = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
options(tibble.print_min = 4, tibble.print_max = 4)
```

# Introduction

There is a lack of publicly available bioinformatic tools to identify
transcription factors (TFs) that regulate cell type-specific regulatory
elements (REs). To address this, we developed the Tracing regulatory Element
Networks using Epigenetic Traits (TENET) R package. TENET uses histone mark and
open chromatin datasets, along with matched DNA methylation and gene expression
data, to identify dysregulated REs and the TFs bound to them in a particular
cell or tissue type in comparison with another. To assist in identifying TFs
and REs linked to a particular cell type, we collected hundreds of epigenomic
datasets from a variety of cell lines, primary cells, and tissues and developed
methods to interrogate findings using motif databases, clinical information,
and other genomic datasets from 10 cancer types. Additionally, many downstream
analysis functions have been included to aid users in analyzing and visualizing
results generated by the TENET workflow.

This vignette provides basic instructions to run the workflow of TENET to
identify key TFs and linked RE DNA methylation sites. We will cover how to
install the package, an overview of the necessary input data to use functions
included in the TENET package, and the use of the TENET step 1-7 functions and
the `easyTENET` and `TCGADownloader` utility functions.

# Acquiring and installing TENET and associated packages

To use TENET, users will need to install the base package as well as its
associated example experiment data package,
[TENET.ExperimentHub](https://github.com/rhielab/TENET.ExperimentHub).
**Note:** TENET.ExperimentHub will install automatically when TENET is
installed from Bioconductor.

TENET also uses annotation datasets hosted in the Bioconductor AnnotationHub
database. These datasets will be automatically loaded from AnnotationHub when
necessary. They are also available separately via the
[TENET.AnnotationHub](https://github.com/rhielab/TENET.AnnotationHub) package.
It is not necessary to install the TENET.AnnotationHub package to use TENET's
functions.

R 4.5 or a newer version is required.

On Ubuntu 22.04, successful installation required several additional packages.
They can be installed by running the following command in a terminal:

`sudo apt-get install r-base-dev libcurl4-openssl-dev libfreetype6-dev libfribidi-dev libfontconfig1-dev libharfbuzz-dev libtiff5-dev libxml2-dev`

No dependencies other than R are required on macOS or Windows.

Two versions of this package are available.

To install the stable version from Bioconductor, start R and run:

```{r eval = FALSE}
## Install BiocManager, which is required to install packages from Bioconductor
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

BiocManager::install(version = "devel")
BiocManager::install("TENET")
```
The development version containing the most recent updates is available from
our GitHub repository (<https://github.com/rhielab/TENET>).

To install the development version from GitHub, start R and run:

```{r eval = FALSE}
## Install prerequisite packages to install the development version from GitHub
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
if (!requireNamespace("remotes", quietly = TRUE)) {
    install.packages("remotes")
}

BiocManager::install(version = "devel")
BiocManager::install("rhielab/TENET.ExperimentHub")
BiocManager::install("rhielab/TENET")
```

# Loading TENET

To load the TENET package, start R and run:

```{r message = FALSE}
library(TENET)
```

To load the TENET.ExperimentHub package, start R and run:

```{r message = FALSE}
library(TENET.ExperimentHub)
```

To load the TENET.AnnotationHub package if it has been installed, start R and
run:

```{r message = FALSE}
library(TENET.AnnotationHub)
```

# Running TENET without internet access

Some TENET features and examples download datasets from the internet if they
have not already been cached. You must run `TENETCacheAllData()` once while
connected to the internet before using these TENET features or examples without
internet access (for example, on HPC cluster nodes). See the documentation for
`TENETCacheAllData` for more information.

# Input data

TENET primarily makes use of a MultiAssayExperiment object containing the
following data:

## Expression SummarizedExperiment object

A SummarizedExperiment object named "expression" containing gene expression
data for genes in the human genome. Although gene expression values can be
given in a variety of forms, TENET has been primarily tested using
log2-transformed, upper-quartile normalized fragments per kilobase of
transcript per million mapped reads (FPKM-UQ) values. Gene expression values
for each gene should be given in the rows, while expression values from each
sample should be included in the columns of the assay object in the
"expression" SummarizedExperiment object. Additionally, IDs for each gene
should be included in the rownames of this object and sample IDs should be
included in the colnames.

The samples within the "expression" SummarizedExperiment object should be
matched with those in the "methylation" SummarizedExperiment object.

```{r message = FALSE}
## Load the MultiAssayExperiment package. This is not strictly necessary, but
## allows the user to avoid typing MultiAssayExperiment:: before its functions.
library(MultiAssayExperiment)
```

```{r}
## Load in the example MultiAssayExperiment dataset from the TENET.ExperimentHub
## package
exampleTENETMultiAssayExperiment <-
    TENET.ExperimentHub::exampleTENETMultiAssayExperiment()

## Examine the SummarizedExperiments that should be contained in a
## MultiAssayExperiment object appropriate for use in TENET analyses, including
## the "expression" object
experiments(exampleTENETMultiAssayExperiment)

class(
    experiments(
        exampleTENETMultiAssayExperiment
    )[["expression"]]
)

## Examine data for the first 6 genes and samples in the assay object of the
## "expression" SummarizedExperiment object
assays(
    experiments(
        exampleTENETMultiAssayExperiment
    )[["expression"]]
)[[1]][
    seq_len(6), seq_len(6)
]
```

Additionally, genomic coordinates for genes can be included in a GRanges object
as the rowRanges of the "expression" SummarizedExperiment object. To properly
use TENET functions, this GRanges object should include gene IDs as names,
which match with the IDs used as the rownames of the expression assay object.
Additionally, it should at least include the chromosome, 1-indexed coordinates,
and strand of each gene, and a metadata column named "geneName" which maps the
IDs for genes to the gene names (as TENET generally assumes the user has used
Ensembl IDs for genes, this allows TENET to map these back to the gene names
for data summary and plots). Additional columns can be included, but are not
used by TENET.

If this rowRanges dataset is not included in the "expression"
SummarizedExperiment object, TENET can still be used, but the user will have to
specify a `geneAnnotationDataset` in subsequent functions from which to pull
information for the genes included in the "expression" SummarizedExperiment
object. As an exception, the rowRanges dataset is required if the user is using
the `easyTENET` function to run the step 1 through step 6 functions with
default arguments.

The rowRanges dataset should contain gene and transcript information, and must
be supplied as a GRanges object (such as one imported by `rtracklayer::import`)
or a path to a GFF3 or GTF file. **Note:** TENET has only been tested with
GENCODE and Ensembl gene annotation datasets. If you are using another dataset,
please ensure that it uses the values "gene" and "transcript" for feature
types, which must be stored in a column named "type". In GFF3 files, feature
types may alternatively be stored in the first colon-separated field of the
"ID" column, the second field of which must be the ID itself. Types stored
there will only be used if the "type" column does not contain the required
types. Gene names must be stored in a column named "geneName" or "Name". GTF
files must contain a "geneId" column, and GFF3 files must contain an "ID"
column. An annotation dataset specified as a GRanges object will be assumed to
be derived from a GFF3 file if it contains an "ID" column, and from a GTF file
otherwise. Ensembl GTF files older than release 75 are not supported.

It is not necessary to include a colData object in the "expression"
SummarizedExperiment object, though a colData object in the outer
MultiAssayExperiment object is required.

```{r}
## Examine the rowRanges of the "expression" SummarizedExperiment object for the
## first genes.
## Note: Additional columns are included here, but only the chromosome
## (seqnames), coordinates (rowRanges), strand, and geneName columns are used.
head(
    rowRanges(
        experiments(
            exampleTENETMultiAssayExperiment
        )[["expression"]]
    )
)

## The names of this object are the gene IDs
head(
    names(
        rowRanges(
            experiments(
                exampleTENETMultiAssayExperiment
            )[["expression"]]
        )
    )
)
```

## Methylation SummarizedExperiment object

A SummarizedExperiment object named "methylation" containing DNA methylation
data for methylation sites. Methylation values should be given in the form of
beta ($\beta$) values ranging from 0 (low methylation) to 1 (high methylation)
with values for each RE DNA methylation site in the rows and data from each
individual sample in the columns of the assay object in the "expression"
SummarizedExperiment object. An ID for each RE DNA methylation site (often the
ID of the corresponding probe in a methylation array) should be included in the
rownames of this assay object. As stated above, the samples within the
"methylation" SummarizedExperiment object should be matched with those in the
"expression" SummarizedExperiment object.

```{r}
## Again, examine the SummarizedExperiments included in the MultiAssayExperiment
## object, focusing on the "methylation" object here
experiments(exampleTENETMultiAssayExperiment)

class(
    experiments(
        exampleTENETMultiAssayExperiment
    )[["methylation"]]
)

## Examine data for the first six RE DNA methylation sites and samples in the
## assay object of the "methylation" SummarizedExperiment object
assays(
    experiments(
        exampleTENETMultiAssayExperiment
    )[["methylation"]]
)[[1]][
    seq_len(6), seq_len(6)
]
```

Like the "expression" object, the "methylation" SummarizedExperiment object can
include a GRanges object in its rowRanges with information on the genomic
coordinates for each RE DNA methylation site included in the assay of the
"methylation" SummarizedExperiment. As previously mentioned, this rowRanges
dataset is required if the user is using the `easyTENET` function.

This GRanges object should include the methylation site IDs as names, and these
should match with the IDs used as the rownames of the methylation assay object.
Otherwise, this GRanges object only needs to include the chromosome and
1-indexed coordinates of each RE DNA methylation site. In contrast with the
"expression" SummarizedExperiment object, strand information and names are not
used. Additional columns can be included, but are not used by TENET.

If this rowRanges dataset is not included in the "methylation"
SummarizedExperiment object, TENET can still be used, but the user will have to
specify a `DNAMethylationArray` in subsequent functions from which to pull
information for the RE DNA methylation sites included in the "methylation"
SummarizedExperiment object. If this is the case, it is assumed the user has
provided methylation beta values from probes in one of the Illumina methylation
arrays supported by the sesameData package (see
`?sesameData::sesameData_getManifestGRanges`).

It is not necessary to include a colData object in the "methylation"
SummarizedExperiment object, though a colData object in the outer
MultiAssayExperiment object is required.

```{r}
## Examine the rowRanges of the "methylation" SummarizedExperiment object for
## the first six RE DNA methylation sites.
## Note: Additional columns are included here, but only the chromosome
## (seqnames) and coordinates (ranges) are used.
head(
    rowRanges(
        experiments(
            exampleTENETMultiAssayExperiment
        )[["methylation"]]
    )[, seq_len(6)]
)

## The names of this object are the RE DNA methylation site IDs (usually probe
## IDs)
head(
    names(
        rowRanges(
            experiments(
                exampleTENETMultiAssayExperiment
            )[["methylation"]]
        )
    )
)
```

## MultiAssayExperiment colData object

The colData object should contain information for each of the patients from
which samples in the "expression" and "methylation" SummarizedExperiment objects
have been derived. The rownames of this object should be the patient IDs,
which can be matched to samples using the MultiAssayExperiment object's
sampleMap discussed subsequently. Columns of data can be included in this
dataset to use some downstream step 7 TENET functions but are not essential for
running most TENET functions. These include "vital_status" and "time" columns,
with information on each sample's survival status and survival time, used in the
`step7TopGenesSurvival` function, as well as a "purity" column and columns with
gene copy number ("...\_CNV") and somatic mutation ("...\_SM") status used in
the `step7ExpressionVsDNAMethylationScatterplots` function. See documentation
for these functions for more information on how these data should be formatted.
Additional columns of information can be included, but are not used by TENET.

```{r}
## Examine the number of rows in the colData of the MultiAssayExperiment object
## compared to the number of samples (columns) in the "expression" and
## "methylation" summarized experiment objects. The number of patient entries
## does not need to match the number of samples included in the "expression" or
## "methylation" objects, as a single "Control" and "Case" sample could be
## derived from the same patient (though ideally, no more than one of each)
nrow(colData(exampleTENETMultiAssayExperiment))

experiments(exampleTENETMultiAssayExperiment)

## Examine some of the rownames, which should contain a unique identifier for
## each patient. These will be used in the MultiAssayExperiment's sampleMap
## object to match them to the samples included in the "expression" and
## "methylation" objects
rownames(colData(exampleTENETMultiAssayExperiment))[seq_len(6)]
```

## MultiAssayExperiment sampleMap object

The final essential component which should be included in the
MultiAssayExperiment object for use in TENET analyses is the sampleMap object.
This object is used to match the samples in the "expression" and "methylation"
data objects to each other and match each of the samples in these datasets to
the data contained in the MultiAssayExperiment colData, as well as to note
which samples are "Control" samples, and which are "Case" samples. This object
should have the following format:

```{r}
## The sampleMap object should contain a row for each of the samples included
## in both the "expression" and "methylation" objects
nrow(sampleMap(exampleTENETMultiAssayExperiment))

experiments(exampleTENETMultiAssayExperiment)

## 4 columns of data should be included in the sampleMap, "assay", "primary",
## "colname", and "sampleType"
colnames(sampleMap(exampleTENETMultiAssayExperiment))

## The "assay" column should be a factor which lists which data object each
## sample is from ("expression", or "methylation)
sampleMap(exampleTENETMultiAssayExperiment)$assay[seq_len(6)]

levels(sampleMap(exampleTENETMultiAssayExperiment)$assay)

table(sampleMap(exampleTENETMultiAssayExperiment)$assay)

## The "primary" column should note which of the patient IDs from the
## MultiAssayExperiment's colData object each sample maps to.
sampleMap(exampleTENETMultiAssayExperiment)$primary[seq_len(6)]

table(
    sampleMap(exampleTENETMultiAssayExperiment)$primary %in%
        rownames(colData(exampleTENETMultiAssayExperiment))
)

## The "colname" column should include the sample names of each of the samples
## as they are listed in the colnames of either the "expression" or
## "methylation" SummarizedExperiments' assay objects
sampleMap(exampleTENETMultiAssayExperiment)$colname[seq_len(6)]

table(
    sampleMap(exampleTENETMultiAssayExperiment)$colname %in% c(
        colnames(
            assays(
                experiments(
                    exampleTENETMultiAssayExperiment
                )[["expression"]]
            )[[1]]
        ),
        colnames(
            assays(
                experiments(
                    exampleTENETMultiAssayExperiment
                )[["methylation"]]
            )[[1]]
        )
    )
)

## Finally, the "sampleType" column should list whether each sample in the
## "expression" or "methylation" SummarizedExperiment objects is a "Case" or
## "Control" sample for the purposes of TENET analyses.
sampleMap(exampleTENETMultiAssayExperiment)$sampleType[seq_len(6)]

table(sampleMap(exampleTENETMultiAssayExperiment)$sampleType)
```

# Overview of main TENET functions

* `easyTENET`: Run the step 1 through step 6 functions with default arguments
* `step1MakeExternalDatasets`: Create a GRanges object representing putative regulatory element regions, based on the data sources selected for inclusion, to be used in later TENET steps
* `step2GetDifferentiallyMethylatedSites`: Identify differentially methylated RE DNA methylation sites
* `step3GetAnalysisZScores`: Calculate Z-scores comparing the mean expression of each gene in the case samples that are hyper- or hypomethylated for each RE DNA methylation site chosen in step 2
* `step4SelectMostSignificantLinksPerDNAMethylationSite`: Select the most significant RE DNA methylation site-gene links to each RE DNA methylation site
* `step5OptimizeLinks`: Find final RE DNA methylation site-gene links using various optimization metrics
* `step6DNAMethylationSitesPerGeneTabulation`: Tabulate the total number of RE DNA methylation sites linked to each of the genes
* `TCGADownloader`: Download TCGA gene expression, DNA methylation, and clinical datasets and compile them into a MultiAssayExperiment object
* `TENETCacheAllData`: Cache all online datasets required by TENET examples and optional features

## `easyTENET`: Run the step 1 through step 6 functions with default arguments

This function runs the main six TENET functions, `step1MakeExternalDatasets`,
`step2GetDifferentiallyMethylatedSites`, `step3GetAnalysisZScores`,
`step4SelectMostSignificantLinksPerDNAMethylationSite`, `step5OptimizeLinks`,
and `step6DNAMethylationSitesPerGeneTabulation`, in sequence on the specified
TENETMultiAssayExperiment object. The operation and arguments of this
function reflect those of its component functions. For additional details on
them, refer to their respective sections in this vignette.

`easyTENET` includes only the most important arguments to its component
functions, thus maintaining core TENET functionality while simplifying its use
and allowing users to run all six key TENET steps with one function call. The
majority of its arguments are those from `step1MakeExternalDatasets`, which
select the types of regulatory elements the user wishes to investigate, as well
as those from `step2GetDifferentiallyMethylatedSites`. If an argument to one of
the component functions is not specified, it will be set to its default value
for that function.

As stated earlier, using `easyTENET` requires the user to include rowRanges
objects in the input MultiAssayExperiment's "expression" and "methylation"
SummarizedExperiments containing the locations of transcription start sites
of genes and DNA methylation sites, respectively.

**Note:** Since this function runs the `step3GetAnalysisZScores` function,
it may take a while to run. It is highly recommended to use multiple cores if
possible, especially when using large datasets.

```{r eval = FALSE}
## This example first creates a dataset of putative enhancer regulatory elements
## from consensus datasets and breast invasive carcinoma-relevant sources
## collected in the TENET.AnnotationHub package, then runs the step 2 through
## step 6 functions analyzing RE DNA methylation sites in potential
## enhancer elements located over 1500 bp from transcription start sites
## listed for genes and transcripts in the GENCODE v36 human genome
## annotations (already contained in the exampleTENETMultiAssayExperiment
## object loaded earlier), using a minimum case sample count of 5 and one
## CPU core to perform the analysis.
exampleObject <- easyTENET(
    TENETMultiAssayExperiment = exampleTENETMultiAssayExperiment,
    extHM = NA,
    extNDR = NA,
    publicEnhancer = TRUE,
    publicNDR = TRUE,
    cancerType = "BRCA",
    ENCODEdELS = TRUE,
    minCaseCount = 5
)

## The exampleObject should have data from the step 2 through step 6 functions,
## as the dataset created by the step1MakeExternalDatasets function is used by
## and saved in the output of the step2GetDifferentiallyMethylatedSites
## function.

## See the types of data that were saved by the step 2 function
names(metadata(exampleObject)$step2GetDifferentiallyMethylatedSites)

## See the GRanges object created by the step 1 function
metadata(
    exampleObject
)$step2GetDifferentiallyMethylatedSites$regulatoryElementGRanges

## See the types of data that were saved by the step 3 function
names(metadata(exampleObject)$step3GetAnalysisZScores)

## See the types of data that were saved by the step 4 function
names(
    metadata(exampleObject)$step4SelectMostSignificantLinksPerDNAMethylationSite
)

## See the types of data that were saved by the step 5 function
names(metadata(exampleObject)$step5OptimizeLinks)

## See the types of data that were saved by the step 6 function
names(
    metadata(
        exampleObject
    )$step6DNAMethylationSitesPerGeneTabulation
)
```

## `step1MakeExternalDatasets`: Create a GRanges object representing putative regulatory element regions, based on the data sources selected for inclusion, to be used in later TENET steps

This function creates a GRanges object containing regions representing putative
regulatory elements, either enhancers or promoters, of interest to the user
based on the presence of specific histone marks and open
chromatin/nucleosome-depleted regions. This function can take input from
user-specified bed-like files (see
<https://genome.ucsc.edu/FAQ/FAQformat.html#format1>) containing regions with
histone modification (via the `extHM` argument) and/or open
chromatin/nucleosome-depleted regions (via the `extNDR` argument), as well as
preprocessed enhancer, promoter, and open chromatin datasets from many
cell/tissue types included in the TENET.ExperimentHub package. The resulting
GRanges object will be returned. GRanges objects created by this function can be
used by the `step2GetDifferentiallyMethylatedSites` function or other downstream
functions.

Regulatory element regions of interest identified by this function represent
those with overlapping histone mark peaks as well as open chromatin regions,
combined with any regions identified in the selected ENCODE SCREEN datasets (as
these regions already represent the overlap of regions with relevant histone
marks as well as with open chromatin).

```{r}
## Create an example GRanges object representing putative enhancer regions for
## BRCA using all available enhancer-relevant BRCA datasets present in the
## TENET.ExperimentHub package. These datasets include consensus enhancer
## histone mark and open chromatin datasets from a wide variety of tissue and
## cell types, enhancer histone mark and open chromatin datasets from a
## variety of BRCA-relevant samples from the Cistrome database and TCGA, as well
## as identified distal enhancer regions from the ENCODE SCREEN project.
step1Output <- step1MakeExternalDatasets(
    consensusEnhancer = TRUE,
    consensusNDR = TRUE,
    publicEnhancer = TRUE,
    publicNDR = TRUE,
    cancerType = "BRCA",
    ENCODEdELS = TRUE
)

## View the first regions in the created GRanges object
head(step1Output)
```

## `step2GetDifferentiallyMethylatedSites`: Identify differentially methylated RE DNA methylation sites

This function identifies DNA methylation sites that mark putative regulatory
elements (REs), including enhancer and promoter regions. These are sites
that lie within regions with specific histone modifications and open
chromatin regions, from a user-supplied GRanges object, such as one created
by the `step1MakeExternalDatasets` function, and which are located at a
user-specified distance relative to the transcription start sites (TSS)
listed in either the rowRanges of the elementMetadata of the "expression"
SummarizedExperiment in the TENETMultiAssayExperiment object, or the
selected `geneAnnotationDataset` (which will be filtered to only genes and
transcripts). After identifying DNA methylation sites representing the
specified REs, the function classifies the RE DNA methylation sites as
methylated, unmethylated, hypermethylated, or hypomethylated based on their
differential methylation between the control and case samples supplied by
the user, defined by cutoff values which are either automatically based
on the mean methylation densities of the identified RE DNA methylation
sites, or manually set by the user. **Note:** Using the algorithm to set
cutoffs is recommended for use with DNA methylation array data, and may not
work for whole-genome DNA methylation data.

To run step 2, the user will need to provide a MultiAssayExperiment object
constructed in the manner described previously, as well as a GRanges object,
such as one created by the `step1MakeExternalDatasets` function. Additionally,
the minimum number of case samples that must exhibit a difference in
methylation for a given RE DNA methylation site to be considered hyper- or
hypomethylated will need to be set by the user.

The output of the `step2GetDifferentiallyMethylatedSites` function, as well as
later TENET functions, is saved in the metadata of the returned
MultiAssayExperiment object.

```{r}
## Identify differentially methylated RE DNA methylation sites using the
## step2GetDifferentiallyMethylatedSites function, using the
## exampleTENETMultiAssayExperiment object loaded previously from the
## TENET.ExperimentHub package as a reference, and the GRanges object that was
## just created using the step1MakeExternalDatasets function. At least 5 case
## samples in the dataset are required to show methylation levels above/below
## the hyper/hypomethylation cutoff for a given RE DNA methylation site to be
## potentially considered differentially methylated.
## All transcription start sites (TSS) included in the rowRanges of the
## elementMetadata of the TENETMultiAssayExperiment object will be considered
## when selecting enhancer DNA methylation sites (which must be at least 1500
## bp from these TSS).
exampleObject <- step2GetDifferentiallyMethylatedSites(
    TENETMultiAssayExperiment = exampleTENETMultiAssayExperiment,
    regulatoryElementGRanges = step1Output,
    minCaseCount = 5
)

## See the data that were saved by the step 2 function
names(metadata(exampleObject)$step2GetDifferentiallyMethylatedSites)

## Since cutoffs were set automatically by the step 2 function in this case,
## we can see what they are set to, using the hypomethylation cutoff as an
## example.
metadata(
    exampleObject
)$step2GetDifferentiallyMethylatedSites$hypomethCutoff

## The identities of all identified RE DNA methylation sites, as well as the
## methylated, unmethylated, and most importantly, hyper- and hypomethylated
## RE DNA methylation sites are also recorded in the siteIdentitiesList. To
## demonstrate this, view the first hypomethylated RE DNA methylation sites
## that were identified.
head(
    metadata(
        exampleObject
    )$step2GetDifferentiallyMethylatedSites$siteIdentitiesList$
        hypomethylatedSites
)
```

## `step3GetAnalysisZScores`: Calculate Z-scores comparing the mean expression of each gene in the case samples that are hyper- or hypomethylated for each RE DNA methylation site chosen in step 2

This function calculates Z-scores comparing the mean expression of each gene in
the case samples that are hyper- or hypomethylated for each RE DNA methylation
site chosen in step 2 to the mean expression of the remaining non hyper- or
hypomethylated case samples. By identifying significant Z-scores, initial RE
DNA methylation site-gene links are identified, as case samples with hyper-
or hypomethylation of a particular RE DNA methylation site also display
particularly high or low expression of specific genes.

TENET supports the use of two different formulas for calculating Z-scores in
this step. By setting the zScoreCalculation argument to "oneSample", a
one-sample Z-score calculation will be used (similar to previous versions of
the TENET package), while a two-sample Z-score calculation will be used if the
zScoreCalculation argument is set to "twoSample".

Also, the sparseResults argument has been included in order to reduce the size
of the MultiAssayExperiment object with TENET results. By setting this to TRUE,
only links with significant Z-scores (as determined by the value of the pValue
argument) are saved in the MultiAssayExperiment object returned by this
function. However, setting this to TRUE affects the function of the subsequent
`step4SelectMostSignificantLinksPerDNAMethylationSite` function if the user
wishes to perform multiple testing correction to select the most significant
links per RE DNA methylation site. Therefore, if you want to use multiple
testing correction instead of just selecting the top n most significant links
per RE DNA methylation site in the
`step4SelectMostSignificantLinksPerDNAMethylationSite`, the sparseResults
argument should be set to FALSE so the multiple testing correction is properly
applied for all results, not just the significant ones.

**Note:** This function takes the longest of all TENET functions to run. It is
highly recommended to use multiple cores if possible, especially when using
large datasets.

```{r}
## Identify significant Z-scores and initial RE DNA methylation site-gene links
## using the exampleTENETMultiAssayExperiment with results from the
## step2GetDifferentiallyMethylatedSites function. For this analysis, we will
## use the one-sample Z-score function, consider only TFs, rather than all
## genes, and save only significant Z-scores, to cut down on computational time
## and reduce the size of the returned MultiAssayExperiment object. Two CPU
## cores will be used if they exist (except on Windows where the
## parallel::mclapply function is not supported).
exampleObject <- step3GetAnalysisZScores(
    TENETMultiAssayExperiment = exampleObject,
    pValue = 0.05,
    TFOnly = TRUE,
    zScoreCalculation = "oneSample",
    hypermethAnalysis = TRUE,
    hypomethAnalysis = TRUE,
    includeControl = FALSE,
    sparseResults = TRUE,
    coreCount = min(
        parallel::detectCores(), ifelse(.Platform$OS.type != "windows", 2, 1)
    )
)

## See the data that were saved by the step 3 function. They are subdivided into
## hypermeth and/or hypometh results based on function options.
names(
    metadata(
        exampleObject
    )$step3GetAnalysisZScores
)

## Since the sparseResults argument was set to TRUE, only
## significant Z-scores are saved, and since the TFOnly argument was also set
## to TRUE, only TF genes were analyzed.
## View the significant Z scores for the first TF genes with links to
## hypomethylated RE DNA methylation sites.
head(
    metadata(
        exampleObject
    )$step3GetAnalysisZScores$hypomethResults
)
```

## `step4SelectMostSignificantLinksPerDNAMethylationSite`: Select the most significant RE DNA methylation site-gene links to each RE DNA methylation site

This function takes the calculated Z-scores for the hyper- or hypomethylated G+
RE DNA methylation site-gene links and selects the most significant links to
each RE DNA methylation site, either up to a number specified by the user, or
based on a significant p-value level set by the user after multiple testing
correction is performed on the Z-scores output by the `step3GetAnalysisZScores`
function per RE DNA methylation site in the RE DNA methylation site-gene pairs.
This helps prioritize individual RE DNA methylation site-gene links where
there are many genes linked to a single RE DNA methylation site.

As described previously, if you wish to use multiple testing correction, the
sparseResults argument in the previous `step3GetAnalysisZScores` function
should have been set to FALSE, otherwise it will affect the generated results
(TENET will display a message warning about this if this is the case) as with
sparseResults, only significant results are saved from step 3 and then used in
the multiple testing which affects the values and number of tests accounted for.

This warning will also occur if multiple testing is done using the example
MultiAssayExperiment object, as the results from step 3 in the object were
created with sparseResults set to TRUE. This is just a warning however, and
results will still be generated by the function and can be used in downstream
functions.

```{r}
## Get the 25 (if present) most significant links per RE DNA methylation site
## identified by the step3GetAnalysisZScores function
exampleObject <- step4SelectMostSignificantLinksPerDNAMethylationSite(
    TENETMultiAssayExperiment = exampleObject,
    hypermethGplusAnalysis = TRUE,
    hypomethGplusAnalysis = TRUE,
    linksPerREDNAMethylationSiteMaximum = 25
)

## See the data that were saved by the step 4 function. They are subdivided into
## hypermeth and/or hypometh results based on function options.
names(
    metadata(exampleObject)$step4SelectMostSignificantLinksPerDNAMethylationSite
)

## View the results for the most significant links to the hypomethylated RE
## DNA methylation sites
head(
    metadata(
        exampleObject
    )$step4SelectMostSignificantLinksPerDNAMethylationSite$hypomethGplusResults
)
```

## `step5OptimizeLinks`: Find final RE DNA methylation site-gene links using various optimization metrics

This function takes the most significant hyper- or hypomethylated G+ RE DNA
methylation site-gene links selected in step 4, and selects optimized links
based on the relative expression of the given gene in hyper- or
hypomethylated case samples compared to control samples, using an unpaired
two-sided Wilcoxon rank-sum test to check that the hyper- or hypomethylated
samples for that given RE DNA methylation site-gene link also show
appropriately higher/lower expression of the linked gene in a number of case
samples greater than or equal to the `minCaseCount` number specified in the
`step2GetDifferentiallyMethylatedSites` function that have maximum/minimum
methylation above/below the `hyperStringency`/`hypoStringency` cutoff values
selected.

This identifies the final RE DNA methylation site-gene links by prioritizing
those that meet the above criteria. The output of this function is used in many
of the downstream TENET functions, and helps users examine the individual RE
DNA methylation sites linked to each gene.

```{r}
## Identify final optimized RE DNA methylation site-gene links
exampleObject <- step5OptimizeLinks(
    TENETMultiAssayExperiment = exampleObject,
    hypermethGplusAnalysis = TRUE,
    hypomethGplusAnalysis = TRUE,
    expressionPvalCutoff = 0.05
)

## See the data that were saved by the step 5 function. They are again
## subdivided into hypermeth and/or hypometh results based on function options.
names(metadata(exampleObject)$step5OptimizeLinks)

## Check the results, which include various metrics used to priortize the
## optimized final RE DNA methylation site-gene links.
## This is a subsection of the data frame detailing all the hypomethylated RE
## DNA methylation site-gene links as an example.
head(
    metadata(
        exampleObject
    )$step5OptimizeLinks$hypomethGplusResults
)
```

## `step6DNAMethylationSitesPerGeneTabulation`: Tabulate the total number of RE DNA methylation sites linked to each of the genes

This function takes the final optimized RE DNA methylation site-gene links
identified in step 5 and tabulates the number of these links per gene. This
tabulation is done separately for both of the hyper- or hypomethylated G+
analysis quadrants, as selected by the user.

This aids in prioritizing the top genes for downstream analyses, as the genes
with the most linked RE DNA methylation sites are the most likely to represent
those with widespread genomic impact.

```{r}
## Prioritize the top genes by adding up the number of RE DNA methylation sites
## linked to each of the genes
exampleObject <- step6DNAMethylationSitesPerGeneTabulation(
    TENETMultiAssayExperiment = exampleObject
)

## See the data that were saved by the step 6 function. They are again
## subdivided into hypermeth and/or hypometh results based on function options.
names(
    metadata(
        exampleObject
    )$step6DNAMethylationSitesPerGeneTabulation
)

## Check the results, which include a count of the RE DNA methylation sites per
## gene, and is organized by decreasing RE DNA methylation site count.
## This is a subsection of the data frame detailing the number of hypomethylated
## RE DNA methylation site links to the top TFs.
head(
    metadata(
        exampleObject
    )$step6DNAMethylationSitesPerGeneTabulation$hypomethGplusResults
)
```

## `TCGADownloader`: Download TCGA gene expression, DNA methylation, and clinical datasets and compile them into a MultiAssayExperiment object

This function downloads and compiles TCGA gene expression and DNA
methylation datasets, as well as clinical data primarily intended for use
with the TENET package. This simplifies the TCGAbiolinks download functions,
identifies samples with matching gene expression and DNA methylation data,
and can also remove duplicate tumor samples taken from the same patient
donor. Data are compiled into a MultiAssayExperiment object, which is
returned and optionally saved in an .rda file at the path specified by the
`outputFile` argument.

```{r eval = FALSE}
## Download a TCGA BRCA dataset with log2-normalized
## FPKM-UQ expression values from tumor and adjacent normal tissue samples
## with matching expression and methylation data and keeping only one tumor
## sample from each patient. Additionally, survival data will be combined
## from three clinical datasets downloaded by TCGAbiolinks. Raw data files
## will be saved to the working directory, and the processed dataset will
## be returned as a variable.
TCGADataset <- TCGADownloader(
    rawDataDownloadDirectory = ".",
    TCGAStudyAbbreviation = "BRCA",
    RNASeqWorkflow = "STAR - FPKM-UQ",
    RNASeqLog2Normalization = TRUE,
    matchingExpAndMetSamples = TRUE,
    clinicalSurvivalData = "combined",
    outputFile = NA
)
```

## `TENETCacheAllData`: Cache all online datasets required by TENET examples and optional features

This function locally caches all online TENET and SeSAMe datasets required
by TENET examples and optional features (TENET.ExperimentHub objects used in
examples, TENET.AnnotationHub datasets used in step 1, and SeSAMe datasets
loaded via the `DNAMethylationArray` argument). The main purpose of this
function is to enable the use of TENET in an HPC cluster environment where
compute nodes do not have internet access. In this case, you must run
`TENETCacheAllData()` once while connected to the internet before using
TENET examples or these optional features.

```{r eval = FALSE}
## Cache all online datasets required by optional TENET features.
## This function takes no arguments and returns NULL.
TENETCacheAllData()
```

# Overview of downstream step 7 functions

The step 7 functions aim to perform downstream analyses based on the identified
RE DNA methylation site-gene links, integrating multi-omic datasets such as
Hi-C, copy number variation, somatic mutation, and patient survival information.

* `step7ExpressionVsDNAMethylationScatterplots`: Create scatterplots displaying the expression of the top genes and the methylation levels of each of their linked RE DNA methylation sites, along with copy number variation, somatic mutation, and purity data, if provided by the user
* `step7LinkedDNAMethylationSiteCountHistograms`: Create histograms displaying the number of genes or transcription factors linked to a given number of RE DNA methylation sites
* `step7LinkedDNAMethylationSitesMotifSearching`: Perform motif searching for transcription factor motifs in the vicinity of RE DNA methylation sites linked to the specified transcription factors
* `step7SelectedDNAMethylationSitesCaseVsControlBoxplots`: Generate boxplots comparing the methylation level of the specified RE DNA methylation sites in case and control samples
* `step7StatesForLinks`: Identify which of the case samples harbor each of the identified regulatory element DNA methylation site-gene links
* `step7TopGenesCaseVsControlExpressionBoxplots`: Create boxplots comparing the expression level of the top genes/transcription factors in case and control samples
* `step7TopGenesCircosPlots`: Generate Circos plots displaying the links between top identified genes and each of the RE DNA methylation sites linked to them
* `step7TopGenesDNAMethylationHeatmaps`: Generate heatmaps displaying the methylation level of all RE DNA methylation sites linked to the top genes/transcription factors, along with the expression of those genes in the column headers, in the case samples within the supplied MultiAssayExperiment object
* `step7TopGenesExpressionCorrelationHeatmaps`: Generate mirrored heatmaps displaying the correlation of the expression values of the top genes/TFs
* `step7TopGenesOverlappingLinkedDNAMethylationSitesHeatmaps`: Generate binary heatmaps displaying which of the top genes/transcription factors share links with each of the unique regulatory element DNA methylation sites linked to at least one top gene/TF
* `step7TopGenesSurvival`: Perform Kaplan-Meier and Cox regression analyses to assess the association of top gene expression and linked RE DNA methylation site methylation with patient survival
* `step7TopGenesTADTables`: Create tables using user-supplied topologically associating domain (TAD) information which identify the topologically associating domains containing each RE DNA methylation site linked to the top genes/transcription factors, as well as other genes in the same topologically associating domain as potential downstream targets
* `step7TopGenesUCSCBedFiles`: Create bed-formatted interact files which can be loaded on the UCSC Genome Browser to display links between top genes and transcription factors and their linked RE DNA methylation sites
* `step7TopGenesUserPeakOverlap`: Identify if RE DNA methylation sites linked to top genes/transcription factors are located within a specific distance of specified genomic regions

Here we will demonstrate the usage of some step 7 functions.

## `step7ExpressionVsDNAMethylationScatterplots`: Create scatterplots displaying the expression of the top genes and the methylation levels of each of their linked RE DNA methylation sites, along with copy number variation, somatic mutation, and purity data, if provided by the user

These scatterplots show the relationship between genes and RE DNA methylation
sites, displaying the expression of the genes in the X-axis and the methylation
of the sites in the Y-axis. The sample type (case or control) is also displayed
in these plots. The shape and size of the points on the scatterplots represent
copy number variation (CNV), somatic mutation (SM) status, and purity for the
samples in the scatterplots.

First, we load the example CNV, SM, and purity data from the
`exampleTENETClinicalDataFrame` object.

```{r}
## Load the exampleTENETClinicalDataFrame object from the TENET.ExperimentHub
## package. It contains copy number variation (CNV), somatic mutation (SM),
## and purity data for the top 10 TFs by linked hypomethylated RE
## DNA methylation sites in the exampleTENETMultiAssayExperiment object.
exampleTENETClinicalDataFrame <-
    TENET.ExperimentHub::exampleTENETClinicalDataFrame()
CNVData <- subset(exampleTENETClinicalDataFrame,
    select = grepl("_CNV$", colnames(exampleTENETClinicalDataFrame))
)
SMData <- subset(exampleTENETClinicalDataFrame,
    select = grepl("_SM$", colnames(exampleTENETClinicalDataFrame))
)
purityData <- subset(exampleTENETClinicalDataFrame, select = "purity")
```

The CNV dataset is a numeric data frame with rownames representing sample names
and colnames representing gene IDs folllowed by "\_CNV", with -2 representing a
loss of both copies, -1 a single copy loss, 0 no copy number change, and
positive integer values representing a gain of that many copies (though changes
of +2 or greater are grouped together in the scatterplots).

```{r}
## Show the CNV data for the first 4 TFs
head(CNVData[, 1:4])
```

The SM dataset is a numeric data frame with rownames representing sample names
and colnames representing gene IDs folllowed by "\_SM", with 0 representing no
mutation and 1 representing mutation.

```{r}
## Show the SM data for the first 4 TFs
head(SMData[, 1:4])
```

The purity dataset in this example is a numeric data frame with rownames
representing sample names and the first column representing the purity, with
the values ranging from 0 (low purity) to 1 (high purity). It can also be a
numeric vector with names representing the sample names.

```{r}
## Show the first few rows of the purity data
head(purityData)
```

The options `CNVData`, `SMData`, and `purityData` are not required. If they are
supplied and `simpleOrComplex` is set to "complex", complex scatterplots will
be created displaying this information. Otherwise, simple scatterplots will be
created. At this time, either all or none of `CNVData`, `SMData`, and
`purityData` must be specified.

```{r}
## Create complex scatterplots using the previously acquired data.
## Since we performed analyses only using TFs in the step 3 function, the
## top genes are all TFs, so a message that separate output for TFs will be
## skipped is displayed.
exampleObject <- step7ExpressionVsDNAMethylationScatterplots(
    TENETMultiAssayExperiment = exampleObject,
    hypermethGplusAnalysis = FALSE,
    hypomethGplusAnalysis = TRUE,
    purityData = purityData,
    SMData = SMData,
    CNVData = CNVData,
    simpleOrComplex = "complex",
    topGeneNumber = 10
)

## Results are included in the metadata of the returned MultiAssayExperiment
## object under the step7ExpressionVsDNAMethylationScatterplots list.
## For each analysis type, results are included in sub-lists, each of which
## contains results for topGenes and topTFs, unless the top genes are
## all TFs, in which case the separate top TFs output is skipped.
names(
    metadata(
        exampleObject
    )$step7ExpressionVsDNAMethylationScatterplots$hypomethGplusResults$topGenes
)

## For each gene, scatterplots are generated showing the expression of that
## gene and the methylation of each RE DNA methylation site linked to it for
## the given analysis type.
head(
    names(
        metadata(
            exampleObject
        )$step7ExpressionVsDNAMethylationScatterplots$hypomethGplusResults$
            topGenes$ENSG00000124664
    )
)
```

As an example, we examine the scatterplot with the expression of the TF SPDEF
(ENSG00000124664) and the methylation of its linked hypomethylated RE DNA
methylation site with the ID cg25731211. Gene expression is given in the X-axis
and methylation is given in the Y-axis. Samples are colored based on whether
they are are cases (red) or controls (blue). The shape and size of the points
are determined by each sample's CNV/SM status and purity, respectively, since
complex scatterplots were selected.

```{r fig.width=10, fig.height=7}
metadata(
    exampleObject
)$step7ExpressionVsDNAMethylationScatterplots$hypomethGplusResults$
    topGenes$ENSG00000124664$cg25731211
```

## `step7LinkedDNAMethylationSiteCountHistograms`: Create histograms displaying the number of genes or transcription factors linked to a given number of RE DNA methylation sites

This function generates histograms displaying the frequency of genes by the
number of RE DNA methylation sites linked to them. These are designed to
highlight the top genes/TFs, which likely have a disproportionately large
number of linked RE DNA methylation sites compared to most genes/TFs.

```{r fig.width=10, fig.height=7}
## Run the step7LinkedDNAMethylationSiteCountHistograms function.
## Since we performed analyses only using TFs in the step 3 function, the
## top genes are all TFs, so a message that separate output for
## TFs will be skipped is displayed.
exampleObject <- step7LinkedDNAMethylationSiteCountHistograms(
    TENETMultiAssayExperiment = exampleObject,
    hypomethGplusAnalysis = TRUE,
    hypermethGplusAnalysis = FALSE
)

## Results are included in the metadata of the returned MultiAssayExperiment
## object under the step7LinkedDNAMethylationSiteCountHistograms list.
## For each analysis type, histograms are included in sub-lists, each of which
## contains results for topGenes and topTFs, unless the top genes are all TFs,
## in which case the separate top TFs output is skipped.

## Display the histogram. Note the relatively small number of top TF genes with
## larger numbers of linked RE DNA methylation sites.
metadata(
    exampleObject
)$step7LinkedDNAMethylationSiteCountHistograms$hypomethGplusResults$topGenes
```

## `step7LinkedDNAMethylationSitesMotifSearching`: Perform motif searching for transcription factor motifs in the vicinity of RE DNA methylation sites linked to the specified transcription factors

To run the motif searching function, it is necessary to provide position weight
matrices (PWMs), which represent DNA binding motifs for the TFs of interest in
a named list, with each of the PWMs named after the TFs they represent in the
TENET dataset.

The easiest way to get PWMs is to use the MotifDb package to search for
available PWMs for a given TF. In this example, we search for PWMs which are
available for the FOXA1 and ESR1 TFs.

```{r}
## View the first few available motifs for the FOXM1 TF
head(names(MotifDb::query(MotifDb::MotifDb, "FOXA1")))

## View the first few available motifs for the ESR1 TF
head(names(MotifDb::query(MotifDb::MotifDb, "ESR1")))
```

Next, we create a named list using the human HOCOMOCO v11 core motif PWMs
available for both TFs, which will be used when running the
`step7LinkedDNAMethylationSitesMotifSearching` function.

```{r}
## The human HOCOMOCO v11 core motif is the 3rd listed for FOXA1, and 4th for
## ESR1
TFMotifList <- list(
    "FOXA1" = MotifDb::query(MotifDb::MotifDb, "FOXA1")[[3]],
    "ESR1" = MotifDb::query(MotifDb::MotifDb, "ESR1")[[4]]
)

TFMotifList
```

Finally, we run the `step7LinkedDNAMethylationSitesMotifSearching` function to
search for the selected FOXA1 and ESR1 motifs in the vicinity of identified
hypomethylated RE DNA methylation sites linked to the RE DNA methylation sites
found to be linked to those TFs in the TENET analyses that were performed
earlier. A threshold of 80% is chosen to assess motif accuracy to surrounding
DNA sequences within 100 base pairs of RE DNA methylation sites. Increasing the
threshold value makes the search more strict, and reduces the number of motifs
found, while decreasing this value does the opposite. Also, longer motifs will
tend to have fewer matches than shorter motifs.

**Note:** Motif searching can take some time. It is highly recommended to run
this function with multiple cores if possible.

```{r eval = FALSE}
exampleObject <- step7LinkedDNAMethylationSitesMotifSearching(
    TENETMultiAssayExperiment = exampleObject,
    TFMotifList = TFMotifList,
    matchPWMMinScore = "80%"
)

## For each analysis type and TF, a seqLogo diagram of the chosen PWM and two
## data frames with information on the found motifs in the vicinity of RE
## DNA methylation sites, and total motifs found per RE DNA methylation site
## linked to those TFs, respectively, are saved in the metadata of the returned
## MultiAssayExperiment object under the
## step7LinkedDNAMethylationSitesMotifSearching list
names(
    metadata(
        exampleObject
    )$step7LinkedDNAMethylationSitesMotifSearching$hypomethGplusResults$FOXA1
)

## View the motif occurrences near hypomethylated RE DNA methylation sites
## linked to the FOXA1 TF
head(
    metadata(
        exampleObject
    )$step7LinkedDNAMethylationSitesMotifSearching$hypomethGplusResults$
        FOXA1$DNAMethylationSiteMotifOccurrences
)

## Also see the total number of motifs found in the vicinity of each RE DNA
## methylation site
head(
    metadata(
        exampleObject
    )$step7LinkedDNAMethylationSitesMotifSearching$hypomethGplusResults$
        FOXA1$totalMotifOccurrencesPerREDNAMethylationSite
)
```

## `step7SelectedDNAMethylationSitesCaseVsControlBoxplots`: Generate boxplots comparing the methylation level of the specified RE DNA methylation sites in case and control samples

We begin this example by identifying some RE DNA methylation sites for which to
generate methylation boxplots.

First, we look at the top genes by number of linked hypomethylated RE DNA
methylation sites.

```{r}
head(
    metadata(
        exampleObject
    )$step6DNAMethylationSitesPerGeneTabulation$hypomethGplusResults
)
```

Next, we retrieve hypomethylated RE DNA methylation sites linked to the FOXA1
(ENSG00000129514) TF. They can be acquired from the output of the
`step5OptimizeLinks` function.

```{r}
DNAMethylationSites <- subset(
    metadata(
        exampleObject
    )$step5OptimizeLinks$hypomethGplusResults,
    geneID == "ENSG00000129514"
)$DNAMethylationSiteID
head(DNAMethylationSites)
```

Finally, we generate boxplots for the selected RE DNA methylation sites.

```{r}
exampleObject <- step7SelectedDNAMethylationSitesCaseVsControlBoxplots(
    TENETMultiAssayExperiment = exampleObject,
    DNAMethylationSites = DNAMethylationSites
)

## Each plot is saved under the ID of the RE DNA methylation site and included
## in the metadata of the returned MultiAssayExperiment object under the
## step7SelectedDNAMethylationSitesCaseVsControlBoxplots list
head(names(
    metadata(
        exampleObject
    )$step7SelectedDNAMethylationSitesCaseVsControlBoxplots
))
```

As an example, we examine the boxplot for the RE DNA methylation site with the
ID cg13776095.

```{r fig.width=10, fig.height=7}
## Note: There may be a warning for "rows containing non-finite values" if there
## are any samples lacking methylation data for the RE DNA methylation site.
metadata(
    exampleObject
)$step7SelectedDNAMethylationSitesCaseVsControlBoxplots$cg13776095
```

## `step7StatesForLinks`: Identify which of the case samples harbor each of the identified regulatory element DNA methylation site-gene links

This function identifies which of the samples provided in the dataset likely
harbor each of the RE DNA methylation site-gene links. This is accomplished by
examining if the RE DNA methylation site in a given link is hyper- or
hypomethylated in a given case sample, and if expression of the gene in the
link for that case sample is significantly less or greater than, respectively,
mean expression of the gene in the control samples. This is potentially helpful
by identifying case samples for further analyses.

```{r}
## Calculate potential link status for case samples for hypomethylated RE DNA
## methylation site-gene links
exampleObject <- step7StatesForLinks(
    TENETMultiAssayExperiment = exampleObject,
    hypomethGplusAnalysis = TRUE,
    hypermethGplusAnalysis = FALSE
)

## Results are included in the metadata of the returned MultiAssayExperiment
## object under the step7StatesForLinks list.
## A single data frame is returned for each analysis type, with the combined RE
## DNA methylation site ID and gene ID from each link in the rows, and each case
## sample in the columns.
dim(
    metadata(
        exampleObject
    )$step7StatesForLinks$hypomethGplusResults
)

## Show the results for the first 6 case samples. 1 indicates that a given case
## sample might harbor a given RE DNA methylation site-gene link, and 0
## indicates that it does not. NA values are shown for samples that lack
## methylation data for the site or expression data for the gene.
head(
    metadata(
        exampleObject
    )$step7StatesForLinks$hypomethGplusResults[
        , seq_len(6)
    ]
)
```

## `step7TopGenesCaseVsControlExpressionBoxplots`: Create boxplots comparing the expression level of the top genes/transcription factors in case and control samples

This function generates boxplots comparing the expression of the top genes/TFs
by number of linked RE DNA methylation sites in the case versus control samples.

```{r}
## Run the step7TopGenesCaseVsControlExpressionBoxplots function.
## Since we performed analyses only using TFs in the step 3 function, the
## top genes are all TFs, so a message that separate output for
## TFs will be skipped is displayed.
exampleObject <- step7TopGenesCaseVsControlExpressionBoxplots(
    TENETMultiAssayExperiment = exampleObject,
    hypomethGplusAnalysis = TRUE,
    hypermethGplusAnalysis = FALSE,
    topGeneNumber = 10
)

## Results are included in the metadata of the returned MultiAssayExperiment
## object under the step7TopGenesCaseVsControlExpressionBoxplots list.
## For each analysis type, results are included in sub-lists, each of which
## contains lists with results for topGenes and topTFs, unless the
## top genes are all TFs, in which case the separate top TFs output is skipped.
## Each boxplot is saved under its gene ID.
names(
    metadata(
        exampleObject
    )$step7TopGenesCaseVsControlExpressionBoxplots$
        hypomethGplusResults$topGenes
)
```

As an example, we examine the boxplot for gene ENSG00000107485 (GATA3).

```{r fig.width=10, fig.height=7}
## Note: There may be a warning for "rows containing non-finite values" if there
## are any samples lacking expression data for the given gene.
metadata(
    exampleObject
)$step7TopGenesCaseVsControlExpressionBoxplots$hypomethGplusResults$
    topGenes$ENSG00000107485
```

## `step7TopGenesCircosPlots`: Generate Circos plots displaying the links between top identified genes and each of the RE DNA methylation sites linked to them

This function generates Circos plots for each of the top genes by number of
linked RE DNA methylation sites showing the links between the gene and each of
its RE DNA methylation sites.

```{r fig.width=8, fig.height=8, out.width="80%", out.height="80%"}
## Run the step7TopGenesCircosPlots function
exampleObject <- step7TopGenesCircosPlots(
    TENETMultiAssayExperiment = exampleObject,
    hypermethGplusAnalysis = FALSE,
    hypomethGplusAnalysis = TRUE,
    topGeneNumber = 10
)

## Results are included in the metadata of the returned MultiAssayExperiment
## object under the step7TopGenesCircosPlots list.
## For each analysis type, results are included in sub-lists, each
## of which contains lists with results for topGenes and topTFs, unless the
## top genes are all TFs, in which case the separate top TFs output is skipped.
## Each Circos plot is saved under its gene ID.

## Note: Since we performed analyses only using TFs in the step 3 function,
## the top genes are all TFs, so topTFs will be NA here.
names(
    metadata(
        exampleObject
    )$step7TopGenesCircosPlots$hypomethGplusResults$topGenes
)

## Display an example Circos plot for ENSG00000091831 (ESR1).
## Note: Plots may take some time to display.
metadata(
    exampleObject
)$step7TopGenesCircosPlots$hypomethGplusResults$topGenes$ENSG00000091831
```

## `step7TopGenesExpressionCorrelationHeatmaps`: Generate mirrored heatmaps displaying the correlation of the expression values of the top genes/TFs

This function generates heatmaps displaying the correlation of the expression
of each of the top genes in the case samples. Each of the top genes is
displayed in both the rows and columns, so the heatmaps are mirrored, with
correlation values of each gene to itself displayed in a diagonal line up the
center of the heatmaps. Red values represent positive correlation and blue
values represent negative correlation, with darker colors representing a
stronger correlation. Dendrograms are included to identify genes which are
closely related in expression correlation.

```{r fig.width=10, fig.height=10.5, out.width="80%", out.height="80%"}
## Run the step7TopGenesExpressionCorrelationHeatmaps function
exampleObject <- step7TopGenesExpressionCorrelationHeatmaps(
    TENETMultiAssayExperiment = exampleObject,
    hypermethGplusAnalysis = FALSE,
    hypomethGplusAnalysis = TRUE,
    topGeneNumber = 10
)

## Results are included in the metadata of the returned MultiAssayExperiment
## object under the step7TopGenesExpressionCorrelationHeatmaps list.
## For each analysis type, results are included in sub-lists, each
## of which contains lists with results for the topGenes and topTFs, unless
## the top genes are all TFs, in which case the separate top TFs output is
## skipped. For each of these, the heatmap is generated along with a data frame
## with the correlation values displayed in the heatmap.

## Display the mirrored heatmap.
## Note: Since we performed analyses only using TFs in the step 3 function,
## the top genes are all TFs, so topTFs will be NA here.
metadata(
    exampleObject
)$step7TopGenesExpressionCorrelationHeatmaps$hypomethGplusResults$
    topGenes$heatmap

## Display the data frame with correlation values
head(
    metadata(
        exampleObject
    )$step7TopGenesExpressionCorrelationHeatmaps$hypomethGplusResults$
        topGenes$correlationMatrix
)
```

## `step7TopGenesDNAMethylationHeatmaps`: Generate heatmaps displaying the methylation level of all RE DNA methylation sites linked to the top genes/transcription factors, along with the expression of those genes in the column headers, in the case samples within the supplied MultiAssayExperiment object

This function creates heatmaps displaying the methylation of unique RE DNA
methylation sites linked to the top genes in the main body of the heatmaps, as
well as a smaller heatmap showing expression of the top genes labeling the
columns. Expression/methylation for each case sample is shown per column, while
expression of each of the top genes or methylation of their linked RE DNA
methylation sites is shown in the rows. Warm colors represent relatively higher
expression/methylation levels, while cold colors represent relatively lower
expression/methylation levels. These are determined per gene/RE DNA methylation
site, and are not comparable between genes/RE DNA methylation sites, only
between samples. Column dendrograms are included to identify subsets of the
case samples which display particular expression or methylation patterns in the
top genes and their linked RE DNA methylation sites.

```{r}
## Run the step7TopGenesDNAMethylationHeatmaps function
exampleObject <- step7TopGenesDNAMethylationHeatmaps(
    TENETMultiAssayExperiment = exampleObject,
    hypermethGplusAnalysis = FALSE,
    hypomethGplusAnalysis = TRUE,
    topGeneNumber = 10
)

## Results are included in the metadata of the returned MultiAssayExperiment
## object under the step7TopGenesDNAMethylationHeatmaps list.
## For each analysis type, results are included in sub-lists, each
## of which contains heatmaps for the topGenes and topTFs, unless the
## top genes are all TFs, in which case the separate top TFs output is skipped.

## Note: Since we performed analyses only using TFs in the step 3 function,
## the top genes are all TFs, so topTFs will be NA here.
names(metadata(
    exampleObject
)$step7TopGenesDNAMethylationHeatmaps$hypomethGplusResults)
```

## `step7TopGenesOverlappingLinkedDNAMethylationSitesHeatmaps`: Generate binary heatmaps displaying which of the top genes/transcription factors share links with each of the unique regulatory element DNA methylation sites linked to at least one top gene/TF

These binary heatmaps provide a visual representation of which of the top genes
each RE DNA methylation site is linked to, as RE DNA methylation sites can be
linked to multiple genes. RE DNA methylation sites are displayed in the
columns, while the top genes are displayed in the rows. Black indicates that a
given RE DNA methylation site is linked to that gene, and white indicates it is
not. Dendrograms are included to identify blocks of RE DNA methylation sites
that are linked to similar genes.

```{r fig.width=20, fig.height=14, out.width="80%", out.height="80%"}
## Run the step7TopGenesOverlappingLinkedDNAMethylationSitesHeatmaps function
exampleObject <- step7TopGenesOverlappingLinkedDNAMethylationSitesHeatmaps(
    TENETMultiAssayExperiment = exampleObject,
    hypermethGplusAnalysis = FALSE,
    hypomethGplusAnalysis = TRUE,
    topGeneNumber = 10
)

## Results are included in the metadata of the returned MultiAssayExperiment
## object under the step7TopGenesOverlappingLinkedDNAMethylationSitesHeatmaps
## list. For each analysis type, results are included in sub-lists, each of
## which contains lists with results for the topGenes and topTFs, unless
## the top genes are all top TFs,  in which case the separate top TFs output is
## skipped. For each of these, the heatmap is generated along with a data frame
## with the correlation values displayed in the heatmap.

## Display the binary heatmap.
## Note: Since we performed analyses only using TFs in the step 3 function,
## the top genes are all TFs, so topTFs will be NA here.
metadata(
    exampleObject
)$step7TopGenesOverlappingLinkedDNAMethylationSitesHeatmaps$
    hypomethGplusResults$topGenes$heatmap

## Display a subset of the data frame noting the presence/absence of links.
## 1 indicates a link, while 0 indicates no link.
head(
    metadata(
        exampleObject
    )$step7TopGenesOverlappingLinkedDNAMethylationSitesHeatmaps$
        hypomethGplusResults$topGenes$linkTable[
        , seq_len(6)
    ]
)
```

## `step7TopGenesSurvival`: Perform Kaplan-Meier and Cox regression analyses to assess the association of top gene expression and linked RE DNA methylation site methylation with patient survival

This function uses the survival status and time for each of the samples in the
dataset to perform survival analyses on the expression of the top genes and
methylation of their linked RE DNA methylation sites. For each gene and RE DNA
methylation site, a Kaplan-Meier survival plot can be generated, and statistics
from both Kaplan-Meier and univariate Cox regression analyses are output.

First, we load the example survival status and survival time data from the
`exampleTENETClinicalDataFrame` object.

```{r}
## Load the exampleTENETClinicalDataFrame object from the TENET.ExperimentHub
## package. It contains the vital_status (survival status) and time (survival
## time) data for each sample in the exampleTENETMultiAssayExperiment
exampleTENETClinicalDataFrame <-
    TENET.ExperimentHub::exampleTENETClinicalDataFrame()
vitalStatusData <- subset(
    exampleTENETClinicalDataFrame,
    select = "vital_status"
)
survivalTimeData <- subset(exampleTENETClinicalDataFrame, select = "time")
```

The vital status dataset is a data frame with rownames representing sample
names and the first column representing the vital status. Sample values are
either "alive" or "dead" (case-insensitive) or 1 or 2, indicating that samples
were collected from a patient who was alive/censored or dead/reached the
outcome of interest, respectively.

Similarly, the survival time dataset is a data frame with rownames representing
sample names and the first column representing the survival time of the patient
the sample was derived from.

```{r}
## Show the vital status data
head(vitalStatusData)

## Show the survival time data
head(survivalTimeData)
```

Next, we perform the survival analysis using the vital status and survival time
data.

```{r fig.width=10, fig.height=7}
## Since we performed analyses only using TFs in the step 3 function, the
## top genes are all TFs, so a message that separate output for
## TFs will be skipped is displayed.
exampleObject <- step7TopGenesSurvival(
    TENETMultiAssayExperiment = exampleObject,
    hypermethGplusAnalysis = FALSE,
    hypomethGplusAnalysis = TRUE,
    vitalStatusData = vitalStatusData,
    survivalTimeData = survivalTimeData,
    topGeneNumber = 10,
    generatePlots = TRUE
)

## Results are included in the metadata of the returned MultiAssayExperiment
## object under the step7TopGenesSurvival list.
## For each analysis type, results are included in sub-lists, each
## of which contains lists with results for topGenes and topTFs, unless the
## top genes are all TFs, in which case the separate top TFs output is skipped.
## Each includes two data frames with the survival statistics for Kaplan-Meier
## and Cox regression survival analyses, and if the generatePlots
## argument is TRUE, topGenesSurvivalPlots and topMethylationSitesSurvivalPlots
## lists are included which contain the Kaplan-Meier survival plots for the top
## genes and each of their unique linked RE DNA methylation sites, respectively.
names(
    metadata(
        exampleObject
    )$step7TopGenesSurvival$hypomethGplusResults$topGenes
)

## The topGenesSurvivalStats and topMethylationSitesSurvivalStats variables are
## data frames containing survival statistics.
## Note: A significant amount of data is output, so selected values are shown
## here.
head(
    metadata(
        exampleObject
    )$step7TopGenesSurvival$hypomethGplusResults$
        topGenes$topGenesSurvivalStats[
        , c(1:2, 15, 17, 22:24, 26)
    ]
)

head(
    metadata(
        exampleObject
    )$step7TopGenesSurvival$hypomethGplusResults$
        topGenes$topMethylationSitesSurvivalStats[
        , c(1, 24, 26, 31:33, 35)
    ]
)

## Show the names of the gene survival plots
names(
    metadata(
        exampleObject
    )$step7TopGenesSurvival$hypomethGplusResults$
        topGenes$topGenesSurvivalPlots
)

## Plot the Kaplan-Meier survival plot for GATA3 (ENSG00000107485) as an example
metadata(
    exampleObject
)$step7TopGenesSurvival$hypomethGplusResults$
    topGenes$topGenesSurvivalPlots$ENSG00000107485
```

## `step7TopGenesTADTables`: Create tables using user-supplied topologically associating domain (TAD) information which identify the topologically associating domains containing each RE DNA methylation site linked to the top genes/transcription factors, as well as other genes in the same topologically associating domain as potential downstream targets

This function requires the user to supply either a path to a directory
containing bed-like files with TAD regions of interest, or a single TAD object
given as a GRanges, data frame, or matrix object, as seen in the example below.
To illustrate the use of this function, we will use an example TAD dataset from
the TENET.ExperimentHub package.

```{r}
## Load the example TAD dataset from the TENET.ExperimentHub package
exampleTADRegions <- TENET.ExperimentHub::exampleTENETTADRegions()

## TAD files for this function must include the chromosome of each TAD region
## in the first column, and the start and end positions of each in the second
## and third columns respectively. Additional columns can be included but
## are not considered in this function.
class(exampleTADRegions)
head(exampleTADRegions)
```

The unique RE DNA methylation sites linked to the top genes, as selected by the
user, will be overlapped with the TAD files, and genes within the same TAD of
each RE DNA methylation site will be recorded (as possible downstream target
genes for the regulatory elements represented by those RE DNA methylation sites,
for further analysis purposes).

```{r}
## Use the example TAD object to perform TAD overlapping.
## Since we performed analyses only using TFs in the step 3 function, the
## top genes are all TFs, so a message that separate output for
## TFs will be skipped is displayed.
exampleObject <- step7TopGenesTADTables(
    TENETMultiAssayExperiment = exampleObject,
    TADFiles = exampleTADRegions,
    hypomethGplusAnalysis = TRUE,
    hypermethGplusAnalysis = FALSE,
    topGeneNumber = 10
)

## Results are included in the metadata of the returned MultiAssayExperiment
## object under the step7TopGenesTADTables list.
## For each analysis type, results are included in sub-lists, each
## of which contains results in the form of a data frame for topGenes and
## topTFs, unless the top genes are all TFs, in which case
## the separate top TFs output is skipped.
class(
    metadata(
        exampleObject
    )$step7TopGenesTADTables$hypomethGplusResults$topGenes
)

## Display results for selected hypomethylated RE DNA methylation sites. A
## variety of data are included for each RE DNA methylation site, including its
## location, the top genes it is linked to, and information on the count and
## identities of other genes found within the same TAD of the RE DNA methylation
## site. Note: A significant amount of data is output, so selected values are
## shown here.
head(
    metadata(
        exampleObject
    )$step7TopGenesTADTables$hypomethGplusResults$topGenes[
        c(1:6, 16:17)
    ]
)
```

## `step7TopGenesUCSCBedFiles`: Create bed-formatted interact files which can be loaded on the UCSC Genome Browser to display links between top genes and transcription factors and their linked RE DNA methylation sites

This function will output the interact file in the specified folder. For the
purposes of this example, we will save the output file in a temporary folder.
This interact file can be uploaded to the UCSC Genome Browser to visualize the
links between each of the top genes/TFs and their linked RE DNA methylation
sites.

```{r}
## Get the path to a temporary directory in which to save the output interact
## file
tempDirectory <- tempdir()

## Run the step7TopGenesUCSCBedFiles function.
## Since we performed analyses only using TFs in the step 3 function, the
## top genes are all TFs, so a message that separate output for
## TFs will be skipped is displayed.
filePaths <- step7TopGenesUCSCBedFiles(
    TENETMultiAssayExperiment = exampleObject,
    outputDirectory = tempDirectory,
    hypomethGplusAnalysis = TRUE,
    hypermethGplusAnalysis = FALSE,
    topGeneNumber = 10
)

## Unlike other functions, this function does not return the
## given TENETMultiAssayExperiment with additional information generated by the
## function in its metadata, but rather returns an object with information on
## where to find the created interact file.

## Get the path of the output file for top genes with hypomethylated G+ links
bedPath <- filePaths$hypoGplus$topGenes

## Read the first few lines of the file.
## The file largely contains information about each RE DNA methylation site-gene
## link, with additional information in the first line which allows it to be
## loaded by the UCSC Genome Browser.
cat(head(readLines(bedPath)), sep = "\n")

## Delete the output file, since this is just an example.
## Do not use this line if running on real data, as it will delete your created
## file.
invisible(file.remove(unlist(bedPath)))
```

## `step7TopGenesUserPeakOverlap`: Identify if RE DNA methylation sites linked to top genes/transcription factors are located within a specific distance of specified genomic regions

The `step7TopGenesUserPeakOverlap` function requires the user to supply either
a path to a directory containing bed-like files with peaks of interest, or a
single peak object given as a GRanges, data frame, or matrix object, as seen in
the example below. To illustrate the use of this function, we will use an
example peak dataset from the TENET.ExperimentHub package.

```{r}
## Load the example peak dataset from the TENET.ExperimentHub package
examplePeakFile <- TENET.ExperimentHub::exampleTENETPeakRegions()

## Peak files for this function must include the chromosome of each peak region
## in the first column, and the start and end positions of each peak in the
## second and third columns respectively. Additional columns can be included,
## but are not considered in this function.
class(examplePeakFile)
head(examplePeakFile)
```

The unique RE DNA methylation sites linked to the top genes will be overlapped
with the peak files, with a specified buffer region added to the RE DNA
methylation sites, so RE DNA methylation sites can be found in the vicinity of
peaks, rather than directly inside of them.

```{r}
## Run the step7TopGenesUserPeakOverlap function.
## Since we performed analyses only using TFs in the step 3 function, the
## top genes are all TFs, so a message that separate output for
## TFs will be skipped is displayed.
exampleObject <- step7TopGenesUserPeakOverlap(
    TENETMultiAssayExperiment = exampleObject,
    peakData = examplePeakFile,
    hypermethGplusAnalysis = FALSE,
    hypomethGplusAnalysis = TRUE,
    topGeneNumber = 10,
    distanceFromREDNAMethylationSites = 100
)

## Results are included in the metadata of the returned MultiAssayExperiment
## object under the step7TopGenesSurvival list.
## For each analysis type, data frames with peak overlap information are
## included in sub-lists, with each data frame saved under the names
## topGenes and topTFs, unless the top genes are all TFs, in which case
## the separate top TFs output is skipped.

## Display the data frame of results for RE DNA methylation sites linked to the
## top TFs. A variety of data are included for each RE DNA methylation site,
## including its location, the coordinates of its search window, the top genes
## it is linked to, and whether it was found within the specified distance to
## any peak in each of the peak files.
head(
    metadata(
        exampleObject
    )$step7TopGenesUserPeakOverlap$hypomethGplusResults$topGenes
)
```

# Datasets included in the TENET package

The following objects are contained in the TENET package. Since LazyData is not
enabled, objects will need to be accessed using the `data()` function, as
demonstrated below.

## `humanTranscriptionFactorList`: Human transcription factor list

A character vector of gene Ensembl IDs which were identified as human TFs by
Lambert SA et al (PMID: 29425488). Candidate proteins were manually examined by
a panel of experts based on available data. Proteins with experimentally
demonstrated DNA binding specificity were considered TFs. Other proteins, such
as co-factors and RNA binding proteins, were classified as non-TFs.
**Citation:** Lambert SA, Jolma A, Campitelli LF, et al. The Human
Transcription Factors. Cell. 2018 Feb 8;172(4):650-665. doi:
10.1016/j.cell.2018.01.029. Erratum in: Cell. 2018 Oct 4;175(2):598-599. PMID:
29425488.

```{r}
## Load the humanTranscriptionFactorList dataset
data("humanTranscriptionFactorList", package = "TENET")
## Display the names of the first few TFs on the list
head(humanTranscriptionFactorList)
```

# Session info

```{r}
sessionInfo()
```