rat-endurance-6m.Rmd

---
title: "Workshop: MoTrPAC R Packages from the Endurance Training in Young Rats Study"
output: 
  rmarkdown::html_document:
    toc: true
    code_folding: show 
    toc_depth: 4
    toc_float:
      collapsed: false
      smooth_scroll: false
    theme: sandstone
    highlight: tango
editor_options: 
  markdown: 
    wrap: 72
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = 'center'
)
```

# Introduction

The Molecular Transducers of Physical Activity Consortium (MoTrPAC)
conducted the first comprehensive, organism-wide, multi-omic
investigation of endurance exercise training effects in rats. This
large-scale study involved six-month-old male and female Fischer 344
rats subjected to progressive treadmill endurance exercise training over
1, 2, 4, or 8 weeks. Tissues were collected from these rats 48 hours
after the last exercise bout to capture the adaptations resulting from
the training, with sex-matched sedentary, untrained rats serving as
controls. Utilizing multiple 'omic' platforms including transcriptomics,
epigenomics, proteomics, and metabolomics—across a wide range of
tissues, the study provides a rich dataset that captures the temporal
effects of endurance exercise training at multiple biological levels.
This dataset offers invaluable insights into the molecular mechanisms
underlying exercise-induced adaptations. Find the open access
[manuscript published on
Nature](https://www.nature.com/nature/volumes/629/issues/8010).

To facilitate access to and analysis of this extensive dataset, two
essential R packages have been developed:

-   [MotrpacRatTraining6moData](https://motrpac.github.io/MotrpacRatTraining6moData/):
    Provides access to the processed data and downstream analysis
    results from the MoTrPAC endurance exercise study.
-   [MotrpacRatTraining6mo](https://motrpac.github.io/MotrpacRatTraining6mo/):
    Offers functions to help retrieve and explore the data, enabling
    users to perform analyses and reproduce key findings from the study.

In this tutorial, we will explore how to utilize the
`MotrpacRatTraining6moData` package to access and navigate the extensive
dataset generated by the endurance exercise study, and employ the
`MotrpacRatTraining6mo` package to perform analyses and gain insights
into the molecular responses to endurance training in young rats.

**Note**: This notebook is based on the vignettes available in
[MotrpacRatTraining6moData](https://motrpac.github.io/MotrpacRatTraining6moData/articles/MotrpacRatTraining6moData.html)
and [MotrpacRatTraining6mo
vignette](https://motrpac.github.io/MotrpacRatTraining6mo/articles/MotrpacRatTraining6mo.html).

## Installations

Install both packages with `devtools`:

``` r
if (!require("devtools", quietly = TRUE)){
  install.packages("devtools")
}
options(timeout=1e5) # extend the timeout for downloading large files

# First the data package, which might take a long time
devtools::install_github("MoTrPAC/MotrpacRatTraining6moData")

