raw_data_processing.Rmd

---
title: "Raw MS data processing using massprocesser"
author:
- name: Xiaotao Shen (https://www.shenxt.info/)
date: "Created on 2020-04-01 and updated on `r Sys.Date()`"
output:
  html_document:
    df_print: paged
    toc: no
  pdf_document:
    toc: no
vignette: >
  %\VignetteIndexEntry{raw_data_processing}
  %\VignettePackage{massprocesser}
  % \VignetteEngine{knitr::rmarkdown}
  % \usepackage[utf8]{inputenc}
  %\VignetteEncoding{UTF-8}
---

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

***

---

# **Download demo data**

---

Raw data can be downloaded from the Google drive. [Link is here](https://drive.google.com/file/d/15OheYiPbqK7Bq8hdEvP9hdA0ER20P-xm/view?usp=sharing)

---

# **Data preparation**

---

*massprocesser* can used to processed the raw MS data for peak detection and alignment, and generate a peak table for next analysis. 

The MS raw data from Mass Spectrometry should be converted to `mzXML` or `mzML`format and then placed in different folders according to their class, for example `Subject`, `QC`, and `Blank` samples.

![](figures/Screen-Shot-2022-01-14-at-11.52.53-PM.png)

You can convert raw MS data to `mzXML` or `mzML` format using [`ProteoWizard` software](http://proteowizard.sourceforge.net/). And the parameter setting is shown the below figure:

![](figures/Screen-Shot-2022-01-14-at-11.54.24-PM.png)
---

# **Data organization**

---

All the `mzXML` format files should be placed in different folder according to sample type, such as `QC`, `Subject`, `Control`, `Case`, and `Blank`. 

![](figures/Screen-Shot-2022-01-14-at-11.52.53-PM.png)

---

# **Run `process_data()`** function

---

This function is used to peak picking and peak alignment.

```{r,eval=FALSE,warning=FALSE, R.options="", message=FALSE}
library(massprocesser)
library(tidyverse)
library(xcms)
library(MSnbase)
library(mzR)
```

## **Data processing**

Next, we use the `process_data()` function for peak detection and alignment.

We first do the positive mode.

```{r,eval=TRUE,warning=FALSE, message=FALSE,R.options="",cache=TRUE, message=FALSE}
massprocesser::process_data(
  path = "massprocesser_demo_data/POS/",
  polarity = "positive",
  ppm = 15,
  peakwidth = c(10, 60),
  snthresh = 5,
  noise = 500,
  threads = 6,
  output_tic = TRUE,
  output_bpc = TRUE,
  output_rt_correction_plot = TRUE,
  min_fraction = 0.5,
  fill_peaks = FALSE,
  group_for_figure = "QC"
)
```

Some important arguments:

* `ppm`: Peak detection ppm. See the `xcms` package.

* `peakwidth`: Peak width. It is dependent on your column and LC system. See the `xcms` package.

* `snthresh`: Signal to noise threshold. See the `xcms` package.

* `noise`: Noise cutoff. See the `xcms` package.

* `threads`: The core number for performing.

* `output_tic`, `output_bpc`, and `group_for_figure`: Do you want to output the TIC or BPC of samples? Some times there are a lot of samples, so you can set the `group_for_figure` as the group name, for example, you can set it as `QC`. 

* `min_fraction`: If one peak appears more than `min.fraction` sample in at least one group samples it will be kept.

* `fill_peaks`: Do you want to fill the missing values (NA)? 

Other parameters you can find here: `process_data()`.

---

# **Output results**

---

After all done, all the results are placed in a new folder named as `Result`.

![](figures/Screen-Shot-2022-01-15-at-8.42.04-AM.png)

## `object`

`object` is a `mass_dataset` class object which can be used for next analysis. See more about [`mass_dataset` class here](https://tidymass.github.io/massdataset/).

```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
load("massprocesser_demo_data/POS/Result/object")
library(massdataset)
object
```

## `Peak_table.csv`

The peak table for next analysis.

```{r,eval=TRUE, echo=FALSE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
data = readr::read_csv("massprocesser_demo_data/POS/Result/Peak_table.csv")
library(kableExtra)
library(magrittr)
kbl(data[1:500,]) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                fixed_thead = TRUE) %>% 
  scroll_box(width = "100%", height = "600px")
```

> Only top 500 features show.


## `Peak_table_for_cleaning.csv`

Some usefulness columns are deleted from the `Peak_table.csv`, this can be directory used for next data cleaning.

```{r,eval=TRUE, echo=FALSE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
data = readr::read_csv("massprocesser_demo_data/POS/Result/Peak_table_for_cleaning.csv")
library(kableExtra)
library(magrittr)
kbl(data[1:500,]) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                fixed_thead = TRUE) %>% 
  scroll_box(width = "100%", height = "600px")
```

> Only top 500 features show.


## `BPC.pdf` 

Base peak chromatogram.

![](figures/Screen-Shot-2022-01-15-at-11.36.11-AM.png)


Default it only shows the samples which are in the group you set in `group_for_figure`. If you want to show other samples, you can try to do use the below script.

```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
###load data
load("massprocesser_demo_data/POS/Result/intermediate_data/xdata2")
##set threads
if(tinytools::get_os() == "windows") {
  bpparam <-
    BiocParallel::SnowParam(workers = 4,
                            progressbar = TRUE)
} else{
  bpparam <- BiocParallel::MulticoreParam(workers = 4,
                                          progressbar = TRUE)
}
```

```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
###extract bpc plot object
bpc.plot <- xcms::chromatogram(object = xdata2,
                               aggregationFun = "max",
                               BPPARAM = bpparam)
```

```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
####plot_chromatogram to show the plot
##set the group_for_figure if you want to show specific groups. And set it as "all" if you want to show all samples.
plot <- plot_chromatogram(object = bpc.plot,
                          title = "BPC",
                          group_for_figure = c("Case", "QC"))
plot
```
```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
plot <- plot_chromatogram(object = bpc.plot,
                          title = "BPC",
                          group_for_figure = c("all"))
plot
```

```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
##if you want to show specific samples, set the `sample_for_figure` parameter
plot <- plot_chromatogram(object = bpc.plot,
                          title = "BPC",
                          sample_for_figure = c("sample_11", "sample_78"))
plot
```

```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
##set `interactive` as TRUE to get interactive plot
plot <- plot_chromatogram(object = bpc.plot,
                          title = "BPC",
                          sample_for_figure = c("sample_11", "sample_78"),
                          interactive = TRUE
                          )
plot
```

## `TIC.pdf` 

Total ion peak chromatogram

![](figures/Screen-Shot.png)


Default it only shows the samples which are in the group you set in `group_for_figure`. If you want to show other samples, you can try to do use the below script.

```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
###load data
load("massprocesser_demo_data/POS/Result/intermediate_data/xdata2")
##set threads
if(tinytools::get_os() == "windows") {
  bpparam <-
    BiocParallel::SnowParam(workers = 4,
                            progressbar = TRUE)
} else{
  bpparam <- BiocParallel::MulticoreParam(workers = 4,
                                          progressbar = TRUE)
}
```

```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
###extract tic plot object
tic.plot <- xcms::chromatogram(object = xdata2,
                               aggregationFun = "sum",
                               BPPARAM = bpparam)
```

```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
####plot_chromatogram to show the plot
##set the group_for_figure if you want to show specific groups. And set it as "all" if you want to show all samples.
plot <- plot_chromatogram(object = tic.plot,
                          title = "TIC",
                          group_for_figure = c("Case", "QC"))
plot
```
```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
plot <- plot_chromatogram(object = tic.plot,
                          title = "TIC",
                          group_for_figure = c("all"))
plot
```

```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
##if you want to show specific samples, set the `sample_for_figure` parameter
plot <- plot_chromatogram(object = tic.plot,
                          title = "TIC",
                          sample_for_figure = c("sample_11", "sample_78"))
plot
```

```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
##set `interactive` as TRUE to get interactive plot
plot <- plot_chromatogram(object = tic.plot,
                          title = "TIC",
                          sample_for_figure = c("sample_11", "sample_78"),
                          interactive = TRUE
                          )
plot
```

## `RT correction plot.pdf`

RT correction.

![](figures/Screen-Shot-1.png)

Default it only shows the samples which are in the group you set in `group_for_figure`. If you want to show other samples, you can try to do use the below script.

```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
###load data
load("massprocesser_demo_data/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.
plot = 
plot_adjusted_rt(object = xdata2, 
                 group_for_figure = c("Case", "QC"))
plot
```
```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
plot = 
plot_adjusted_rt(object = xdata2, 
                 group_for_figure = "all")
plot
```

```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
plot = 
plot_adjusted_rt(object = xdata2, 
                sample_for_figure = c("sample_11", "sample_78"))
plot
```

```{r,eval=TRUE, echo=TRUE,warning=FALSE, R.options="", message=FALSE, cache=TRUE}
plot = 
plot_adjusted_rt(object = xdata2, 
                sample_for_figure = "all", 
                interactive = TRUE)
plot
```

## `intermediate_data` 

This is a folder which contains some intermediate data. If you want to re-run the data processing, please delete this folder first.

![](figures/Screen-Shot-2022-01-15-at-9.03.43-AM.png)

---

# **Session information**

---

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