```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
```



```{r}
library(tidyverse)
library(readxl)
```


## Dataset and description

[**KOMP Plasma Metabolomics Dataset**](https://www.kaggle.com/desertman/komp-plasma-metabolomics-dataset)

Mouse knockouts facilitate the study of gene functions. Often, multiple abnormal phenotypes are induced when a gene is inactivated. The International Mouse Phenotyping Consortium (IMPC) has generated thousands of mouse knockouts and cataloged their phenotype data. We have acquired metabolomics data from 220 plasma samples from 30 unique mouse gene knockouts and corresponding wild-type mice from the IMPC. To acquire comprehensive metabolomics data, we have used liquid chromatography (LC) combined with mass spectrometry (MS) for detecting polar and lipophilic compounds in an untargeted approach. We have also used targeted methods to measure bile acids, steroids and oxylipins. In addition, we have used gas chromatography GC-TOFMS for measuring primary metabolites. The metabolomics dataset reports 832 unique structurally identified metabolites from 124 chemical classes as determined by ChemRICH software

In this demo we will use only the metabolites quantified in the targeted assays. The belong to the chemical classes of _bile acids, steroids and oxylipins_


Here we get the data
```{r}
load("../data/KOMP_data_targeted.RData") ## the initial data were saved in RData format
```


We have now three files

* the description of the metabolites
* the description of the samples
* the actual data matrix

## Metabolite Description

The amount of info in the table is uncommon

```{r}
colnames(metabolite_meta)
```

What is interesting for us is to see how much effort has been made to precisely identify the compounds. The problem of "naming" is typical of analytical chemistry and metabolomics.

- common chemical names are ambiguous
- unique identifiers are impractical (look to the InChiKey for example)

```{r}
metabolite_meta$InChiKey[1]
```
Some of the identifier refer to publicly available databases. This is the general trend of this scientific discipline.

**KEGGID** refers to _kegg_ so it can be used to link the metabolite to its position within metabolic networks.


## Sample Description

The second tibble contains the description of the sample metadata 

```{r}
head(sample_meta_data)
```

The field are auto explanatory. The MouseID is the code we need to associate the sample to the class

First of all we give a look to the experimental design, just to understand the numbers we are dealing with

```{r}
table(sample_meta_data$Gender,sample_meta_data$Genotype)
```

* Is true that we have many animals ... but the number of samples per class is really limited
* Luckily enough the design is "orthogonal" (in the sense that for each genotype we have both genders) and balanced 
* The small number of samples per class is in a sense ruling out the possibility of using with a good efficacy non parametric tests ...at least to compare the genotypes.

## The Data Matrix

```{r}
head(DM[,1:20])
```


* it is switched compared to the "standard" shape. This is often done because the number of variables is larger then the number of samples. Since R wants the reverse we transpose ...

Since our mouse ID is a number without the leading "KOMP_" we fix this on the fly


```{r}
DM <- DM %>% 
  column_to_rownames("CompoundID") %>% 
  t(.) %>% 
  as_tibble(rownames = "Sampname") %>% 
  separate(Sampname, c("One","MouseID"),"_") %>% ## split the column
  dplyr::select(-One) 
  
```


Now the matrix looks like this 

```{r}
head(DM[,1:20])
```

Now I fix the sample name in the `sample_meta`

```{r}
sample_meta_data <- sample_meta_data %>% 
  mutate(MouseID = as.character(MouseID))
```


* As expected we have missing data, which most likely show up for metabolites with very low concentrations. 
* This statement is not always true because the detection limit con vary between different compounds.


To check that we organize the matrix in a different way also joining the concentration data with the sample metadata

```{r}
DM_nest <- DM %>% 
  left_join(sample_meta_data) %>% 
  pivot_longer(starts_with("CPD"), names_to = "CompoundID", values_to = "I") %>% 
  nest(data = -CompoundID)
```

This data structure highlights the role of each metabolite 

```{r}
head(DM_nest)
```

the intensities of each metabolite are stored in the "data" column

```{r}
head(DM_nest$data[[1]])
```

The new tibble can be efficiently used to inspect the data in a univariate perspective ...


## Missing Values

For example, the number of NAs in relation of the median intensity can be visualized as

```{r}
DM_nest %>% 
  mutate(nnas = map_dbl(data,function(p) sum(is.na(p$I))),                  ## calculate the univariate statistics "on the fly"
         medI = map_dbl(data,~median(.x$I, na.rm = TRUE))) %>% 
  filter(nnas > 0) %>% ## we keep only features with at least one NA!
  ggplot() +   ## plot them
  geom_point(aes(x = medI, y = nnas)) + 
  scale_x_log10() + 
  theme_light()
  
```

Which actually shows what we expected. High number of missing values are visible only with low median intensities ... note the log scale in x. 

This is telling us that there is a distinction of metabolites in two groups showing a marked difference in concentration

Is this effect dependent on the class of metabolite? 


```{r}
library(plotly)


ggplotly(DM_nest %>% 
  left_join(metabolite_meta) %>% 
  mutate(nnas = map_dbl(data,~sum(is.na(.x$I))),                  ## calculate the univariate statistics "on the fly"
         medI = map_dbl(data,~median(.x$I, na.rm = TRUE))) %>% 
  filter(nnas > 0) %>% ## we keep only features with at leas on NA!
  ggplot(aes(x = medI, y = nnas, col = Assay, text = CompoundName)) +   ## plot them
  geom_point() + 
  scale_x_log10() + 
  scale_color_brewer(palette = "Set2") + 
  theme_light(), tooltip="text")
```

That's strange ... the largest number of NAs is corresponding to lower intensities for the Bile Acids_Steroids and not for Oxylipins. This is curious 

It is worth highlighting the name of the metabolites showing a number of NAs larger than 50 ...

```{r}
DM_nest %>% 
  left_join(metabolite_meta) %>% 
  mutate(nnas = map_dbl(data,~sum(is.na(.x$I))),                  ## calculate the univariate statistics "on the fly"
         medI = map_dbl(data,~median(.x$I, na.rm = TRUE))) %>% 
  filter(nnas > 0) %>% ## we keep only features with at leas on NA!
  ggplot(aes(x = medI, y = nnas, col = Assay)) +   ## plot them
  geom_point(size = 2) + 
  geom_text(aes(label=ifelse(nnas>50,CompoundName,'')), size = 4, nudge_y = -3) + 
  scale_x_log10() + 
  scale_color_brewer(palette = "Set2") + 
  theme_light()
```


As we discussed in the lecture, it is difficult to decide which metabolites have too many missing values ...

When we have a "limited" number of variables the best approach is to plot the data. With the `DM_nest` we can do that quite efficiently

To do that I create a custom plotting function which takes as input the tibble of the data and a potential title

```{r fig.height=3, fig.width=8}
myplot <- function(t,l = ""){
  p <- t %>% 
  ggplot() + 
  geom_jitter(aes(x = Genotype, y = I, color = Gender, shape = Zygosity), width = 0.1) + 
  scale_shape_manual(values = c(19,3,1)) + 
  scale_y_log10() + 
  scale_color_brewer(palette = "Set1") + 
  theme_light() + 
  ggtitle(l) + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
  
  p
}


myplot(DM_nest$data[[1]], "test")

```

And I recursively apply it to the nested data tibbles ...


```{r}
DM_nest <- DM_nest %>%
  left_join(metabolite_meta %>% dplyr::select(CompoundID,CompoundName)) %>% 
  mutate(rawplots = map2(data,CompoundName, ~myplot(.x,.y)))

## adding compound names as names oif the list of plots

names(DM_nest$rawplots) <- DM_nest$CompoundName

```


Now I can visualize (or save) the plot I want ...


This one, in particular show the trend of one of the features with a large fraction of NAs
```{r fig.height=3, fig.width=8}
DM_nest$rawplots[[5]]
```

As it can be expected, progesterone is lower in males, and, most likely close to the detection limit. It is worth remarking that progesterone was the steroid showing the lower intensity. Its behavior is then not odd.

Let's now look to PGD2

```{r fig.height=3, fig.width=8}
DM_nest$rawplots$PGD2
```
The trend of this metabolite is really wired. Its large fraction of missing values do not seems to be associated to the intensity of the signal


As far as I can see (but you could check better ... ;-) ), the distribution of the missing data does not seems to be highly associated with the study factors. 

