---
title: "Package datasets"
author: "Tyler Grimes"
date: "`r Sys.Date()`"
output:
  html_vignette:
    number_sections: yes
    self_contained: yes
    toc: true 
vignette: >
  %\VignetteIndexEntry{Package datasets}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r setup, message=FALSE, warning=FALSE, echo=FALSE}
reset_options_scipen <- getOption("scipen")
reset_options_digits <- getOption("digits")
options(scipen = 1, digits = 5)
library(dnapath)
```

This vignette shows how the package datasets were created. The package data includes a list of cancer-related Reactome pathways, `p53_pathways`, and an example gene expression dataset, `glioma`. These are provided to be used in examples and for illustrative purposes.

# Reactome pathways 

A list of Reaectome pathways is provided for Homo sapiens (human). This list was initialized using the `dnapath::get_reactome_pathways()` function.

```{r}
pathway_list <- dnapath::get_reactome_pathways("human")
```

This list contains the gene sets for `r length(pathway_list)` pathways. Here are the first two pathways:
```{r}
pathway_list[1:2]
```

The list names correspond to the Reactome pathway name. Each pathway provides the entrezgene IDs of all genes in that pathway. For example, from the output above, we see that the first pathway in the list is the "purine salvage" pathway, and it contains the genes: 100, 132, 161823, etc.

In practice, all of these pathways can be used in the `dnapath()` analysis. However, in some studies we may only be interested in a subset of this list. For example, if we are studying a cancer-related gene expression dataset, we may have a specific interest in cancer-related pathways. There are two ways to obtain a subset from the pathway list: (1) using the pathway names or (2) using the pathway genes. These two approaches are demonstrated in turn.



## Subsetting pathways using pathway names

Suppose we are interested in pathways that are related to "p53" or "pi3k" regulatory processes (these are known cancer-related pathways). One option is to subset pathways based on whether or not the pathway names include these keywords. The `grep()` function is used to do this:

```{r}
# Obtain index of pathways containing "p53" or "pi3k" in their name.
index_pathways <- grep("((p53)|(pi3k))", names(pathway_list),
                       ignore.case = TRUE)
