---
title: "Case study for tidyMass"
author:
- name: Hong Yan (https://ysph.yale.edu/profile/hong_yan/)
- name: Xiaotao Shen (https://www.shenxt.info/)
- name: Chuhu Wang (https://www.linkedin.com/in/chuchu-wang-71331190/)
date: "Created on 2022-03-03 and updated on `r Sys.Date()`"
output:
  html_document:
    number_sections: true
    toc_depth: 3
    df_print: paged
    toc: true
    toc_float: true
    smooth_scroll: true
    toc_collapsed: true
  pdf_document:
    toc: no
vignette: >
  %\VignetteIndexEntry{tidymass_case_study}
  %\VignettePackage{tidymass}
  % \VignetteEngine{knitr::rmarkdown}
  % \usepackage[utf8]{inputenc}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE, echo=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE,
  message = TRUE,
  out.width = "100%"
)
```

# Introduction

Data introduction. We have RPLC (positive and negative mode), HILIC (positive and negative mode). More information can be found here.

[Sex Differences in Colon Cancer Metabolism Reveal A Novel Subphenotype](https://www.nature.com/articles/s41598-020-61851-0)

> It may take a while to run all the code.

# Download data

Mass spectrometry raw data (mzXML) for the case study in this paper is accessible on `MetaboLights` with MTBLS1122 (HILIC positive), MTBLS1124 (HILIC negative), MTBLS1122 (RPLC positive) and MTBLS1130 (RPLC negative). The MS2 data (mgf) and processed data (“mass_dataset” class) from the `massProcesser` package are available on the tidyMass project website (https://tidymass.github.io/case_study_data/).

> Note: In some steps, we just skip them to save time, try to change `eval=FALSE` to `eval=TRUE` so you can run this code chunk. For the annotation step, because we can't share the in-house database, so we just provide the `mass_dataset` class after annotation. We apologize for the inconvenience.

# Data preparation

Please download all the data and [code](https://drive.google.com/file/d/1vop-hU1J5Ykw7ta0dOb6eAtQjA0R3iSx/view?usp=sharing), and put them in one folder named as `case_study`.

![](figures/fig1.png)

# Install packages

Please install the pacakge we need in this analysis.

## `tidymass`

Please refer [this document](https://tidymass.github.io/tidymass/articles/tidymass_install.html).

```{r,eval=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE, results = 'hide'}
if(!require(remotes)){
  install.packages("remotes")
}

if(!require(remotes)){
  remotes::install_gitlab("jaspershen/tidymass")
}
```

## `tidyverse`

```{r,eval=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE, results = 'hide'}
if(!require(tidyverse)){
  install.packages("tidyverse")
}
```

## Other packages

```{r,eval=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE, results = 'hide'}
if(!require(BiocManager)){
  install.packages("BiocManager")
}

if(!require(ComplexHeatmap)){
  BiocManager::install("ComplexHeatmap")
}

if(!require(ggraph)){
  install.packages("ggraph")
}

if(!require(tidygraph)){
  install.packages("tidygraph")
}

if(!require(extrafont)){
  install.packages("extrafont")
}

if(!require(shadowtext)){
  install.packages("shadowtext")
}
```

# Raw data processing

`massprocesser` package is used to do the raw data processing. Please refer this [website](file:///Users/xiaotaoshen/tidymass/massprocesser/docs/index.html) for more information.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
library(tidymass)
library(tidyverse)
```

## RPLC positive mode

The code used to do raw data processing.

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
process_data(
  path = "mzxml_ms1_data/RPLC/POS/",
  polarity = "positive",
  ppm = 20,
  peakwidth = c(5, 30),
  threads = 6,
  output_tic = FALSE,
  output_bpc = FALSE,
  output_rt_correction_plot = FALSE,
  min_fraction = 0.5,
  group_for_figure = "QC"
)
```

All the results will be placed in the folder named `case_study/data/mzxml_ms1_data/POS/Result`. More information about that can be found [here](https://tidymass.github.io/massprocesser/articles/raw_data_processing.html).

You can just load the `object`, which is a `mass_dataset` class object.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
load("mzxml_ms1_data/RPLC/POS/Result/object")
object
dim(object)
```

We can see that there are 14,585 metabolic features in positive mode.

RT correction plot

load data

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
load("mzxml_ms1_data/RPLC/POS/Result/intermediate_data/xdata2")
```

Set the group_for_figure if you want to show specific groups. And set it as "all" if you want to show all samples.

We can use the `plot_adjusted_rt()` function to get the interactive plot.

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
####set the group_for_figure if you want to show specific groups. 
####And set it as "all" if you want to show all samples.
massprocesser::plot_adjusted_rt(object = xdata2, 
                 group_for_figure = "QC", 
                 interactive = TRUE)
```

## RPLC negative mode

The processing of negative mode is same with positive mode data.

Same with positive mode, change `polarity` to `negative`.

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
process_data(
  path = "mzxml_ms1_data/RPLC/NEG/",
  polarity = "negative",
  ppm = 20,
  peakwidth = c(5, 30),
  threads = 6,
  output_tic = FALSE,
  output_bpc = FALSE,
  output_rt_correction_plot = FALSE,
  min_fraction = 0.5,
  group_for_figure = "QC"
)
```

Same with positive mode.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
load("mzxml_ms1_data/RPLC/NEG/Result/object")
object
dim(object)
```

We can see that there are 4,685 metabolic features in negative mode.

## HILIC positive mode

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
process_data(
  path = "mzxml_ms1_data/HILIC/POS/",
  polarity = "positive",
  ppm = 20,
  peakwidth = c(5, 30),
  threads = 6,
  output_tic = FALSE,
  output_bpc = FALSE,
  output_rt_correction_plot = FALSE,
  min_fraction = 0.5,
  group_for_figure = "QC"
)
```

## HILIC negative mode

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
process_data(
  path = "mzxml_ms1_data/HILIC/NEG/",
  polarity = "negative",
  ppm = 20,
  peakwidth = c(5, 30),
  threads = 6,
  output_tic = FALSE,
  output_bpc = FALSE,
  output_rt_correction_plot = FALSE,
  min_fraction = 0.5,
  group_for_figure = "QC"
)
```

# Explore data

After the [`raw data processing`](https://tidymass.github.io/tidymass/articles/raw_data_processing.html), peak tables for positive and negative mode will be generated.

## RPLC positive mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
load("mzxml_ms1_data/RPLC/POS/Result/object")
object_rplc_pos <- object
object_rplc_pos
```

Read sample information.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
sample_info_pos <- 
  readr::read_csv("sample_info/RPLC/sample_info_pos_1.csv")
head(sample_info_pos)
```

Add `sample_info_pos` to `object_rplc_pos`

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_pos %>% 
  extract_sample_info() %>% 
  head()

object_rplc_pos <- 
  object_rplc_pos %>% 
  activate_mass_dataset(what = "sample_info") %>%
  dplyr::select(-c("group", "class", "injection.order"))
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_pos <-
  object_rplc_pos %>%
  activate_mass_dataset(what = "sample_info") %>%
  left_join(sample_info_pos,
            by = "sample_id")

object_rplc_pos %>% 
  extract_sample_info() %>% 
  head()
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
dim(object_rplc_pos)
object_rplc_pos %>% 
  activate_mass_dataset(what = "sample_info") %>% 
  dplyr::count(class)

object_rplc_pos %>% 
  activate_mass_dataset(what = "sample_info") %>% 
  dplyr::count(group)

object_rplc_pos %>% 
  activate_mass_dataset(what = "sample_info") %>% 
  dplyr::count(batch)
```

Next, we can get the peak distributation plot of positive mode.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_pos %>%
  `+`(1) %>% 
  log(10) %>%
  show_mz_rt_plot() +
  scale_size_continuous(range = c(0.01, 2))
```

We can explore the missing values in positive mode data.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
get_mv_number(object = object_rplc_pos)
```

Missing value number in each sample.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
get_mv_number(object = object_rplc_pos, by = "variable") %>% 
  head()
```

Missing value number in each variable.

We can use the figure to show the missing value information.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
library(grid)
show_missing_values(object = object_rplc_pos, 
                    show_column_names = FALSE, 
                    percentage = TRUE)
```

Show the mvs in samples.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
show_sample_missing_values(object = object_rplc_pos, 
                           percentage = TRUE) +
  scale_size_continuous(range = c(0.1, 1)) +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank())
```

Show the mvs in variables

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
show_variable_missing_values(object = object_rplc_pos, 
                             percentage = TRUE, 
                             show_x_text = FALSE, 
                             show_x_ticks = FALSE) +
  scale_size_continuous(range = c(0.01, 1))
```

## RPLC Negative mode

Load `object`.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
load("mzxml_ms1_data/RPLC/NEG/Result/object")
object_rplc_neg <- object
object_rplc_neg
```

Read sample information.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
sample_info_neg <- readr::read_csv("sample_info/RPLC/sample_info_neg_1.csv")
head(sample_info_neg)
```

Add `sample_info_neg` to `object_rplc_neg`

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_neg <- 
  object_rplc_neg %>% 
  activate_mass_dataset(what = "sample_info") %>%
  dplyr::select(-c("group", "class", "injection.order"))
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_neg <-
  object_rplc_neg %>%
  activate_mass_dataset(what = "sample_info") %>%
  left_join(sample_info_neg,
            by = "sample_id")
```

## HILIC positive mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
load("mzxml_ms1_data/HILIC/POS/Result/object")
object_hilic_pos <- object
object_hilic_pos
```

Read sample information.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
sample_info_pos <- readr::read_csv("sample_info/HILIC/sample_info_pos_1.csv")
head(sample_info_pos)
```

Add `sample_info_pos` to `object_hilic_pos`

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_pos <- 
  object_hilic_pos %>% 
  activate_mass_dataset(what = "sample_info") %>%
  dplyr::select(-c("group", "class", "injection.order"))
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_pos <-
  object_hilic_pos %>%
  activate_mass_dataset(what = "sample_info") %>%
  left_join(sample_info_pos,
            by = "sample_id")
```

## HILIC Negative mode

Load `object`.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
load("mzxml_ms1_data/HILIC/NEG/Result/object")
object_hilic_neg <- object
object_hilic_neg
```

Read sample information.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
sample_info_neg <- readr::read_csv("sample_info/HILIC/sample_info_neg_1.csv")
head(sample_info_neg)
```

Add `sample_info_neg` to `object_hilic_neg`

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_neg <- 
  object_hilic_neg %>% 
  activate_mass_dataset(what = "sample_info") %>%
  dplyr::select(-c("group", "class", "injection.order"))
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_neg <-
  object_hilic_neg %>%
  activate_mass_dataset(what = "sample_info") %>%
  left_join(sample_info_neg,
            by = "sample_id")
```

# Data cleaning

## Data quality assessment before data cleaning

Change batch to character.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_pos <- 
  object_rplc_pos %>% 
  activate_mass_dataset(what = "sample_info") %>% 
  dplyr::mutate(batch = as.character(batch))

object_rplc_neg <- 
  object_rplc_neg %>% 
  activate_mass_dataset(what = "sample_info") %>% 
  dplyr::mutate(batch = as.character(batch))

object_hilic_pos <- 
  object_hilic_pos %>% 
  activate_mass_dataset(what = "sample_info") %>% 
  dplyr::mutate(batch = as.character(batch))

object_hilic_neg <- 
  object_hilic_neg %>% 
  activate_mass_dataset(what = "sample_info") %>% 
  dplyr::mutate(batch = as.character(batch))
```

Save them for subsequent analysis.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
save(object_rplc_pos, file = "data_cleaning/RPLC/POS/object_rplc_pos")
save(object_rplc_neg, file = "data_cleaning/RPLC/NEG/object_rplc_neg")

save(object_hilic_pos, file = "data_cleaning/HILIC/POS/object_hilic_pos")
save(object_hilic_neg, file = "data_cleaning/HILIC/NEG/object_hilic_neg")
```


Here, we can use the `massqc` package to [assess the data quality](https://tidymass.github.io/massqc/).

We can just use the `massqc_report()` function to get a html format quality assessment report.

### RPLC Positive mode

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
massqc::massqc_report(object = object_rplc_pos, 
                      path = "data_cleaning/RPLC/POS/data_quality_before_data_cleaning")
```

A html format report and pdf figures will be placed in the folder `data/data_cleaning/POS/data_quality_before_data_cleaning/Report`.

The html report is below:

![](figures/fig2.png)

### RPLC negative mode.

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
massqc::massqc_report(object = object_rplc_neg, 
                      path = "data_cleaning/RPLC/NEG/data_quality_before_data_cleaning")
```

### HILIC positive mode

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
massqc::massqc_report(object = object_hilic_pos, 
                      path = "data_cleaning/HILIC/POS/data_quality_before_data_cleaning")
```

### HILIC negative mode

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
massqc::massqc_report(object = object_hilic_neg, 
                      path = "data_cleaning/HILIC/NEG/data_quality_before_data_cleaning")
```

## Remove noisy metabolic features

Remove variables which have MVs in more than 20% QC samples and in at lest 50% samples in control group or case group.

### RPLC positive mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_pos %>% 
  activate_mass_dataset(what = "sample_info") %>% 
  dplyr::count(group)
```

MV percentage in QC samples.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
show_variable_missing_values(object = object_rplc_pos %>% 
                               activate_mass_dataset(what = "sample_info") %>% 
                               filter(class == "QC"), 
                             percentage = TRUE) +
  scale_size_continuous(range = c(0.01, 2))
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
qc_id <-
  object_rplc_pos %>%
  activate_mass_dataset(what = "sample_info") %>%
  filter(class == "QC") %>%
  pull(sample_id)

control_id <-
 object_rplc_pos %>%
  activate_mass_dataset(what = "sample_info") %>%
  filter(stringr::str_detect(group, "control")) %>%
  pull(sample_id)

case_id <-
  object_rplc_pos %>%
  activate_mass_dataset(what = "sample_info") %>%
  filter(stringr::str_detect(group, "LCC|RCC")) %>%
  pull(sample_id)

object_rplc_pos <-
  object_rplc_pos %>%
  mutate_variable_na_freq(according_to_samples = qc_id) %>%
  mutate_variable_na_freq(according_to_samples = control_id) %>%
  mutate_variable_na_freq(according_to_samples = case_id)

head(extract_variable_info(object_rplc_pos))
```

Remove variables.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_pos <-
  object_rplc_pos %>%
  activate_mass_dataset(what = "variable_info") %>%
  filter(na_freq < 0.2 &
           (na_freq.1 < 0.5 |
              na_freq.2 < 0.5))
object_rplc_pos
```

Only 6,132 variables left.

### RPLC negative mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
qc_id <-
  object_rplc_neg %>%
  activate_mass_dataset(what = "sample_info") %>%
  filter(class == "QC") %>%
  pull(sample_id)

control_id <-
 object_rplc_neg %>%
  activate_mass_dataset(what = "sample_info") %>%
  filter(stringr::str_detect(group, "control")) %>%
  pull(sample_id)

case_id <-
  object_rplc_neg %>%
  activate_mass_dataset(what = "sample_info") %>%
  filter(stringr::str_detect(group, "LCC|RCC")) %>%
  pull(sample_id)

object_rplc_neg =
  object_rplc_neg %>%
  mutate_variable_na_freq(according_to_samples = qc_id) %>%
  mutate_variable_na_freq(according_to_samples = control_id) %>%
  mutate_variable_na_freq(according_to_samples = case_id)
```

Remove variables.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_neg <-
  object_rplc_neg %>%
  activate_mass_dataset(what = "variable_info") %>%
  filter(
    na_freq < 0.2 &
      (
        na_freq.1 < 0.5 |
          na_freq.2 < 0.5
      )
  )
object_rplc_neg
```

2,719 features left.

### HILIC positive mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
qc_id <-
  object_hilic_pos %>%
  activate_mass_dataset(what = "sample_info") %>%
  filter(class == "QC") %>%
  pull(sample_id)

control_id <-
 object_hilic_pos %>%
  activate_mass_dataset(what = "sample_info") %>%
  filter(stringr::str_detect(group, "control")) %>%
  pull(sample_id)

case_id <-
  object_hilic_pos %>%
  activate_mass_dataset(what = "sample_info") %>%
  filter(stringr::str_detect(group, "LCC|RCC")) %>%
  pull(sample_id)

object_hilic_pos =
  object_hilic_pos %>%
  mutate_variable_na_freq(according_to_samples = qc_id) %>%
  mutate_variable_na_freq(according_to_samples = control_id) %>%
  mutate_variable_na_freq(according_to_samples = case_id)
```

Remove variables.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_pos <-
  object_hilic_pos %>%
  activate_mass_dataset(what = "variable_info") %>%
  filter(
    na_freq < 0.2 &
      (
        na_freq.1 < 0.5 |
          na_freq.2 < 0.5
      )
  )
object_rplc_neg
```

### HILIC negative mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
qc_id <-
  object_hilic_neg %>%
  activate_mass_dataset(what = "sample_info") %>%
  filter(class == "QC") %>%
  pull(sample_id)

control_id <-
 object_hilic_neg %>%
  activate_mass_dataset(what = "sample_info") %>%
  filter(stringr::str_detect(group, "control")) %>%
  pull(sample_id)

case_id <-
  object_hilic_neg %>%
  activate_mass_dataset(what = "sample_info") %>%
  filter(stringr::str_detect(group, "LCC|RCC")) %>%
  pull(sample_id)

object_hilic_neg =
  object_hilic_neg %>%
  mutate_variable_na_freq(according_to_samples = qc_id) %>%
  mutate_variable_na_freq(according_to_samples = control_id) %>%
  mutate_variable_na_freq(according_to_samples = case_id)
```

Remove variables.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_neg <-
  object_hilic_neg %>%
  activate_mass_dataset(what = "variable_info") %>%
  filter(
    na_freq < 0.2 &
      (
        na_freq.1 < 0.5 |
          na_freq.2 < 0.5
      )
  )
object_hilic_neg
```

## Filter outlier samples

We can use the `detect_outlier()` to find the outlier samples.

More information about how to detect outlier samples can be found [here](https://privefl.github.io/blog/detecting-outlier-samples-in-pca/).

### RPLC positive mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
##change class to subject
object_rplc_pos <-
  object_rplc_pos %>%
  activate_mass_dataset(what = "sample_info") %>%
  dplyr::mutate(class = case_when(class == "QC" ~ class,
                                  TRUE ~ "Subject"))
```

Detect outlier samples.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
outlier_samples <-
  object_rplc_pos %>%
  `+`(1) %>% 
  log(2) %>%
  scale() %>%
  detect_outlier(na_percentage_cutoff = 0.5)

outlier_samples
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
outlier_table <-
extract_outlier_table(outlier_samples)

outlier_table %>% 
  apply(1, function(x){
    sum(x)
  }) %>% 
  `>`(0) %>% 
  which()
```

Here, we don't remove outlier samples.

### RPLC negative mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
##change class to subject
object_rplc_neg <-
  object_rplc_neg %>%
  activate_mass_dataset(what = "sample_info") %>%
  dplyr::mutate(class = case_when(class == "QC" ~ class,
                                  TRUE ~ "Subject"))
```

Detect outlier samples.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
outlier_samples <-
  object_rplc_neg %>%
  `+`(1) %>% 
  log(2) %>%
  scale() %>%
  detect_outlier()

outlier_samples
```

### HILIC positive mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
##change class to subject
object_hilic_pos <-
  object_hilic_pos %>%
  activate_mass_dataset(what = "sample_info") %>%
  dplyr::mutate(class = case_when(class == "QC" ~ class,
                                  TRUE ~ "Subject"))
```

Detect outlier samples.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
outlier_samples <-
  object_hilic_pos %>%
  `+`(1) %>% 
  log(2) %>%
  scale() %>%
  detect_outlier()

outlier_samples
```

### HILIC negative mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
##change class to subject
object_hilic_neg <-
  object_hilic_neg %>%
  activate_mass_dataset(what = "sample_info") %>%
  dplyr::mutate(class = case_when(class == "QC" ~ class,
                                  TRUE ~ "Subject"))
```

Detect outlier samples.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
outlier_samples <-
  object_hilic_neg %>%
  `+`(1) %>% 
  log(2) %>%
  scale() %>%
  detect_outlier()

outlier_samples
```

## Missing value imputation

```{r,eval=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE, results = 'hide'}
object_rplc_pos <- 
  impute_mv(object = object_rplc_pos, 
            method = "knn", colmax = 0.9)
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE, results = 'hide'}
object_rplc_neg <- 
  impute_mv(object = object_rplc_neg, method = "knn")
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE, results = 'hide'}
object_hilic_pos <- 
  impute_mv(object = object_hilic_pos, method = "knn")
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE, results = 'hide'}
object_hilic_neg <- 
  impute_mv(object = object_hilic_neg, method = "knn")
```

## Data normalization and integration

### RPLC positive mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_pos <- 
  normalize_data(object_rplc_pos, method = "svr")

object_rplc_pos2 <- 
  integrate_data(object_rplc_pos, method = "subject_median")
```

Save it for subsequent analysis.

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
save(object_rplc_pos2, file = "data_cleaning/RPLC/POS/object_rplc_pos2")
```

### RPLC negative mode


```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_neg <- 
  normalize_data(object_rplc_neg, method = "svr")

object_rplc_neg2 <- 
  integrate_data(object_rplc_neg, method = "subject_median")
```

Save it for subsequent analysis.

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
save(object_rplc_neg2, file = "data_cleaning/RPLC/NEG/object_rplc_neg2")
```

### HILIC positive mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_pos <- 
  normalize_data(object_hilic_pos, method = "svr")

object_hilic_pos2 <- 
  integrate_data(object_hilic_pos, method = "subject_median")
```

Save it for subsequent analysis.

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
save(object_hilic_pos2, file = "data_cleaning/HILIC/POS/object_hilic_pos2")
```

### HILIC negative mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_neg <- 
  normalize_data(object_hilic_neg, method = "svr")

object_hilic_neg2 <- 
  integrate_data(object_hilic_neg, method = "subject_median")
```

Save it for subsequent analysis.

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
save(object_hilic_neg2, file = "data_cleaning/HILIC/NEG/object_hilic_neg2")
```

## Data quality assessment after data cleaning

We can use the `massqc_report()` function to get a html format quality assessment report.

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
massqc::massqc_report(object = object_rplc_pos2, 
                      path = "data_cleaning/RPLC/POS/data_quality_after_data_cleaning")
```

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
massqc::massqc_report(object = object_rplc_neg2, 
                      path = "data_cleaning/RPLC/NEG/data_quality_after_data_cleaning")
```


```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
massqc::massqc_report(object = object_hilic_pos2, 
                      path = "data_cleaning/HILIC/POS/data_quality_after_data_cleaning")
```

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
massqc::massqc_report(object = object_hilic_neg2, 
                      path = "data_cleaning/HILIC/NEG/data_quality_after_data_cleaning")
```

# Metabolite annotation

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
library(tidymass)
library(tidyverse)
```

## Add MS2 spectra to object

Download the [MS2 data here](https://drive.google.com/file/d/1fR75tWp9jqp8vuWq85GvDGJ864Pb6AH2/view?usp=sharing).

Place it in the `data` folder and uncompress it.

### RPLC positive mode

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_pos2 <-
   mutate_ms2(
    object = object_rplc_pos2,
    column = "rplc",
    polarity = "positive",
    ms1.ms2.match.mz.tol = 5,
    ms1.ms2.match.rt.tol = 30,
    path = "mgf_ms2_data/RPLC/POS/"
  )
```

### RPLC negative mode

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_neg2 <-
   mutate_ms2(
    object = object_rplc_neg2,
    column = "rplc",
    polarity = "negative",
    ms1.ms2.match.mz.tol = 5,
    ms1.ms2.match.rt.tol = 30,
    path = "mgf_ms2_data/RPLC/NEG/"
  )
```

### HILIC positive mode

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_pos2 <-
   mutate_ms2(
    object = object_hilic_pos2,
    column = "hilic",
    polarity = "positive",
    ms1.ms2.match.mz.tol = 5,
    ms1.ms2.match.rt.tol = 30,
    path = "mgf_ms2_data/HILIC/POS/"
  )
```

### HILIC negative mode

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_neg2 <-
   mutate_ms2(
    object = object_hilic_neg2,
    column = "hilic",
    polarity = "negative",
    ms1.ms2.match.mz.tol = 5,
    ms1.ms2.match.rt.tol = 30,
    path = "mgf_ms2_data/HILIC/NEG/"
  )
```

## Metabolite annotation using `metid` pacakge.

Metabolite annotation is based on the [`metid` package](https://tidymass.github.io/metid/).

### Download database

We need to download MS1 and MS2 database from `metid` [website](https://tidymass.github.io/metid/articles/public_databases.html).

The in-house database is constructed with MS1 and RT information.

We also download the public databases, `snyder_database_hilic0.0.3`. And place them in a new folder named as `metabolite_annotation`.

### RPLC positive mode

#### Annotate features using in-house database.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
load("metabolite_annotation/database/RPLC.database")
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_pos3 <- 
  annotate_metabolites_mass_dataset(object = object_rplc_pos2, 
                                    ms1.match.ppm = 15, 
                                    rt.match.tol = 30, 
                                    polarity = "positive",
                                    database = RPLC.database)
```
#### Annotate features using `snyder_database_rplc0.0.3`.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_pos3 <- 
  annotate_metabolites_mass_dataset(object = object_rplc_pos3, 
                                    ms1.match.ppm = 15, 
                                    rt.match.tol = 1000000,
                                    polarity = "positive",
                                    database = snyder_database_rplc0.0.3)
```

Save it for subsequent analysis.

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
save(object_rplc_pos3, file = "metabolite_annotation/RPLC/POS/object_rplc_pos3")
```

### RPLC negative mode

#### Annotate features using in-house database.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_neg3 <- 
  annotate_metabolites_mass_dataset(object = object_rplc_neg2, 
                                    ms1.match.ppm = 15, 
                                    rt.match.tol = 30, 
                                    polarity = "negative",
                                    database = RPLC.database)
```
#### Annotate features using `snyder_database_rplc0.0.3`.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_neg3 <- 
  annotate_metabolites_mass_dataset(object = object_rplc_neg3, 
                                    ms1.match.ppm = 15, 
                                    rt.match.tol = 1000000,
                                    polarity = "negative",
                                    database = snyder_database_rplc0.0.3)
```

Save it for subsequent analysis.

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
save(object_rplc_neg3, file = "metabolite_annotation/RPLC/NEG/object_rplc_neg3")
```

### HILIC positive mode

#### Annotate features using in-house database.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
load("metabolite_annotation/database/HILIC.database")
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_pos3 <- 
  annotate_metabolites_mass_dataset(object = object_hilic_pos2, 
                                    ms1.match.ppm = 15, 
                                    rt.match.tol = 30, 
                                    polarity = "positive",
                                    database = HILIC.database)
```
#### Annotate features using `snyder_database_hilic0.0.3`.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
load("metabolite_annotation/database/snyder_database_hilic0.0.3")
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_pos3 <- 
  annotate_metabolites_mass_dataset(object = object_hilic_pos3, 
                                    ms1.match.ppm = 15, 
                                    rt.match.tol = 1000000,
                                    polarity = "positive",
                                    database = snyder_database_hilic0.0.3)
```

Save it for subsequent analysis.

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
save(object_hilic_pos3, file = "metabolite_annotation/HILIC/POS/object_hilic_pos3")
```

### HILIC negative mode

#### Annotate features using in-house database.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_neg3 <- 
  annotate_metabolites_mass_dataset(object = object_hilic_neg2, 
                                    ms1.match.ppm = 15, 
                                    rt.match.tol = 30, 
                                    polarity = "negative",
                                    database = HILIC.database)
```
#### Annotate features using `snyder_database_hilic0.0.3`.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_neg3 <- 
  annotate_metabolites_mass_dataset(object = object_hilic_neg3, 
                                    ms1.match.ppm = 15, 
                                    rt.match.tol = 1000000,
                                    polarity = "negative",
                                    database = snyder_database_hilic0.0.3)
```

Save it for subsequent analysis.

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
save(object_hilic_neg3, file = "metabolite_annotation/HILIC/NEG/object_hilic_neg3")
```
## Annotation result

The annotation result can be got using the `extract_annotation_table` function.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
head(extract_annotation_table(object = object_rplc_pos3))
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
variable_info_pos <-
  extract_variable_info(object = object_rplc_pos3)
head(variable_info_pos)
```

# Statistical analysis and pathway analysis

## Remove the features without annotations

We only remain the annotations with level 2.

### RPLC positive mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_pos3 <- 
  object_rplc_pos3 %>% 
  activate_mass_dataset(what = "annotation_table") %>% 
  filter(!is.na(Level)) %>% 
  filter(Level == 1 | Level == 2)
```

### RPLC negative mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc_neg3 <- 
  object_rplc_neg3 %>% 
  activate_mass_dataset(what = "annotation_table") %>% 
  filter(!is.na(Level)) %>% 
  filter(Level == 1 | Level == 2)
```

### HILIC positive mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_pos3 <- 
  object_hilic_pos3 %>% 
  activate_mass_dataset(what = "annotation_table") %>% 
  filter(!is.na(Level)) %>% 
  filter(Level == 1 | Level == 2)
```

### RPLC negative mode

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic_neg3 <- 
  object_hilic_neg3 %>% 
  activate_mass_dataset(what = "annotation_table") %>% 
  filter(!is.na(Level)) %>% 
  filter(Level == 1 | Level == 2)
```

## Merge data

We need to merge RPLC, HILIC positive and negative mode data as one `mass_dataset` class.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
head(colnames(object_rplc_pos3))
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
head(colnames(object_rplc_neg3))
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
head(colnames(object_hilic_pos3))
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
head(colnames(object_hilic_neg3))
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_rplc <-
  merge_mass_dataset(
    x = object_rplc_pos3,
    y = object_rplc_neg3,
    sample_direction = "inner",
    variable_direction = "full",
    sample_by = "sample_id",
    variable_by = c("variable_id", "mz", "rt")
  )
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object_hilic <-
  merge_mass_dataset(
    x = object_hilic_pos3,
    y = object_hilic_neg3,
    sample_direction = "inner",
    variable_direction = "full",
    sample_by = "sample_id",
    variable_by = c("variable_id", "mz", "rt")
  )
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object <-
  merge_mass_dataset(
    x = object_rplc,
    y = object_hilic,
    sample_direction = "inner",
    variable_direction = "full",
    sample_by = "sample_id",
    variable_by = c("variable_id", "mz", "rt")
  )
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object
```

## Remove redundant metabolites

Remove the redundant annotated metabolites bases on `Level` and score.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object <-
  object %>% 
  activate_mass_dataset(what = "annotation_table") %>% 
  group_by(Compound.name) %>% 
  filter(Level == min(Level)) %>% 
  filter(Total.score == max(Total.score)) %>% 
  slice_head(n = 1)
```

Save it for subsequent analysis.

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
save(object, file = "statistical_analysis/object")
```

## Trace processing information of object

Then we can use the `report_parameters()` function to trace processing information of object.

All the analysis results will be placed in a folder named as `statistical_analysis`.

```{r,eval=FALSE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
report_parameters(object = object, path = "statistical_analysis")
```

A html format file named as `parameter_report.html` will be generated.

![](figures/fig3.png)

## Clustering of all samples

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
load("statistical_analysis/object")
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE, results = 'hide'}
temp_object <- 
object %>% 
  activate_mass_dataset(what = "sample_info") %>% 
  filter(group != "QC") %>% 
  `+`(1) %>% 
  log(2) %>% 
  scale()
library(ComplexHeatmap)
h1 <- HeatmapAnnotation(class = extract_sample_info(temp_object)$group)
massstat::Heatmap(matrix = temp_object, 
                  name = "z-score", 
                  row_names_gp = gpar(cex = 0.2),
                  column_names_gp = gpar(cex = 0.2),
                  top_annotation = h1)
```

## Differential expression metabolites

We want to see if the dysregulated metabolites and pathways are different in males and female cancer tissues.

### Calculate the fold changes.

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
femaleLCC <-
  object %>%
  activate_mass_dataset(what = "sample_info") %>%
  dplyr::filter(group == "female LCC") %>%
  dplyr::pull(sample_id)

femaleRCC <-
  object %>%
  activate_mass_dataset(what = "sample_info") %>%
  dplyr::filter(group == "female RCC") %>%
  dplyr::pull(sample_id)

maleLCC <-
  object %>%
  activate_mass_dataset(what = "sample_info") %>%
  dplyr::filter(group == "male LCC") %>%
  dplyr::pull(sample_id)

maleRCC <-
  object %>%
  activate_mass_dataset(what = "sample_info") %>%
  dplyr::filter(group == "male RCC") %>%
  dplyr::pull(sample_id)

femalecontrol <-
  object %>%
  activate_mass_dataset(what = "sample_info") %>%
  dplyr::filter(group == "female control") %>%
  dplyr::pull(sample_id)

malecontrol <-
  object %>%
  activate_mass_dataset(what = "sample_info") %>%
  dplyr::filter(group == "male control") %>%
  dplyr::pull(sample_id)
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object <-
  mutate_fc(object = object, 
            control_sample_id = malecontrol, 
            case_sample_id = maleLCC, 
            mean_median = "mean")
object <- 
  object %>% 
  activate_mass_dataset(what = "variable_info") %>% 
  dplyr::rename(fc_malelcc_vs_malecontrol = fc)
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object <-
  mutate_fc(object = object, 
            control_sample_id = femalecontrol, 
            case_sample_id = femaleLCC, 
            mean_median = "mean")
object <- 
  object %>% 
  activate_mass_dataset(what = "variable_info") %>% 
  dplyr::rename(fc_femalelcc_vs_femalecontrol = fc)
```

### Calculate p values

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object <-
  mutate_p_value(
    object = object,
    control_sample_id = malecontrol,
    case_sample_id = maleLCC,
    method = "t.test",
    p_adjust_methods = "BH"
  )
object <- 
  object %>% 
  activate_mass_dataset(what = "variable_info") %>% 
  dplyr::rename(p_value_malelcc_vs_malecontrol = p_value,
         p_value_adjust_malelcc_vs_malecontrol = p_value_adjust)
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
object <-
  mutate_p_value(
    object = object,
    control_sample_id = femalecontrol,
    case_sample_id = femaleLCC,
    method = "t.test",
    p_adjust_methods = "BH"
  )
object <- 
  object %>% 
  activate_mass_dataset(what = "variable_info") %>% 
  dplyr::rename(p_value_femalelcc_vs_femalecontrol = p_value,
         p_value_adjust_femalelcc_vs_femalecontrol = p_value_adjust)
```

### Volcano plot

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
plot <- 
volcano_plot(object = object,
             fc_column_name = "fc_femalelcc_vs_femalecontrol",
             p_value_column_name = "p_value_adjust_femalelcc_vs_femalecontrol",
             add_text = TRUE,
             text_from = "Compound.name", 
             fc_up_cutoff = 1, 
             fc_down_cutoff = 1,
             p_value_cutoff = 0.05,
             point_size_scale = "p_value_adjust_femalelcc_vs_femalecontrol") +
  scale_size_continuous(range = c(0.5, 5))
plot
ggsave(plot, filename = "statistical_analysis/female_LCC_vs_femal_control.pdf", width = 8, height = 7)
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
plot <- 
volcano_plot(object = object,
             fc_column_name = "fc_malelcc_vs_malecontrol",
             p_value_column_name = "p_value_adjust_malelcc_vs_malecontrol",
             add_text = TRUE,
             text_from = "Compound.name", 
             fc_up_cutoff = 1, 
             fc_down_cutoff = 1,
             p_value_cutoff = 0.05,
             point_size_scale = "p_value_adjust_malelcc_vs_malecontrol") +
  scale_size_continuous(range = c(0.5, 5))
plot
ggsave(plot, filename = "statistical_analysis/male_LCC_vs_male_control.pdf", width = 8, height = 7)
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
diff_metabolites_malelcc_vs_malecontrol <-
  object %>% 
  activate_mass_dataset(what = "variable_info") %>% 
  filter(p_value_adjust_malelcc_vs_malecontrol < 0.05) %>% 
  extract_variable_info()

readr::write_csv(diff_metabolites_malelcc_vs_malecontrol, 
                 file = "statistical_analysis/differential_metabolites.csv")
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
diff_metabolites_femalelcc_vs_femalecontrol <-
  object %>% 
  activate_mass_dataset(what = "variable_info") %>% 
  filter(p_value_adjust_femalelcc_vs_femalecontrol < 0.05) %>% 
  extract_variable_info()

readr::write_csv(diff_metabolites_femalelcc_vs_femalecontrol, 
                 file = "statistical_analysis/differential_metabolites.csv")
```

### Heatmap

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
temp_object <- 
object %>% 
  activate_mass_dataset(what = "sample_info") %>% 
  filter(sample_id %in% c(malecontrol, maleLCC)) %>% 
  activate_mass_dataset(what = "variable_info") %>% 
  filter(variable_id %in% diff_metabolites_malelcc_vs_malecontrol$variable_id) %>% 
  `+`(1) %>% 
  log(2) %>% 
  scale()
library(ComplexHeatmap)
h1 <- HeatmapAnnotation(class = extract_sample_info(temp_object)$group)
massstat::Heatmap(matrix = temp_object, 
                  name = "z-score", 
                  row_names_gp = gpar(cex = 0.4),
                  column_names_gp = gpar(cex = 0.2),
                  top_annotation = h1,
                  row_labels = extract_variable_info(temp_object)$Compound.name,
                  border = TRUE)
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
temp_object <- 
object %>% 
  activate_mass_dataset(what = "sample_info") %>% 
  filter(sample_id %in% c(femalecontrol, femaleLCC)) %>% 
  activate_mass_dataset(what = "variable_info") %>% 
  filter(variable_id %in% diff_metabolites_femalelcc_vs_femalecontrol$variable_id) %>% 
  `+`(1) %>% 
  log(2) %>% 
  scale()
library(ComplexHeatmap)
h1 <- HeatmapAnnotation(class = extract_sample_info(temp_object)$group)
massstat::Heatmap(matrix = temp_object, 
                  name = "z-score", 
                  row_names_gp = gpar(cex = 0.4),
                  column_names_gp = gpar(cex = 0.2),
                  top_annotation = h1, 
                  row_labels = extract_variable_info(temp_object)$Compound.name,
                  border = TRUE)
```


### Pathway enrichment analysis

Next, we want to know what pathways are enriched for male and female LCC tumor, respectively.

#### Load `KEGG` human pathway database

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
data("kegg_hsa_pathway", package = "metpath")
```

Remove the disease pathways

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
#get the class of pathways
pathway_class <- 
  metpath::pathway_class(kegg_hsa_pathway)
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
remain_idx <-
  pathway_class %>%
  unlist() %>%
  stringr::str_detect("Disease") %>%
  `!`() %>%
  which()
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
pathway_database <-
  kegg_hsa_pathway[remain_idx]
```

#### Pathway enrichment

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
kegg_id <-
  diff_metabolites_femalelcc_vs_femalecontrol$KEGG.ID 
kegg_id <-
  kegg_id[!is.na(kegg_id)]
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
result_female <- 
enrich_kegg(query_id = kegg_id, 
            query_type = "compound", 
            id_type = "KEGG",
            pathway_database = pathway_database, 
            p_cutoff = 0.05, 
            p_adjust_method = "BH", 
            threads = 3)
save(result_female, file = "pathway_enrichment/result_female")
```


```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
kegg_id <-
  diff_metabolites_malelcc_vs_malecontrol$KEGG.ID 
kegg_id <-
  kegg_id[!is.na(kegg_id)]
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
result_male <- 
enrich_kegg(query_id = kegg_id, 
            query_type = "compound", 
            id_type = "KEGG",
            pathway_database = pathway_database, 
            p_cutoff = 0.05, 
            p_adjust_method = "BH", 
            threads = 3)
save(result_male, file = "pathway_enrichment/result_male")
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
plot <- 
enrich_bar_plot(object = result_male,
                x_axis = "p_value",
                cutoff = 0.05)
plot
ggsave(plot,
       filename = "pathway_enrichment/male_pathway_bar.pdf",
       width = 7,
       height = 7)
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
plot <- 
enrich_scatter_plot(object = result_male,
                    y_axis = "p_value",
                    y_axis_cutoff = 0.05)
plot
ggsave(plot,
       filename = "pathway_enrichment/male_pathway_scatter.pdf",
       width = 7,
       height = 7)
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
plot <- 
enrich_bar_plot(object = result_female,
                x_axis = "p_value",
                cutoff = 0.05)
plot
ggsave(plot,
       filename = "pathway_enrichment/female_pathway_bar.pdf",
       width = 7,
       height = 7)
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
plot <- 
enrich_scatter_plot(object = result_female,
                    y_axis = "p_value",
                    y_axis_cutoff = 0.05)
plot
ggsave(plot,
       filename = "pathway_enrichment/female_pathway_scatter.pdf",
       width = 7,
       height = 7)
```

# Correlation network analysis

Next, we want to know if the correlation network are different between male and female tumor tissues.

## Female tumor

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
temp_object <- 
object %>% 
  activate_mass_dataset(what = "sample_info") %>% 
  filter(sample_id %in% c(femalecontrol)) %>% 
  activate_mass_dataset(what = "variable_info") %>% 
  filter(variable_id %in% diff_metabolites_femalelcc_vs_femalecontrol$variable_id) %>% 
  `+`(1) %>% 
  log(2) %>% 
  scale()

library(ggraph)
library(tidygraph)
graph_data_female <-
  convert_mass_dataset2graph(
    object = temp_object,
    margin = "variable",
    cor_method = "spearman",
    p_adjust_cutoff = 1,
    p_value_cutoff = 0.05,
    pos_cor_cutoff = 0,
    neg_cor_cutoff = 0
  ) %>% 
  mutate(Degree = centrality_degree(mode = 'all'))
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE, results = 'hide'}
library(extrafont)
loadfonts()
plot <- 
ggraph(graph = graph_data_female, layout = "kk") +
  geom_edge_fan(aes(color = correlation,
                    width = -log(p_value, 10)),
                show.legend = TRUE) +
  geom_node_point(aes(size = Degree)) +
  shadowtext::geom_shadowtext(aes(x = x, y = y,
                                  label = Compound.name),
                              bg.colour = "white",
                              colour = "black")+
  theme_graph() +
  scale_edge_color_gradient2(low = "darkblue", mid = "white", high = "red")
plot
ggsave(plot, filename = "pathway_enrichment/female_cor_network.pdf", width = 7, height = 7)
```

## Male tumor

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
temp_object <- 
object %>% 
  activate_mass_dataset(what = "sample_info") %>% 
  filter(sample_id %in% c(malecontrol)) %>% 
  activate_mass_dataset(what = "variable_info") %>% 
  filter(variable_id %in% diff_metabolites_malelcc_vs_malecontrol$variable_id) %>% 
  `+`(1) %>% 
  log(2) %>% 
  scale()

graph_data_male <-
  convert_mass_dataset2graph(
    object = temp_object,
    margin = "variable",
    cor_method = "spearman",
    p_adjust_cutoff = 1,
    p_value_cutoff = 0.05,
    pos_cor_cutoff = 0,
    neg_cor_cutoff = 0
  ) %>% 
  mutate(Degree = centrality_degree(mode = 'all'))
```

```{r,eval=TRUE,warning=FALSE, R.options="", message=TRUE, cache=TRUE}
plot <- 
ggraph(graph = graph_data_male, layout = "kk") +
  geom_edge_fan(aes(color = correlation,
                    width = -log(p_value, 10)),
                show.legend = TRUE) +
  geom_node_point(aes(size = Degree)) +
  shadowtext::geom_shadowtext(aes(x = x, y = y,
                                  label = Compound.name),
                              bg.colour = "white",
                              colour = "black")+
  theme_graph() +
  scale_edge_color_gradient2(low = "darkblue", mid = "white", high = "red")
plot

ggsave(plot, filename = "pathway_enrichment/male_cor_network.pdf", width = 7, height = 7)
```

# Session information

```{r,eval=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
sessionInfo()
```