# Then the analysis package:
devtools::install_github("MoTrPAC/MotrpacRatTraining6mo")
```

Once we install this package, we can load these and other libraries used
in this tutorial:

```{r setup}
library(MotrpacRatTraining6moData)
library(MotrpacRatTraining6mo)
library(ggplot2) # for plots in this tutorial 
library(dplyr)
library(impute)  # For KNN imputation
library(tidyverse)
library(lubridate)
library(pheatmap)
library(viridis)
```

# Study design: the PHENO data object

Details of the experimental design can be found in the [supplementary
methods of the Nature
paper](https://static-content.springer.com/esm/art%3A10.1038%2Fs41586-023-06877-w/MediaObjects/41586_2023_6877_MOESM1_ESM.pdf).

Briefly, 6-month-old young adult rats were subjected to progressive
endurance exercise training for 1, 2, 4, or 8 weeks, with tissues
collected 48 hours after the last training bout. Sex-matched sedentary,
untrained rats were used as controls. Whole blood, plasma, and 18 solid
tissues were analyzed using genomics, proteomics, metabolomics, and
protein immunoassay technologies, with most assays performed in a subset
of these tissues. Depending on the assay, between 3 and 6 replicates per
sex per time point were analyzed.

The **PHENO** object represents phenotypic data from the MoTrPAC
endurance exercise training study in six-month-old rats. It is a
comprehensive dataset containing **5,955 rows and 510 variables**, with
each row corresponding to a unique sample identified by `viallabel`. The
dataset captures detailed information on the animals, their training
regimens, specimen collection, and various physiological metrics. Here
is a summary of the major components:

### Categories of Variables in the PHENO Dataset

1.  **Identifiers**: Unique identifiers such as `pid`, `bid`,
    `viallabel`, and `labelid` are used to link samples to individual
    animals and their specimen labels.

2.  **Animal Information**: Variables such as `registration.d_birth`,
    `registration.sex`, and `registration.weight` provide basic details
    about the rats, including birth date, sex, and weight.

3.  **Registration and Housing**: Data about the arrival, housing
    conditions, registration (`registration.d_arrive`,
    `registration.cagenumber`), and light conditions
    (`registration.d_reverselight`) are included to document the
    environmental conditions in which animals were kept.

4.  **Intervention and Randomization**: Information about intervention
    groups (e.g., control or training) is captured in variables like
    `key.intervention`, `key.anirandgroup`, and `key.protocol`,
    indicating how each rat was assigned to different treatment
    conditions.

5.  **Familiarization and Training Data**: Detailed records of treadmill
    familiarization (`familiarization.d_treadmillbegin`) and progressive
    endurance training are included. Training data cover up to **40
    days**, with variables for each day’s treadmill speed, incline, and
    exercise time (`training.dayX_treadmillspeed`,
    `training.dayX_timeontreadmill`).

6.  **VO2 Max and NMR Testing**: Information from **VO2 max tests** and
    **NMR body composition tests** is included. Variables such as
    `vo2.max.test.vo2_max` and `nmr.testing.nmr_fat` capture maximum
    oxygen uptake, percent body fat, and other key metabolic measures.

7.  **Specimen Collection**: Specimen collection details are provided,
    including dates (`specimen.collection.d_visit`), times of collection
    (`specimen.processing.t_collection`), and types of tissues collected
    (`specimen.processing.sampletypedescription`).

8.  **Terminal Measures**: Terminal metrics, such as body weight at
    sacrifice and tissue weights (`terminal.weight.bw`,
    `terminal.weight.sol`), document the final physical characteristics
    of the animals after completion of the study.

9.  **Calculated Variables**: Derived variables, such as changes in body
    composition (`calculated.variables.pct_body_fat_change`), lactate
    levels, and `vo2_max`, provide insights into physiological
    adaptations resulting from the exercise intervention.

10. **Custom Variables**: Variables like `sacrificetime`,
    `intervention`, and `time_to_freeze` have been added for simplified
    grouping and analysis, reflecting key time points, interventions,
    and time-to-processing of samples.

### Key Variables for Analysis

-   **`pid`**: Participant Identifier. A unique 8-digit identifier
    assigned to each animal subject (rats in this study). All samples
    (`viallabel`) coming from the same animal have the same `pid`. This
    variable is crucial when combining results from multiple assays with
    phenotypic data, allowing for animal-specific longitudinal analysis.
-   **`viallabel`**: Vial Label ID. A unique 11-digit code assigned to
    each sample vial. This ID is present across all related results and
    metadata files, serving as a key to link the quantitative results
    with the phenotypic data.
-   **`sex`**: Sex of the animal. Represented as `"female"` for Female
    and `"male"`. This variable is critical for any analysis that aims
    to determine sex-based differences.
-   **`group`**: A simplify version of `study_group_timepoint` that only
    includes the group and time points: `control` (sedentary), and `1w`,
    `2w`, `4w`, and `8w` (exercise)
-   **`tissue`**: tissue name (short)

In summary, the **PHENO** object provides a comprehensive overview of
each animal's involvement in the study, including the conditions under
which they were kept, the specific training they underwent, and the
physiological changes observed as a result of the endurance training
intervention. This dataset allows for in-depth analysis of the effects
of physical activity at various biological levels.

### Visualizations of Phenotypes

#### Animal counts

First, let's select the ids of the rats used in this study:

```{r}
# Load the proteomics data, for example:
data_prot <- MotrpacRatTraining6mo::combine_normalized_data(
  assays = c("PROT","UBIQ","PHOSPHO","ACETYL"), 
  exclude_outliers = TRUE, )

# Get the unique pid (rat id) values from the proteomics data
rat_omics <- setdiff(names(data_prot), c("feature", "feature_ID", "tissue", "assay"))
```

```{r}

# Count unique 'pid' by 'group' and 'sex'
count_data <- PHENO %>%
  group_by(group, sex) %>%
  summarise(unique_pid_count = n_distinct(pid),
            .groups = 'drop') %>%
  ungroup()

# Create the grouped bar chart
p <- ggplot(count_data, aes(x = group, y = unique_pid_count, fill = sex)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Animal count",
       subtitle = "All animals processed at the Animal Site",
       x = "Study Group Time Point",
       y = "Count of Unique Animals (pid)",
       fill = "Sex") +
  theme_linedraw() +
  coord_flip()

print(p)

#  Animals used in this study

# Filter by pid (i.e., the animal id)
count_data <- PHENO %>%
  filter(pid %in% rat_omics) %>%
  group_by(group, sex) %>%
  summarise(unique_pid_count = n_distinct(pid),
            .groups = 'drop') %>%
  ungroup()

# Create the grouped bar chart
p2 <- ggplot(count_data, aes(x = group, y = unique_pid_count, fill = sex)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Animal Count",
       subtitle = "Only animals used for molecular profile",
       x = "Study Group Time Point",
       y = "Count of Unique Animals (pid)",
       fill = "Sex") +
  theme_linedraw() +
  coord_flip()

print(p2)