```

`grep` uses the regular expression "((p53)|(pi3k))", which says to identify all pathway names that include the phrase "p53" or "pi3k". However, we don't want to leave out pathways that use capital "P53" or "PI3K", so we set `ignore.case = TRUE` to ignore the case of the pathway names.

```{r}
# Subset onto the "p53" and "pi3k" pathways.
cancer_pathways <- pathway_list[index_pathways]
cancer_pathways[1:3] # Print out the first three
```

Subsetting on these cancer-related pathways results in `r length(cancer_pathways)` pathways. 

In general, the pathway list obtained from the `dnapath::get_reactome_pathways()` function can be subset based on the pathway names, as illustrated in this example. Here, we focused on cancer pathways, but the regular expression "((p53)|(pi3k))" can be replaced to search for any desired pathways. To learn more about the `grep()` function and regular expression, see [chapter 17](https://bookdown.org/rdpeng/rprogdatascience/regular-expressions.html) of R Programming for Data Science by Roger Peng.




## Subsetting pathways using genes

As in the previous example, suppose we are interested in pathways that are related to "p53" or "pi3k". However, this time we want to identify *all* pathways that involve the genes, not just pathways that are named after them. (Note that some pathways may contain the p53 gene but not have "p53" in its name.)

Since the pathway list uses entrezgene IDs, the first step is to find the entrezgene ID for each gene of interest. The ID for P53 is [7157](https://www.ncbi.nlm.nih.gov/gene?cmd=Retrieve&dopt=full_report&list_uids=7157). PI3K actually refers to a family of enzymes; one class of genes in this family include PIK3CA, PIK3CB, PIK3CD, and PIK3CG, which have the entrezgene IDs
[5290](https://www.ncbi.nlm.nih.gov/gene?cmd=Retrieve&dopt=full_report&list_uids=5290), 
[5291](https://www.ncbi.nlm.nih.gov/gene?cmd=Retrieve&dopt=full_report&list_uids=5291),
[5293](https://www.ncbi.nlm.nih.gov/gene?cmd=Retrieve&dopt=full_report&list_uids=5293), and
[5294](https://www.ncbi.nlm.nih.gov/gene?cmd=Retrieve&dopt=full_report&list_uids=5294).

```{r}
# Make a vector of all entrezgene IDs that we are interested in.
genes_of_interest <- c(7157, 5290, 5291, 5293, 5294)
# Obtain index of pathways containing any genes of interest.
index_pathways <- which(sapply(pathway_list, function(pathway) {
  any(genes_of_interest %in% pathway)
}))
```

The `sapply()` function is used to loop over each pathway in the pathway list, and for each we determine if any of the genes of interest are in the pathway's gene set. We end up with a set of indices identifying which pathways contain at least one of the five targeted genes.

```{r}
# Subset onto pathways containing p53 gene or pi3k family of genes.
cancer_pathways <- pathway_list[index_pathways]
cancer_pathways[1:3] # Print out the first three
```

Subsetting results in `r length(cancer_pathways)` pathways that contain at least one of the five cancer-related genes of interest. In general, any set of entrezgene IDs can be used to subset the pathways obtained from the `dnapath::get_reactome_pathways()` function, as illustrated in this example. 



## `p53_pathways` data

The `dnapath` package provides an example pathway list named `p53_pathways`. This is a small list provided for illustrative purposes and for use in examples. It includes T53 Reactome pathways, which was created using the steps described previously.

```{r}
# Obtain index of pathways containing "p53" in their name.
index_pathways <- grep("(p53)", names(pathway_list),ignore.case = TRUE)
p53_pathways <- pathway_list[index_pathways]
```

The dataset was compressed and saved for the package using the code:
```R
# Not run: Only for building the package.
usethis::use_data(p53_pathways)
tools::resaveRdaFiles(file.path("data"), compress = "auto")
tools::checkRdaFiles(file.path("data"))
```

# Gene expression example

A small gene expression dataset is provided in `dnapath`. This example data is a Mesothelioma dataset containing gene expression and clinical data that were generated by [The Cancer Genome Atlas](https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga) (TCGA) and is downloaded using the [LinkedOmics portal](http://www.linkedomics.org/data_download/TCGA-GBMLGG/). 

## `meso` data

The first step is to download and process the clinical and gene expression dataset from LinkedOmics.
```{r message = FALSE}
file_clinical <- paste0("http://linkedomics.org/data_download/TCGA-MESO/Human__",
                        "TCGA_MESO__MS__Clinical__Clinical__01_28_2016__BI__",
                        "Clinical__Firehose.tsi")
file_rnaseq <- paste0("http://linkedomics.org/data_download/TCGA-MESO/Human__",
                      "TCGA_MESO__UNC__RNAseq__HiSeq_RNA__01_28_2016__BI__Gene",
                      "__Firehose_RSEM_log2.cct.gz")

clinical <- readr::read_tsv(file_clinical)
rnaseq <- readr::read_tsv(file_rnaseq)
```

### Clinical dataset

The differential network analysis performed using `dnapath()` compares the gene-gene association structure between two groups (populations). In the Mesothelioma dataset, the clinical data can be used to identify two subgroups of interest. In this example dataset, the tumor stage (classified as stage i, ii, iii, or iv) will be used to partition the data. In practice, `dnapath()` could be used make pair-wise comparisons across these tumor stages, or perhaps some stages can be combined together. But for the purpose of creating a small, example dataset, we will only look at tumor stages ii and iv.

The first step is to process the downloaded data so that it contains one column of sample IDs and one column of tumor stage. 

```{r}
# First, we can glimpse at the raw data to see how it is structured.
clinical[1:5, 1:3]
```

```{r}
# The columns contain samples, and rows contain different variables
# We want the transpose of this.
variable_names <- clinical$attrib_name # Store the variable names.

# Transpose the matrix so that columns correspond to variables.
clinical <- t(clinical[, -1]) # Don't include the "attrib_name" column as a row.
colnames(clinical) <- variable_names # Rename columns using variable names.

# Glimpse at the transposed dataset:
clinical[1:4, 1:4]
```

```{r}
# Now that the rows correspond to subjects and columns to variables,
# the next step is to add a new column that contains the sample IDs. 

clinical <- cbind(ID = rownames(clinical), clinical)
rownames(clinical) <- NULL # Remove the row names.

