---
title: "Mass calculation"
output: 
  bookdown::html_document2:
    base_format: rmarkdown::html_vignette
    toc: true
    number_sections: true
    always_allow_html: yes
pkgdown: 
    as_is: true
vignette: >
  %\VignetteIndexEntry{Mass calculation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# Comparison of results with weighted mean and without

As it is described in `vignette("datafiles")` article, the data within each replicate of the experiment is aggregated using weighted mean. Below we present the effect it has on the results.

The whole workflow is described in the mentioned article. Here the focus is on the first aspect of the processing.

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(dev = "png", dev.args = list(type = "cairo-png"))


library(tidyr)
library(dplyr)
library(HaDeX2)
library(gridExtra)
library(ggplot2)

```

```{r echo = FALSE}

f_calculate_mass <- function(dat){

  proton_mass <- 1.00727647

  dat %>%
    mutate(exp_mass = Center*z - z*proton_mass) %>%
    select(-Center, -z, -Modification) %>%
    group_by(Sequence, Start, End, MHP, MaxUptake, State, Exposure, Protein, File) %>%
    summarize(avg_exp_mass = weighted.mean(exp_mass, Inten, na.rm = TRUE)) %>%
    ungroup(.) %>%
    group_by(Sequence, Start, End, MHP, MaxUptake, State, Exposure, Protein) %>%
    summarize(mass = mean(avg_exp_mass, na.rm = TRUE),
              err_mass = coalesce(sd(avg_exp_mass, na.rm = TRUE)/sqrt(sum(!is.na(avg_exp_mass))), 0),
              num = (sum(!is.na(avg_exp_mass)))) %>%
    ungroup(.) %>%
    arrange(Start, End, Start, Exposure) %>%
    as.data.frame()

}

f_calculate_mass_no_inten <- function(dat){

  proton_mass <- 1.00727647

  dat %>%
    mutate(exp_mass = Center*z - z*proton_mass) %>%
    select(-Center, -z, -Modification) %>%
    group_by(Sequence, Start, End, MHP, MaxUptake, State, Exposure, Protein, File) %>%
    summarize(avg_exp_mass = mean(exp_mass, na.rm = TRUE)) %>%
    ungroup(.) %>%
    group_by(Sequence, Start, End, MHP, MaxUptake, State, Exposure, Protein) %>%
    summarize(mass = mean(avg_exp_mass, na.rm = TRUE),
              err_mass = coalesce(sd(avg_exp_mass, na.rm = TRUE)/sqrt(sum(!is.na(avg_exp_mass))), 0)) %>%
    ungroup(.) %>%
    arrange(Start, End, Start, Exposure) %>%
    as.data.frame()

}

```

For the analysis, we use the example file from HaDeX, measurements of alpha subunit from elongation factor eEF1.

As described, the first step is to transform the $Center$ value (geometric centroid of the isotopic envelope for given peptide in a given state in given time point) and then to aggregate the values measured for different charge values. This is all done within each replicate.

The aggregation in our workflow is a weighted mean, with $Inten$ (intensity) values as weights, as shown below.

```{r eval=FALSE}
avg_exp_mass = weighted.mean(exp_mass, Inten, na.rm = TRUE)
```

Later, the results from replicates are aggregated using the mean, and their uncertainty is calculated as the standard deviation of the mean. 

How the weighted mean change the result in comparison with the simple mean? Let's see.

```{r include=FALSE}
dat_no_weight <- f_calculate_mass_no_inten(alpha_dat) %>%
  arrange(Start, End, State, Exposure) %>%
  mutate(source = "no_weight") %>%
  select(Sequence, Start, End, State, Exposure, mass, err_mass, source)

dat_weight <- f_calculate_mass(alpha_dat) %>%
  arrange(Start, End, State, Exposure) %>%
  mutate(source = "weight") %>%
  select(Sequence, Start, End, State, Exposure, mass, err_mass, source)

tmp <- bind_rows(dat_no_weight, dat_weight) %>%
  gather(type, value, -Sequence, -Start, -End, -source, -State, -Exposure) %>%
  spread(source, value) %>%
  mutate(diff = (weight-no_weight))
```

Below are calculated values of mass in two approaches and the difference between them for the example peptide. 

```{r echo=FALSE}
tmp %>%
  filter(type == "mass") %>%
  select(-type) %>%
  filter(Sequence == "GFGDLKSPAGL")
```

Below are calculated values of the uncertainty of mass in two approaches and the difference between them for the example peptide. 

```{r echo=FALSE}
tmp %>%
  filter(type == "err_mass") %>%
  select(-type) %>%
  filter(Sequence == "GFGDLKSPAGL")
```

And for the whole set of peptides, the average mass difference is:

```{r echo=FALSE}
mean(tmp[ tmp[["type"]] == "mass" , "diff"])
```
and uncertainty difference is:

```{r echo=FALSE}
mean(tmp[ tmp[["type"]] == "err_mass" , "diff"])
```
Below are the histograms of these differences.

First, differences between mass:

```{r message=FALSE, warning=FALSE, echo=FALSE}
ggplot(filter(tmp, type == "mass"), aes(diff)) +
  geom_histogram() +
  labs(title = "Differences between mass",
       x = "Difference")
```

Secondly, the differences between uncertainties:


```{r message=FALSE, warning=FALSE, echo=FALSE}
ggplot(filter(tmp, type == "err_mass"), aes(diff)) +
  geom_histogram() +
  labs(title = "Differences between uncertainties of mass",
       x = "Difference")
```