```

#### Sample counts by tissue

```{r}
# Count total number of unique 'viallabel' by 'tissue'
viallabel_data <- PHENO %>%
  filter(pid %in% rat_omics) %>% # only include animals used in the omics data
  group_by(tissue, group) %>%
  summarise(total_viallabel_count = n_distinct(viallabel), .groups = 'drop') %>%
  ungroup()

p <- ggplot(viallabel_data, aes(x = tissue, y = total_viallabel_count, fill = group)) +
  geom_bar(stat = "identity") +
  labs(title = "Total Number of Samples collected for Each Tissue",
       x = "Tissue",
       y = "Total Count of Vial Labels",
       fill = "Group") +
  theme_minimal() +
  theme(legend.position = "right") + coord_flip()

print(p)

```

#### Weight distributions

```{r}
# Select variables from PHENO that start with "training" and end with "_weight"
weight_cols <- grep("^training.*_weight$", names(PHENO), value = TRUE)

# Select those columns along with "pid", "sex", and "group"
selected_cols <- c("pid", "sex", "group", weight_cols)

# Get unique values
unique_pheno <- unique(PHENO[, selected_cols])

# Reshape data to long format for plotting
pheno_long <- unique_pheno %>%
  dplyr::filter(pid %in% rat_omics) %>%
  tidyr::pivot_longer(cols = starts_with("training"), names_to = "training_session", values_to = "weight")

# Define the desired order of training sessions
training_session_order <- c("training.day1_weight", 
                            "training.day6_weight", 
                            "training.day11_weight", 
                            "training.day16_weight", 
                            "training.day21_weight", 
                            "training.day26_weight", 
                            "training.day31_weight", 
                            "training.day36_weight")

# Convert training_session to a factor with the specified order
pheno_long$training_session <- factor(pheno_long$training_session, 
                                      levels = training_session_order)

# Plot the jittered plot for each training session, faceted by sex and colored by group
p1 <- ggplot(pheno_long, aes(x = training_session, y = weight, color = group)) +
  geom_jitter(width = 0.2, size = 2, alpha = 0.6, na.rm = TRUE) +
  facet_wrap(~sex) +
  theme_bw() +
  labs(title = "Weight by Training Session",
       x = "Training Session",
       y = "Weight") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

print(p1)

p2 <- ggplot(pheno_long, aes(x = training_session, y = weight, color = group)) +
  geom_boxplot(alpha = 0.3, outlier.shape = NA, na.rm = TRUE, linewidth = 0.5) +
  facet_wrap(~sex, scales = "free") +
  theme_bw() +
  labs(title = "Weight by Training Session",
       x = "Training Session",
       y = "Weight") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

print(p2)

p3 <- ggplot(pheno_long, aes(x = training_session, y = weight, group = group, color = group, fill = group)) +
  stat_summary(fun = median, geom = "point", size = 4, shape = 21, alpha = 0.7, na.rm = TRUE) +
  stat_summary(fun = median, geom = "line", size = 1, alpha = 0.7, na.rm = TRUE) +
  facet_wrap(~sex, scales = "free") +
  theme_linedraw() +
  labs(title = "Mean Weight by Training Session",
       x = "Training Session",
       y = "Mean Weight") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

print(p3)

```

#### VO2max

VO2max measurements were taken at the start for all time points, and at the end for the 4-week, 8-week, and control training sessions.

```{r}
# colnames(PHENO)[grep("vo2*d_visit", colnames(PHENO))]

# Select those columns along with "pid", "sex", and "group"
selected_cols <- c("pid", 
                   "sex", 
                   "group", 
                   "vo2.max.test.vo2_max_1", 
                   "vo2.max.test.vo2_max_2")

# Get unique values
unique_pheno <- unique(PHENO[, selected_cols])

# Reshape data to long format for plotting
pheno_long <- unique_pheno %>%
  dplyr::filter(pid %in% rat_omics) %>%
  tidyr::pivot_longer(cols = starts_with("vo2"), 
                      names_to = "training_session", 
                      values_to = "vo2max")

p1 <- ggplot(pheno_long, 
             aes(x = training_session, y = vo2max, color = group)) +
  geom_boxplot(alpha = 0.3, outlier.shape = NA, na.rm = TRUE) +
  # geom_jitter(width = 0.2, size = 2, alpha = 0.6, na.rm = TRUE) +
  facet_wrap(~sex, scales = "free") +
  theme_bw() +
  labs(title = "VO2max distributions at first and last training session",
       x = "Training Session",
       y = "VO2max") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

print(p1)

p2 <- ggplot(pheno_long, aes(x = training_session, y = vo2max, group = group, color = group, fill = group)) +
  stat_summary(fun = mean, geom = "point", size = 4, shape = 21, alpha = 0.7, na.rm = TRUE) +
  stat_summary(fun = mean, geom = "line", size = 1, alpha = 0.7, na.rm = TRUE) +
  facet_wrap(~sex) +
  theme_linedraw() +
  labs(title = "Mean VO2max at first and last training session",
       x = "Training Session",
       y = "Mean VO2max") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