# Glimpse at the data.
clinical[1:4, 1:4]
```

```{r}
# For the example dataset, we are only interested in ID and tumor stage.
clinical <- clinical[, c("ID", "pathologic_stage")]
# The two groups that will be compared are stage 2 and stage 4.
# Subset rows onto those that have stage 2 or stage 4.
clinical <- clinical[clinical[, "pathologic_stage"] %in% c("stageii", "stageiv"), ]
clinical[1:5, ]
```
```{r}
# We are left with only stage 2 and stage 4 samples:
table(clinical[, 2])
```

The final dataset contains 16 samples for group 1 (tumor stage ii) and 16 samples for group 2 (tumor stage iv). Note that the groups don't need equal sample sizes, but they happen to be equal in these data.

### RNA-seq dataset

The second step is to process the gene expression dataset. The LinkedOmics portal provides normalized counts from the RNA-seq experiments (the data are not raw reads, which must be aligned and annotated in order to obtain gene expression counts). For the purpose of creating a small example dataset, these data will be subset onto the 32 samples obtained in the previous section, and onto the 160 genes in the `cancer_pathway` list.

```{r}
# First, we can glimps at the raw data to see how it is structured.
rnaseq[1:5, 1:5]
```

```{r}
# As with the clinical data, we need to transpose the raw data so that
# samples are in the rows.
gene_symbols <- rnaseq$attrib_name # Store the gene names

rnaseq <- t(rnaseq[, -1])
colnames(rnaseq) <- gene_symbols

rnaseq[1:5, 1:5]
```

```{r}
# Check if all 32 clinical samples contain gene expression data.
all(clinical[, "ID"] %in% rownames(rnaseq))

# Subset the rnaseq data onto those in the clinical data
rnaseq <- rnaseq[rownames(rnaseq) %in% clinical[, "ID"], ]

# Finally, we must make sure the two dataset are aligned.
# There are many ways to do this, here is one:
#   For each ID in the clinical dataset, find the corresponding row in rnaseq.
#   The rows are then reordered to match the IDs in clinical.
rnaseq <- rnaseq[sapply(clinical[, "ID"], 
                        function(ID) which(rownames(rnaseq) == ID)), ]

# Check that the IDs match:
all(rownames(rnaseq) == clinical[, "ID"])

rnaseq[1:5, 1:5]
```

The final step of processing the gene expression data is to subset the columns to the genes contained in the `cancer_pathway` list. However, the raw data provide gene symbols rather than entrezgene IDs (which the Reactome pathways use). So, we must either relabel the gene symbols in the `rnaseq` dataset to entrezgene IDs, or relabel the entrezgene IDs in `cancer_pathway` to contain gene symbols. Either of these operations can be easily carried out using the `symbol_to_entrez()`, `entrez_to_symbol()`, and `rename_genes()` functions provided in the package. These two approaches are briefly demonstrated next.

> **Note:** Internet connection is required to connect to biomaRt, which the `entrez_to_symbol` and `symbol_to_entrez` methods use to map entrezgene IDs and gene symbols. The `dir_save` argument can be set when running these to save the ID mapping obtained from biomaRt in the specified directory. This way, the mapping can be obtained once (using the internet connection) and accessed from memory in all future calls of these functions. The temporary directory `tempdir()` is used in this example.

First, we show how to convert gene symbols in the gene expression dataset to entrezgene IDs.

```{r}
# Convert gene symbols -> entrezgene IDs
gene_symbols <- colnames(rnaseq) # Extract the gene symbols to relabel.
gene_mat <- symbol_to_entrez(gene_symbols, "human",
                             dir_save = tempdir()) # Obtain mapping.