Considering that we have 220 samples the level of missingness is never extremely severe. 


**Something For you**

1. play around trying to visualize the other metabolites showing high levels of missingness


## Variable Distribution

The need of displaying the intensity in log scale is already suggesting that the distribution of the intensities is far from being normal. This is not unexpected in metabolomics. In real life what I do is to work almost invariably with log-transformed datasets.

For sure

1. The low number of samples for each genotypes makes a class based assessment completely hopeless
2. What we could do is to look to the different variables highlighting the gender and the zygosity.


We can follow the spirit of the previous plot by defining a function and using it to produce a full list of plots ...

An alternative could be to `unnest` the data and rely on faceting ...


```{r fig.height=14, fig.width=8}
DM_nest %>% 
  dplyr::select(CompoundName,data) %>% 
  unnest(data) %>% 
  ggplot(aes(sample = I)) +
  stat_qq(aes(col = Gender, shape = Zygosity)) + 
  scale_shape_manual(values = c(19,3,1)) + 
  scale_color_brewer(palette = "Set1") + 
  stat_qq_line() + 
  facet_wrap(~CompoundName, scales = "free_y", ncol = 5) +
  theme_light() + 
  theme(aspect.ratio = 1, 
        legend.position = "bottom")
```


Ok, the plot is "ugly" (maybe creating a set of plots as we did before could have been more pleasant), but is really informative.