print(p2)
```

#### Treadmill speed over time

```{r, fig.width=8, fig.height=10}

# Filter the data for training rats
PHENO_training <- PHENO %>%
  filter(intervention == "training")

# Pivot the date columns
PHENO_long_date <- PHENO_training %>%
  select(pid, group, matches("training\\.day\\d+date$")) %>%
  pivot_longer(
    cols = -c(pid, group),
    names_to = "day",
    names_pattern = "training\\.day(\\d+)date$",
    values_to = "date"
  ) %>%
  mutate(
    day = as.integer(day),
    date = dmy(date)
  )

# Pivot the treadmillspeed columns
PHENO_long_speed <- PHENO_training %>%
  select(pid, group, matches("training\\.day\\d+_treadmillspeed$")) %>%
  pivot_longer(
    cols = -c(pid, group),
    names_to = "day",
    names_pattern = "training\\.day(\\d+)_treadmillspeed$",
    values_to = "treadmillspeed"
  ) %>%
  mutate(
    day = as.integer(day),
    treadmillspeed = as.numeric(treadmillspeed)
  )

# Merge the two datasets
PHENO_long <- suppressWarnings(left_join(PHENO_long_date, PHENO_long_speed, by = c("pid", "group", "day")) %>% distinct())

# Remove rows where both 'date' and 'treadmillspeed' are NA
PHENO_long <- PHENO_long %>%
  filter(!(is.na(date) & is.na(treadmillspeed)))

# Create a new Id combining group and pid
PHENO_long <- PHENO_long %>%
  mutate(group_pid = paste0(group, "_", pid))

# # Identify duplicates
# duplicates <- PHENO_long %>%
#   select(group_pid, date, treadmillspeed) %>%
#   group_by(group_pid, date) %>%
#   filter(n() > 1)
# 
# # View duplicates
# print(duplicates)

# Reshape the data to wide format
PHENO_wide <- PHENO_long %>%
  pivot_wider(
    names_from = date,
    values_from = treadmillspeed
  )

# Aggregate duplicates by taking the mean of treadmillspeed
PHENO_long <- PHENO_long %>%
  group_by(group_pid, date) %>%
  summarise(
    treadmillspeed = mean(treadmillspeed, na.rm = TRUE),
    .groups = 'drop'
  )

# Reshape the data to wide format
PHENO_wide <- PHENO_long %>%
  pivot_wider(
    names_from = date,
    values_from = treadmillspeed
  )

# Convert the data frame to a matrix
PHENO_matrix <- PHENO_wide %>%
  column_to_rownames(var = "group_pid") %>%
  as.matrix()

# Replace NA values with zeros
PHENO_matrix[is.na(PHENO_matrix)] <- 0

# Create the heatmap
pheatmap(
  PHENO_matrix,
  cluster_rows = FALSE,
  cluster_cols = FALSE,
  show_rownames = TRUE,
  show_colnames = TRUE,
  scale = "none",
  fontsize_row = 4,
  fontsize_col = 4,
  angle_col = 45,
  color = viridis::viridis(100),
  main = "Treadmill Speed Over Time for Training Rats"
)