gene_mat[1:5, ]
```

```{r warning = FALSE}
# The rename_genes() method is a multipurpose function that can be used to rename 
# the genes in a vector (the current case), a pathway list, or the 'dnapath_list'
# object that is returned from dnapath() after performing the differential
# network analysis.
gene_entrez <- rename_genes(gene_symbols, gene_mat) # Rename the symbols.
colnames(rnaseq) <- gene_entrez # Update the columns using the entrezgene IDs.
rnaseq[1:5, 1:5] # Columns now contain entrezgene IDs
```

This is the approach used to create the `meso` data that is provided in this package. In general, it is advised to use entrezgene IDs throughout the analysis pipeline, and **only convert to symbols at the very end**. Entrezgene IDs are unique identifiers for each gene, whereas gene symbols can sometimes be ambiguous. This ambiguity can lead to information loss if used early on the analysis. For example, notice the third column in above output, the "A2BP1" gene, was not mapped to an entrezgene ID. This is because "A2BP1" is an alias for the RBFOX1 gene. The `symbol_to_entrez()` function obtains a mapping for RBFOX1 but not for the alias A2BP1. Because an entrezgene ID was not identified for A2BP1, it will not be identified in any of the pathways during the differential network analysis. 

Gene symbols should only be considerd at the end of the analysis. It can be conventient to rename the IDs into symbols, since symbols are more recognizable for viewing and interpreting the results. 

**Important note:** In the Mesothelioma dataset, we are utilizing summarized gene expression counts as opposed to raw RNA-sequencing reads. If the raw reads are unavailable, then the user may not be able to control whether or not the genes are entrezgene IDs or gene symbols. In this case, we must be aware of the limitations and take caution if many gene symbols fail to map to entrezgene IDs. However, if raw reads are available, then it is recommended to align and annotate the reads using entrezgene IDs when obtaining the expression counts.

To finish the processing, the columns of `rnaseq` are subset onto the genes in `p53_pathways`.
```{r}
# Subset the columns onto only those genes contained in 'p53_pathways'
index_genes <- which(colnames(rnaseq) %in% get_genes(p53_pathways))
rnaseq <- rnaseq[, index_genes]
dim(rnaseq) # 32 samples with 150 genes.
```

For completeness, we will show how to convert the entrezgene IDs in the pathway list into gene symbols.

```{r}
# Convert entrezgene IDs -> gene symbols
gene_entrez <- get_genes(p53_pathways) # Extract genes from pathway list.
gene_mat <- entrez_to_symbol(gene_entrez, "human", 
                             dir_save = tempdir()) # Obtain mapping.
gene_mat[1:5, ]
```

```{r}
# Convert the entrezgene IDs into gene symbols in the pathway list.
new_pathway_list <- rename_genes(p53_pathways, gene_mat) 
new_pathway_list[1:2] # Print the first two pathways.
```

### Finalizing

The `dnapath()` function requires either a list of two gene expression datasets (each corresponding to a different group), or a single gene expression data and a 'group' vector that indicates which rows belong to which group. The `meso` data is formatted using the second option. The current `rnaseq` dataset that we have processed contains samples from both groups; the group information can be obtained from the `clinical` dataset. These two pieces of information can be conveniently stored as a list.

```{r}
meso <- list(gene_expression = rnaseq,
             groups = clinical[, "pathologic_stage"])
```

The list was compressed and saved for the package using the code:
```R
# Not run: Only for building the package.
usethis::use_data(meso)
tools::resaveRdaFiles(file.path("data"), compress = "auto")
tools::checkRdaFiles(file.path("data"))
```

```{r echo=FALSE}
options(scipen = reset_options_scipen, digits = reset_options_digits)
```

# biomaRt 

The `biomaRt` package is used to link gene symbols with entrezgene IDs. The `entrez_to_symbol()` and `symbol_to_entrez()` functions require an internet connection to connect with the biomaRt database.

This package contains an instance of the mapping between HGNC gene symbols and entregene IDs, which will be used if internet access is not available or if the `biomaRt` package is not installed. However, this dataset will not be applicable to species other than *hsapiens* and may not contain the most up-to-date mapping. A warning is given when it is used.

```R
# Not run: Internet access and `biomaRt` package is required.
mart <- biomaRt::useMart(biomart = "ensembl", 
                         dataset = "hsapiens_gene_ensembl")

biomart_hsapiens <- biomaRt::getBM(attributes = c("entrezgene_id", "hgnc_symbol"),
                                   mart = mart)

library(dplyr) # For pipe operator "%>%".
biomart_hsapiens %>%
  dplyr::filter(!is.na(entrezgene_id)) %>%
  dplyr::group_by(entrezgene_id) %>%
  dplyr::summarise(across(everything(), dplyr::first)) ->
  biomart_hsapiens
```

This mapping was compressed and saved for the package using the code:

```R
# Not run: Only for building the package.
usethis::use_data(biomart_hsapiens)
tools::resaveRdaFiles(file.path("data"), compress = "auto")
tools::checkRdaFiles(file.path("data"))
```