* No sign of normal distribution
* For some metabolites M and F are really different ... and for 12(13)Ep-9-KODE the trend is really interesting.

Let's focus on it
```{r}
DM_nest %>% 
  dplyr::select(CompoundName,data) %>% 
  unnest(data) %>% 
  filter(CompoundName == "12(13)Ep-9-KODE") %>% 
  ggplot(aes(sample = I)) +
  stat_qq(aes(col = Gender, shape = Zygosity)) + 
  scale_shape_manual(values = c(19,3,1)) + 
  scale_color_brewer(palette = "Set1") + 
  stat_qq_line() + 
  facet_wrap(~CompoundName, scales = "free_y", ncol = 5) +
  theme_light() + 
  theme(aspect.ratio = 1, 
        legend.position = "bottom")
```


```{r fig.height=3, fig.width=8}
DM_nest$rawplots$`12(13)Ep-9-KODE`
```

The interpretation for that would, most likely require a chat with a good biochemist...

Let's look to the global effect of a log transformation on the q-q plots

```{r fig.height=14, fig.width=8}
DM_nest %>% 
  ungroup() %>% 
  dplyr::select(CompoundName,data) %>% 
  unnest(data) %>% 
  ggplot(aes(sample = log10(I))) +
  stat_qq(aes(col = Gender, shape = Zygosity)) + 
  scale_color_brewer(palette = "Set1") + 
  stat_qq_line() + 
  facet_wrap(~CompoundName, scales = "free_y", ncol = 5) +
  theme_light() + 
  theme(aspect.ratio = 1, 
        legend.position = "bottom")
```

The situation largely improves, even if some compounds are still showing strange behaviors ...

**Something For you**

1. Could you look to the metabolites showing outlying behaviors? Is this a matter of genotype? 


```{r}
DM_nest$rawplots$`Taurodeoxycholic acid`
```

Get out the samples wiuth intensity higher than 2000


```{r}
names(DM_nest$data) <- DM_nest$CompoundName
```

```{r}
critical_1 <- DM_nest$data$`Taurodeoxycholic acid` %>% 
  filter(I > 1500) %>% 
  pull(MouseID)

critical_1
```



## Multivariate Visualization - PCA

The previous univariate journey should have been already pointing out some interesting aspects of the dataset.

The subsequent step is to look to visualize the multivariate structure of the data to spot larger scale patterns.

In order to do that it is necessary to impute the data matrix, substituting NAs with meaningful numbers.

1. The number should have a variability
2. The number should not know anything about the experimental design
3. The number should have a reasonable analytical meaning