```


# DATA

**Tip:** To learn more about any data object, use `?` to retrieve the
documentation, e.g., `?METAB_FEATURE_ID_MAP`. Note that
`MotrpacRatTraining6moData` must be installed and loaded with
`library()` for this to work.

## Tissue and assay abbreviations

It is important to be aware of the tissue and assay abbreviations
because they are used to name many data objects. The vectors of
abbreviations are also available in `TISSUE_ABBREV` and `ASSAY_ABBREV`.

### Tissues

-   **ADRNL:** adrenal gland
-   **BAT:** brown adipose tissue
-   **BLOOD:** whole blood
-   **COLON:** colon
-   **CORTEX:** cerebral cortex
-   **HEART:** heart
-   **HIPPOC:** hippocampus
-   **HYPOTH:** hypothalamus
-   **KIDNEY:** kidney
-   **LIVER:** liver
-   **LUNG:** lung
-   **OVARY:** ovaries (female gonads)
-   **PLASMA:** plasma from blood
-   **SKM-GN:** gastrocnemius (leg skeletal muscle)
-   **SKM-VL:** vastus lateralis (leg skeletal muscle)
-   **SMLINT:** small intestine
-   **SPLEEN:** spleen
-   **TESTES:** testes (male gonads)
-   **VENACV:** vena cava
-   **WAT-SC:** subcutaneous white adipose tissue

### Assays/omes

-   **ACETYL:** acetylproteomics; protein site acetylation
-   **ATAC:** chromatin accessibility, ATAC-seq data
-   **IMMUNO:** multiplexed immunoassays (cytokines and hormones)
-   **METAB:** metabolomics and lipidomics
-   **METHYL:** DNA methylation, RRBS data
-   **PHOSPHO:** phosphoproteomics; protein site phosphorylation
-   **PROT:** global proteomics; protein abundance
-   **TRNSCRPT:** transcriptomics, RNA-Seq data
-   **UBIQ:** ubiquitylome; protein site ubiquitination

## Summary of data

Here is a brief summary of the kinds of data included in the data
package:

-   Assay, tissue, sex, and training group abbreviations, codes, colors,
    and order used in plots
-   Phenotypic data, `PHENO`
-   Mapping between various feature identifiers, i.e.,
    `FEATURE_TO_GENE`, `RAT_TO_HUMAN_GENE`
-   Ome-specific feature annotation, i.e., `METAB_FEATURE_ID_MAP`,
    `METHYL_FEATURE_ANNOT` (GCP only), `ATAC_FEATURE_ANNOT` (GCP only),
    `PROT_FEATURE_ANNOT` (GCP only), `PHOSPHO_FEATURE_ANNOT` (GCP only),
    `UBIQ_FEATURE_ANNOT` (GCP only), `ACETYL_FEATURE_ANNOT` (GCP only),
    `TRNSCRPT_FEATURE_ANNOT` (GCP only)
-   Ome-specific sample-level metadata, i.e., `TRNSCRPT_META`,
    `ATAC_META`, `METHYL_META`, `IMMUNO_META`, `PHOSPHO_META`,
    `PROT_META`, `ACETYL_META`, `UBIQ_META`
-   Raw counts for RNA-Seq (TRNSCRPT), ATAC-Seq (ATAC), and RRBS
    (METHYL) data, e.g., `TRNSCRPT_LIVER_RAW_COUNTS`. Note that
    epigenetic data (ATAC and METHYL) must be downloaded from Google
    Cloud Storage. See more details
    [here](https://motrpac.github.io/MotrpacRatTraining6moData/index.html#access-epigenomics-data-through-google-cloud-storage).
-   Normalized sample-level data, e.g., `TRNSCRPT_SKMGN_NORM_DATA`
-   Differential analysis results, e.g., `HEART_PROT_DA`
-   Sample outliers excluded from differential analysis, `OUTLIERS`
-   Table of training-regulated features at 5% FDR,
    `TRAINING_REGULATED_FEATURES`
-   Bayesian graphical analysis inputs and results
-   Pathway enrichment of main graphical clusters, `GRAPH_PW_ENRICH`

A list of all the available data objects and a brief description are
available
[here](https://motrpac.github.io/MotrpacRatTraining6moData/index.html#explore-available-data-objects).

### Samples used in this study

```{r}
# List all data objects in the package
all_objects <- data(package = "MotrpacRatTraining6moData")$results[, "Item"]

# Filter the objects that contain the string "NORM" and exclude specific unwanted objects
exclude_objects <- c("IMMUNO_NORM_DATA_FLAT", "IMMUNO_NORM_DATA_NESTED", 
                     "METAB_NORM_DATA_FLAT", "METAB_NORM_DATA_NESTED", 
                     "TRAINING_REGULATED_NORM_DATA", "TRAINING_REGULATED_NORM_DATA_NO_OUTLIERS")

norm_objects <- all_objects[grepl("NORM_DATA", all_objects) & !(all_objects %in% exclude_objects)]

```

## Load the data

Use `load_sample_data()` to load sample-level data for a specific ome
and tissue. Here we fetch various forms of the sample-level RNA-Seq
(TRNSCRPT) data for brown adipose tissue (BAT) as an example.

```{r load sample data, eval=F}
# Load RNA-seq raw counts for brown adipose tissue 
data <- load_sample_data("BAT", "TRNSCRPT", normalized = FALSE)

# Load the normalized RNA-seq data for brown adipose tissue instead
data <- load_sample_data("BAT", "TRNSCRPT")

# Load the normalized RNA-seq data for brown adipose tissue, but exclude sample outliers 
data <- load_sample_data("BAT", "TRNSCRPT", exclude_outliers = TRUE)

# Load the normalized RNA-seq data for brown adipose tissue for training-regulated features only
data <- load_sample_data("BAT", "TRNSCRPT", training_regulated_only = TRUE)
```

`load_sample_data()` will tell you if the specified dataset doesn't
exist.

```{r load sample data 2}
data = load_sample_data("VENACV", "PROT")
```

`load_sample_data()` will also download epigenetic data from Google
Cloud Storage.

::: yellow
**Note:** Epigenetic data require substantially more memory than other
omes.
:::

```{r load epigen sample data, eval=F}
# Load ATAC-seq raw counts for hippocampus, excluding outliers 
data = load_sample_data("HIPPOC", 
                        "ATAC", 
                        exclude_outliers = TRUE, 
                        normalized = FALSE, 
                        scratchdir = "/tmp")
```

We can plot the normalized sample-level data for a single feature using
`plot_feature_normalized_data()`. All of the following examples are
different ways to plot the same feature.

```{r plot_feature_normalized_data}
plot_feature_normalized_data(feature = "PROT;SKM-GN;NP_786937.1",
                             add_gene_symbol = TRUE)

# Two different ways to plot the same results for one feature
plot_feature_normalized_data(feature = "PROT;SKM-GN;NP_786937.1",
                             exclude_outliers = TRUE,
                             scale_x_by_time = FALSE)

plot_feature_normalized_data(assay = "PROT",
                             tissue = "SKM-GN",
                             feature_ID = "NP_786937.1",
                             exclude_outliers = TRUE,
                             scale_x_by_time = FALSE)

plot_feature_normalized_data(assay = "PROT",
                             tissue = "SKM-GN",
                             feature_ID = "NP_786937.1",
                             facet_by_sex = TRUE)
```

`combine_normalized_data()` is a wrapper for `load_sample_data()` that
returns combined sample-level normalized data for multiple datasets.
Note that the sample-specific vial labels used as column names for most
of the sample-level data are replaced with rat-specific participant IDs
(PIDs) to allow measurements from multiple datasets for the same animal
to be concatenated.

```{r combine normalized data, eval=F}
# Return all normalized RNA-seq data
data_tra <- combine_normalized_data(assays = "TRNSCRPT")

# Return all normalized proteomics data. Exclude outliers 
data_prot <- combine_normalized_data(assays = c("PROT","UBIQ","PHOSPHO","ACETYL"),
                               exclude_outliers = TRUE)

ratid_tra <- setdiff(names(data_tra), c("feature", "feature_ID", "tissue", "assay"))


ratid_prot <- setdiff(names(data_prot), c("feature", "feature_ID", "tissue", "assay"))


# Return normalized ATAC-seq data for training-regulated features 
data <- combine_normalized_data(assays = "ATAC", 
                               training_regulated_only = TRUE)

# Return all non-epigenetic data
# Note that the "include_epigen" argument is FALSE by default
data <- combine_normalized_data()
```

Similarly, `combine_da_results()` concatenates differential analysis
results for multiple datasets.

```{r combine da results, eval=F}
# Return all non-epigenetic differential analysis results, 
# including meta-regression results for metabolomics
res <- combine_da_results()

# Return METHYL and ATAC differential analysis results for gastrocnemius 
# res = combine_da_results(tissues="SKM-GN", 
#                          assays=c("ATAC","METHYL"),
#                          include_epigen=TRUE)
```

Finally, `list_available_data()` returns a list of all of the available
data objects in the specified package.

```{r list available data}
list_available_data("MotrpacRatTraining6moData")
```

If the
[MotrpacRatTraining6moData](https://motrpac.github.io/MotrpacRatTraining6moData/)
library is attached, you can learn more about any of these data objects
with `?`, e.g.,
[`?TISSUE_COLORS`](https://motrpac.github.io/MotrpacRatTraining6moData/reference/TISSUE_COLORS.html),
and you can load data objects into your environment using `data()`,
e.g., `data(TRAINING_REGULATED_FEATURES)`.

# ANALYSIS

## Principal component analysis

Function to perform a PCA on any of the available NORM datasets

```{r}
# Function to perform PCA and plot the first two PCs
library(impute)  # For KNN imputation