_My choice here is to work variable wise and replace the NAs with a random number drawn from an uniform distribution spanning from 0 to the  minimum value measured for that variable_

The rationale behind this choice is that 

* a concentration is a positive number
* everything I cannot measure is equally likely


In my practical example imputation could be performed both in the `DM` and in its nested version `DM_nest`. Since we are going to do PCA, I'll work on the unnested data matrix.

First of all I need a function to perform the imputation of a vector

```{r}
myimputer <- function(v){
  if (sum(is.na(v)) == 0) {  ## If I do not have NAs please leave the vector unchanged
    return(v)
  } else {
    napos <- which(is.na(v))  ## position of the NAs in the vector
    newval <- runif(length(napos),0,min(v, na.rm = TRUE))  ## calculate the random numbers I should put in place of the NAs
    out <- v
    out[napos] <- newval
    return(out)
  }
  
}
```

Now we apply it to the full set of columns

```{r}
DM_i <- DM %>% 
  mutate(across(starts_with("CPD"), ~ myimputer(.x)))
```


As we have seen, log transformation is helpful ...


```{r}
DM_i <- DM_i %>% 
  mutate(across(starts_with("CPD"), ~log10(.x)))
```

And now PCA!

```{r}
library(FactoMineR)
library(factoextra)
```

```{r}
myPCA <- PCA(DM_i %>% 
               column_to_rownames("MouseID") %>% 
               dplyr::select(starts_with("CPD")), 
             graph = FALSE, scale.unit = FALSE   ## to avoid factmineR producing a plot by default
             )
```


```{r}

library(plotly)

fviz_pca_ind(myPCA, 
             habillage = factor(sample_meta_data$Gender), ## factor which will be used to color the dots
             # geom = "point", 
             pointsize = 2,
             axes = c(1,2)) + 
  scale_color_brewer(palette = "Set1")
```
Even if the amount of variability captured in this representation is only the 30%, the previous plot shows the presence of a separation between males and females.

It would be nice to have a similar plot highlighting the role of the different genetic backgrounds, but we do not have enough colors to put everything in the same plot. 



```{r}
## we get out the scores
PC_12_scores <- facto_summarize(myPCA, "ind", 
                                axes = c(1,2),
                                result = "coord"
                                ) %>%    ## get the data from the PCA object
  mutate(name = as.character(name)) %>%  
  as_tibble() %>% 
  ## transform the sample name to a character
  dplyr::rename("MouseID" = "name") %>%     ## align the name with the correct name
  left_join(sample_meta_data)        ## join with the sample data

PC_12_scores
```

```{r fig.height=6, fig.width=8}
PC_12_scores %>% 
  ggplot() + 
  geom_point(data = PC_12_scores %>% select(Dim.1,Dim.2), 
             mapping = aes(x= Dim.1, y = Dim.2), col = "gray80", size = 0.2) + 
  geom_hline(yintercept = 0, lty = 2)+ 
  geom_vline(xintercept = 0, lty = 2)+ 
  geom_point(aes(x = Dim.1, y = Dim.2, col = Gender)) + 
  facet_wrap(~Genotype, ncol = 6) + 
  theme_bw() + 
  theme(aspect.ratio = 1)
```

## Something for You

Suppose that I'm reall yinterested in the comparison between Pipox and Null (wt) ...
Could you make a PCA plot only of these classes of samples?






So, at least in the PCA, no clear separation between the genotypes

This is not unexpected considering the results of our univariate investigation


What about the variables?

```{r}
fviz_pca_biplot(myPCA, 
             habillage = factor(sample_meta_data$Gender), 
             geom = "point", 
             pointsize = 2,
             axes = c(1,2)) + 
  scale_color_brewer(palette = "Set1")
```

The biplot here shows that females show an higher level of the large majority of the metabolites. 

Since the partial separation is visible along PC1 we could visualize the variables which show the higher contribution to that direction


```{r}
fviz_contrib(myPCA, "var", axes = 1)
```

These plots are really informative, let's give a look inside the PCA object to see how these number could be get out. 