do_pca_plot <- function(norm_data_object, object_name = NULL, impute = FALSE, remove_outliers = FALSE) {
  # Check if the object name is provided or extract it
  if (is.null(object_name)) {
    object_name <- deparse(substitute(norm_data_object))
  }
  
  # Check if the object name contains 'NORM'
  if (!grepl("NORM", object_name)) {
    stop("The provided object does not contain 'NORM' in its name. Please provide a valid NORM_DATA object.")
  }
  
  # Extract the sample IDs from the list
  sample_ids <- setdiff(names(norm_data_object), c("feature", "feature_ID", "tissue", "assay"))
  
  # Subset the data using sample_ids
  data_subset <- norm_data_object[, sample_ids]
  
  # Remove columns with more than 80% missing values
  valid_cols <- colMeans(is.na(data_subset)) < 0.8
  data_subset <- data_subset[, valid_cols]
  
  # Ensure there are enough columns left for PCA
  if (ncol(data_subset) < 2) {
    stop("Not enough columns remaining after filtering for PCA.")
  }
  
  # Impute or clean data based on 'impute' argument
  if (impute) {
    # Impute missing values using KNN, suppressing messages
    sink(file = nullfile())
    data_imputed <- impute::impute.knn(as.matrix(data_subset))$data
    sink()
  } else {
    # Remove rows with NA or Inf values
    data_imputed <- data_subset[stats::complete.cases(data_subset) & apply(data_subset, 1, function(x) all(is.finite(x))), ]
  }
  
  # Ensure there are enough rows left for PCA
  if (nrow(data_imputed) < 2) {
    stop("Not enough rows remaining after data cleaning for PCA.")
  }
  
  # Update the sample_ids based on the final data used in PCA
  valid_sample_ids <- colnames(data_imputed)
  
  # Check if sample_ids match 'pid' or 'viallabel' in PHENO and use them accordingly
  if (all(valid_sample_ids %in% PHENO$pid)) {
    df_sample_ids <- match(valid_sample_ids, PHENO$pid)
  } else if (all(valid_sample_ids %in% PHENO$viallabel)) {
    df_sample_ids <- match(valid_sample_ids, PHENO$viallabel)
  } else {
    stop("Sample IDs do not match either 'pid' or 'viallabel' in PHENO.")
  }
  
  # Perform PCA on the transposed imputed/cleaned data
  pca_result <- stats::prcomp(t(data_imputed))
  
  # Calculate the percentage of variance explained by each PC
  variance_explained <- summary(pca_result)$importance[2, ] * 100
  
  # Create a data frame with phenotypic data and the first 3 PCs
  df <- data.frame(
    group = PHENO[df_sample_ids, "group"],
    sex = PHENO[df_sample_ids, "sex"],
    pca_result$x[, 1:3] # take the first three principal components
  )
  
  # Remove extreme outliers based on IQR for PC1 and PC2 if 'remove_outliers' is TRUE
  if (remove_outliers) {
    for (pc in c("PC1", "PC2")) {
      Q1 <- quantile(df[[pc]], 0.25, na.rm = TRUE)
      Q3 <- quantile(df[[pc]], 0.75, na.rm = TRUE)
      IQR_value <- Q3 - Q1
      
      # Define the boundaries
      lower_bound <- Q1 - 3 * IQR_value
      upper_bound <- Q3 + 3 * IQR_value
      
      # Filter the data to remove outliers
      df <- df[df[[pc]] >= lower_bound & df[[pc]] <= upper_bound, ]
    }
  }
  
  # Ensure that there are still data points left to plot
  if (nrow(df) < 1) {
    stop("No data points left to plot after removing outliers.")
  }
  
  # Plot the first two PCs with variance explained in axis labels
  p <- ggplot2::ggplot(df, ggplot2::aes(x = PC1, y = PC2, fill = group, shape = sex)) + 
    ggplot2::geom_point(size = 3, colour = "black") +
    ggplot2::scale_fill_manual(values = GROUP_COLORS) +
    ggplot2::scale_shape_manual(values = c(male = 21, female = 24)) +
    ggplot2::theme_bw() +
    ggplot2::guides(fill = ggplot2::guide_legend(override.aes = list(shape = 21))) +
    ggplot2::labs(title = paste(object_name),
         x = paste0("PC1 (", round(variance_explained[1], 2), "% Variance)"),
         y = paste0("PC2 (", round(variance_explained[2], 2), "% Variance)")) +
    ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
  
  return(p)
}

```

### One tissue at a time

```{r}
# Example usage of the function
do_pca_plot(ATAC_LIVER_NORM_DATA_05FDR)
do_pca_plot(ATAC_LIVER_NORM_DATA_05FDR, remove_outliers = TRUE)

do_pca_plot(TRNSCRPT_LIVER_NORM_DATA)
do_pca_plot(TRNSCRPT_LIVER_NORM_DATA, remove_outliers = TRUE)

do_pca_plot(PROT_LIVER_NORM_DATA)
do_pca_plot(PHOSPHO_LIVER_NORM_DATA)
```


### All tissues / Assays (except Metabolomics)

Perform PCA and plot for all objects in the list of norm objects
generated above

```{r}
for (object_name in norm_objects) {
  norm_data_object <- get(object_name)
  p <- do_pca_plot(norm_data_object, object_name, impute = FALSE)
  print(p)
}
```

### Metabolomics

The `NORM` data object for metabolomics contains all tissues and assays.
To perform a PCA, it requires previous filtering and assay selection.
Let's massage the data to get a PCA for every assay available for
`LIVER`. Some of the metabolomics assays require imputation due to a 
larger number of missing values.

```{r}
METAB_NORM_ONE_TISSUE <- METAB_NORM_DATA_FLAT %>% filter(tissue == "HEART")

all_assays_one_tissue <- unique(METAB_NORM_ONE_TISSUE$dataset)

for(current_assay in all_assays_one_tissue) {
  METAB_NORM_ONE_ASSAY <- METAB_NORM_ONE_TISSUE %>% filter(dataset == current_assay)
  METAB_NORM_ONE_ASSAY <- METAB_NORM_ONE_ASSAY %>% select(-dataset)
  p <- do_pca_plot(norm_data_object = METAB_NORM_ONE_ASSAY, 
                   object_name = paste0(toupper(current_assay),"_NORM"),
                   impute = TRUE)
  print(p)
}
```

## Differential analysis

More details about the differential analysis methods are available in
the [supplementary information of our *Nature*
publication](https://www.nature.com/articles/s41586-023-06877-w#Sec14).

Simply put: 

- the *training* differential analysis considers all
training groups for each sex (sedentary controls and 4 training time
points) to determine if the analyte significantly changes in either sex
at any point during the training time course. The adjusted p-values from
this analysis were used to determine the set of analytes that are
regulated by endurance exercise training at 5% FDR, referred to as the
*training-regulated features*. 
- The *timewise* differential analysis
performs pairwise contrasts between trained animals at each time point
(1, 2, 4, or 8 weeks) and the sex-matched sedentary control animals.
This gives us sex- and time- specific p-values and effect sizes,
referred to as the *timewise summary statistics*.

### Volcano plots

```{r, fig.width=8, fig.height=10}
# Function to generate a volcano plot for a given data object
create_volcano_plot <- function(data_object) {
  # Extract the name of the data object
  object_name <- deparse(substitute(data_object))
  
  # Filter out outliers with absolute logFC greater than 10
  filtered_data <- data_object %>%
    filter(abs(logFC) <= 10)
  
  p <- ggplot(filtered_data, aes(x = logFC, y = -log10(adj_p_value), color = adj_p_value < 0.05)) +
    geom_point(alpha = 0.6, size = 1.5) +
    scale_color_manual(values = c("TRUE" = "red", "FALSE" = "grey")) +
    facet_grid(comparison_group~sex) +
    theme_linedraw() +
    labs(
      title = paste(object_name),
      x = "logFC",
      y = "-log10(Adjusted P-value)",
      color = "Significant (FDR < 0.05)"
    ) +
    theme(
      plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
      axis.title = element_text(size = 14),
      axis.text = element_text(size = 12),
      strip.text = element_text(size = 12),
      legend.position = "bottom"
    )
  
  return(p)
}

# Example usage
create_volcano_plot(TRNSCRPT_BAT_DA)
```

### Replicating differential analysis results

We provide the following functions to replicate our *timewise* and
*training* differential analysis results for each dataset. Look at the
corresponding documentation for each function for more details, e.g.
`?transcript_timewise_da`.

-   Proteomics (ACETYL, PHOSPHO, PROT, UBIQ)
    -   `proteomics_timewise_da()`\
    -   `proteomics_training_da()`\
-   ATAC
    -   `atac_training_da()`\
    -   `atac_timewise_da()`\
-   IMMUNO
    -   `immuno_timewise_da()`\
    -   `immuno_training_da()`\
-   METAB
    -   `metab_timewise_da()`\
    -   `metab_training_da()`\
    -   `metab_meta_regression()`\
-   METHYL
    -   `rrbs_differential_analysis()`\
-   TRNSCRPT
    -   `transcript_timewise_da()`\
    -   `transcript_training_da()`

Here, we replicate the protein acetylation timewise differential
analysis results provided in
[MotrpacRatTraining6moData](https://motrpac.github.io/MotrpacRatTraining6moData/).

```{r replicate prot}
timewiseList = list()
for (tissue in c("HEART","LIVER")){
  timewiseList[[tissue]] = proteomics_timewise_da("ACETYL", tissue)
}
timewise = do.call("rbind", timewiseList)

# merge with version of results in MotrpacRatTraining6moData
original = combine_da_results(assays="ACETYL")
merged = merge(original, timewise, 
               by=c("feature_ID","assay","assay_code","tissue","tissue_code","sex","comparison_group"),
               suffixes=c("_orginal","_reproduced"))

# plot 
plot(-log10(merged$p_value_orginal), 
     -log10(merged$p_value_reproduced), 
     type="p", cex=0.5,
     xlab="Original timewise p-value (-log10)",
     ylab="Reproduced timewise p-value (-log10)",
     main="Timewise DA results for ACETYL data")
```

### Plot DA for a single feature

We can plot the differential analysis results for a single feature using
`plot_feature_logfc()`. All of the following examples are different ways
to plot results for the same feature.

```{r}
plot_feature_logfc(feature = "PROT;SKM-GN;NP_786937.1",
                   add_gene_symbol = TRUE)

plot_feature_logfc(assay = "PROT",
                   tissue = "SKM-GN",
                   feature_ID = "NP_786937.1",
                   scale_x_by_time = FALSE)

plot_feature_logfc(assay = "PROT",
                   tissue = "SKM-GN",
                   feature_ID = "NP_786937.1",
                   facet_by_sex = TRUE,
                   add_adj_p = FALSE)
```

## Bayesian graphical clustering

For a detail tutorial on the graphical clustering analysis, please,
visit the original vignette part of the MotrpacRatTraining6mo package
[here](https://motrpac.github.io/MotrpacRatTraining6mo/articles/MotrpacRatTraining6mo.html#bayesian-graphical-clustering).

## Getting help

For questions, bug reporting, and data requests for the
`MotrpacRatTraining6moData` package, please [submit a new issue
here](https://github.com/MoTrPAC/MotrpacRatTraining6moData/issues) and
include as many details as possible.

If the issue is related to functions provided in the
`MotrpacRatTraining6mo` package, please submit an issue
[here](https://github.com/MoTrPAC/MotrpacRatTraining6mo/issues).

## Acknowledgements

MoTrPAC is supported by the National Institutes of Health (NIH) Common
Fund through cooperative agreements managed by the National Institute of
Diabetes and Digestive and Kidney Diseases (NIDDK), National Institute
of Arthritis and Musculoskeletal Diseases (NIAMS), and National
Institute on Aging (NIA).

## Session Info

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