This list contains the information about the samples (the scores) so it can be used to reconstruct the score plot. What one needs are the coordinates 

```{r}
str(myPCA$ind)
```

The `dist` and `cos2` can be used to assess how far a sample is from the PCA plane (or projection). This is important to identify potential outliers.


The list

```{r}
str(myPCA$var)
```

instead contains the info about the variables, their loadings and the correlation with the principal components.

To extract the previous information in a more programmatic way, `FactoMineR` has a specific function

```{r}
head(facto_summarize(myPCA, "var", axes = 1))
```
which returns as a table the content of the PCA object.

By using the content of the myPCA object it is possible to "reconstruct" the PCA plots with `ggplot` and this can be handy if specific representation of the data are needed.

Suppose, for example, that you would like to construct a more informative variable importance plot which shows the correct metabolite name and also color the bars according to the chemical/metabolic class of the compound.


```{r fig.height=8, fig.width=6}
facto_summarize(myPCA, "var", axes = 1) %>% 
  as_tibble() %>% 
  dplyr::select(name, `Dim.1`) %>% 
  left_join(metabolite_meta, by = c("name" = "CompoundID")) %>% 
  arrange(`Dim.1`) %>% 
  mutate(CompoundName = factor(CompoundName, levels = CompoundName)) %>% 
  ggplot() + 
  geom_point(aes(x = CompoundName,  y= `Dim.1`, col = Assay), size = 2) + 
  geom_segment(aes(x = CompoundName,  yend = `Dim.1`, col = Assay, xend = CompoundName, y = 0)) + 
  scale_color_brewer(palette = "Set2") + 
  coord_flip() + 
  theme_light()
```

Which can be of help for the interpretation of the results

**Something For you**

1. Could you try a different imputation strategy? Is this affecting the results of your exploratory data analysis?
2. Just play and ask!


## Biomarker Discovery

What we did so far is called **exploratory data analysis**. To write the paper we should also find the variables which are showing a significant effect of the design factors. 

**Notes**

1. One could look for biomarkers with univariate and multivariate methods. 
2. In a univariate perspective, to avoid the risks of non matching the assumption of many statistical tests it would be safer to rely on non-parametric approaches
3. As far as the difference in genotype is concerned, the number of samples is limited. Most likely, there we should use parametric tests
4. Since we are performing multiple tests, we should remember that false positives will be present!

For illustrative purposes let's focus on the Female-Male Comparison.

```{r}
## here i nest the imputed data matrix, most likely I could do the testing also in the "non" imputed one
## we also add the gender of the animal
DM_i_nest <- DM_i %>% 
  pivot_longer(starts_with("CPD"), names_to = "CompoundID", values_to = "I") %>% 
  left_join(sample_meta_data) %>% 
  nest(data = -c(CompoundID))
```

Now we run the tests on the metabolites

```{r}
DM_i_nest <- DM_i_nest %>% 
  mutate(wilcoxon = map(data, ~ wilcox.test(I ~ Gender, data = .x, exact = FALSE)))
```


If we want to see one of the tests ... 

```{r}
DM_i_nest$wilcoxon[[10]]
```

Now we extract the p-values

```{r}
DM_i_nest <- DM_i_nest %>% 
  mutate(ps = map_dbl(wilcoxon, ~.x$p.value))
```


since we have potentially false positive it is always a good idea to look to the distribution of the p-values

```{r}
hist(DM_i_nest$ps, breaks = 20)
```

We see a clear enrichment of the "low" p values. One could ask what would be the distribution of the p-values if there would be no difference between males and females.

Let's try so simulate this case. To do that I scramble the intensities for each metabolite, so I'll destroy the assiciation between the intensity and the gender

```{r}
DM_i_nest <- DM_i_nest %>% 
  mutate(permuted_data = map(data, ~ .x %>% mutate(I = sample(I))))
  
```

Let's now run a wikoxon test on the permuted data


```{r}
DM_i_nest <- DM_i_nest %>% 
  mutate(wilcoxon_p = map(permuted_data, ~ wilcox.test(I ~ Gender, data = .x, exact = FALSE))) %>% 
  mutate(ps_p = map_dbl(wilcoxon_p, ~.x$p.value))
```

And now we plot the histogram

```{r}
hist(DM_i_nest$ps_p, breaks = 20)
```

This is clearly demonstrating that there is a difference between males and females and that this difference is "statistically significant"


Now we would like to rank the variables on the base of their contrast on the two classes ... but as we know it is not correct to use the _p-value_ for that. 

We need to calculate the **effect size**! Now I'll use the Cohen's d. A word of caution here, this type of effect size is parametric ... but I've been stressing that non parametric analysis here should be preferred. 
Right. A non parametric measure of the effect size should be better. A possibility could be to calculate the ratio between the difference between the two medians and divide it by the larger interquartile range. You could be creative there ...


Let's stick to Cohen's d. In the same time I'll also calculate the p-values corrected for multiplicity

```{r}
library(effsize)

DM_i_nest <- DM_i_nest %>%  
  mutate(pcorr = p.adjust(ps, method = "BH")) %>%    ## here I'm calculating the corrected p-values
  mutate(cohend = map_dbl(data, ~cohen.d(I ~ Gender, data = .x)$estimate))  ## here I'm calculating the effect size
```

Now we use the effect size and the p value to create a volcano plot


```{r}
library(plotly)
ggplotly(DM_i_nest %>% 
           left_join(metabolite_meta %>% select(CompoundID,CompoundName, Assay)) %>% 
           ggplot(aes(text=CompoundName)) + 
           geom_hline(yintercept = 1.3, lty = 2, col = "red") + 
           geom_point(aes(x = cohend, y = -log10(ps), fill = Assay), size = 3, pch = 21) + 
           theme_bw(), tooltip="text")
```


Now we focus on the compound which are significant at the 0.01 level and we sort the results by effect size ...


```{r}
markers <- DM_i_nest %>% 
  filter(pcorr < 0.01) %>% 
  arrange(desc(abs(cohend))) %>% 
  left_join(metabolite_meta)

head(markers)
```


This is the list of "biomarkers", let's plot one of them ...

```{r fig.height=3, fig.width=8}
DM_nest$rawplots[[54]]
```
Which exactly shows what we expect.

**Notes**

1. This is the metabolite showing the larger "contrast" so our interpretation could start from there. 
2. This is also the metabolite showing up in the PCA ;-)


## Variable Correlation

The association between the intensity of two metabolites in the samples is often considered an interesting evidence in the perspective of the interpretation. Metabolites are indeed linked in metabolic pathways and their concentrations are not independent

As we have seen yesterday, however, correlation can be there only by chance so care have to be taken to jump from the observation of an association to the conclusion of a causal relation.

These caveats are particularly true in the case of untargeted metabolomics assays and can be somehow taken more easily as far as targeted metabolomics is concerned.

The best tool to inspect variable correlation in R is the `pairs` plot. With `pairs` one really sees the experimental data and it is easy to spot subgroups and outlying samples.

Unfortunately, the pairs visualization is manageable only with a small number of variables (something like 15). For this reason it is a resource which can be exploited after the identification of the most relevant _biomarkers_

```{r}
## here I get the names of the first 5
focuscp <- markers %>% 
  slice(1:5) %>% 
  pull(CompoundID)
```


```{r fig.height=6, fig.width=6}
## This is a littel of magic ...
mypalette <- c("#F24D2970","#1C366B70")
names(mypalette) <- unique(sample_meta_data$Gender)

par(pty="s")

DM_i %>% 
  dplyr::select(all_of(focuscp)) %>% 
  as.matrix(.) %>% 
  pairs(., col = mypalette[sample_meta_data$Gender], pch = 19)


```

**Notes**

1. We see correlations between the compounds. The biochemist will be happy. In particular we see an indication of a change in correlation between different oxylipins. This should be checked against metabolic pathways. 
2. It would be interesting to use the previous approach to see if some of the outlying samples are associated to a specific genotype


**Something For you**

1. Could you try to use ANOVA to look to metabolites potentially different in the genotypes?
2. Do you think that the ANOVA assumption here are fulfilled?
3. Run a similar analysis on the rubus targeted data

