MTXQCvX_part2_template_params.Rmd

---
#title: "MTXQCvX2 Part2: Post-Processing -
#thanks: "Template MTXQCvX2 part 2 provided by Christin Zasada, Kempa Lab"
author:
- affiliation: LAB, LAB LOCATION 
  name: NAME ONE 
- affiliation: LAB, LAB LOCATION
  name: NAME TWO
biblio-style: apsr
output:
  pdf_document:
    citation_package: natbib
    fig_caption: yes
    keep_tex: yes
    toc: yes
    latex_engine: pdflatex
    template: config_files/textemplate.tex
  html_document: default
fontsize: 9pt
geometry: margin=1in
keywords: MTXQCvX, post-processing, data visualisation
date: "`r format(Sys.time(), '%V %d, %Y')`"
fontfamily: mathpazo
params:
  subs:
    input: text
    label: "Run on MTXQC project folder:"
    value: ""
  analysis:
    input: select
    choices:
      - stringent
      - less_stringent
    value: stringent
  ca_fac:
    label: "Internal Standard"
    input: select
    choices: 
      - QC
      - all
    value: QC
  quant:
    input: select
    label: 'Perform post-processing based on:'
    choices:
    - PeakArea
    - absconc
    - corr_absconc
    - Conc_pmio
    - Conc_microM
    - sumA_Conc
    - sumA_Conc_pmio
    - sumA_Conc_microM
    - IntStd_Conc
    - IntStd_Conc_pmio
    - IntStd_Conc_microM
    - IntStd_sumA_Conc
    - IntStd_sumA_Conc_pmio
    value: IntStd_Conc_pmio
  par1:
    input: text
    label: "Statistics: Define parameter 1 (plot: x-axis)."
    value: ""
  par2:
    input: text
    label: "Statistics: Define parameter 2 (plot: subgrouping, fill-value)."
    value: ""
  par3:
    input: text
    label: "Statistics: Define parameter 3 (plot: sub-figure panel)."
    value: ""
  par4:
    input: text
    label: "Statistics: Define parameter 4."
    value: ""
  updated:
    input: select
    choices: [none, PeakArea, Incorporation, both]
    label: "Incorporation manually validated data:"
    value: none
  mypal:
    input: select
    choices:
      - main
      - cheerful
      - hot
      - greens
      - sleek
      - clean
      - cool
      - art
    label: "Select your color scheme."
    value: greens
abstract: This document provides the post-processing of the output files derived from MTXQC part1; including the determination of statistics and visualisation.
  Here you have to work the most to adapt this file to your own needs and the
  experimental setup. Take care of decent plots. Most probably you're going to
  reuse this script for other projects. Feel free to extend this script as much
  as you want and to modify or define colors. Check carefully if you defined the 
  grouping for the statistics parts correctly!
---

---
title: `r paste("MTXQCvX2 part2 - PostProcessing ", params$analysis)`
---
# Summary

*Summarise the most important findings or questions.*


#Project initiation

```{r init_doc, echo=FALSE, warning=FALSE}
#mypal-def for color
mypal = params$mypal

#Source R-Files
source('R/MTXQC_fcn.R')
source('R/MTXQC_fcn_metabolicprofile.R')
source('R/MTXQC_config.R')
source('R/MTXQC_colors.R')
source("R/MTXQC_pck.R")

set_input = "input/"
set_output = "output/"

## subfolder for postprocessing
#Defining pp-folder
set_pp = paste0(Sys.Date(), "_PostProc_", params$analysis, "_",params$ca_fac, "_" ,params$quant, "/")

#directory definition and figure_name definition
if (params$subs == "") {
  path_setup = ""
  set_fig = paste0(path_setup, set_output, set_pp,'figure/MTXQCp2_', params$quant, '-')
} else {
  path_setup = paste0(params$subs, "/")
  set_fig = paste0(path_setup, set_output, set_pp,'figure/MTXQCp2_', params$tissue, '_', params$quant, '-')
}

#create folder: PostProcessing if not present
ifelse(!dir.exists(file.path(paste0(path_setup, set_output, set_pp))), 
			 dir.create(paste0(path_setup, set_output, set_pp)), 
			 sprintf('The folder exists already! Be careful overwriting files and figures!'))

dir.create(paste0(path_setup, set_output, set_pp, "figure/"))

#set path for figure export and size
knitr::opts_chunk$set(fig.width = 7, fig.align = 'center', fig.height = 8, 
                      fig.path = set_fig, 
                      echo = FALSE,  #show R code
                      warning = FALSE, #show warnings
                      message = TRUE) #show messages
                      
```                    

## Report parameter
```{r status_messages, tidy=TRUE, echo=FALSE}

  message("R markdwon report processed the last time: ", Sys.time())  
  message("The following analysis is performed based on that kind of data: ", params$quant)
  message("Analysis has been performed in the following manner: ", params$analysis)
  
  message("Internal standard values have been limited to: ", params$ca_fac)
  
  message("Generated output-files saved in: ", paste0(path_setup, set_output, set_pp))
  
  par = as.data.frame(do.call(rbind, params))
  colnames(par) = c("Value")
  knitr::kable(par, caption = "Parametrized report values", format = "markdown")
  
```

# Post-processing of the data
## Project settings

```{r init_project, tidy=TRUE}

  setup_params = read.csv(paste0(path_setup,"MTXQC_params.csv"), TRUE)

  # kind of experiment
  idx_exp = as.character(setup_params[which(setup_params$Parameter == "data"), "Value"])
  
  #### Annotation file
    ann_idx = as.character(setup_params[which(setup_params$Parameter == "ann"), "Value"])
    ann  =  read.csv(paste0(path_setup, set_input, ann_idx), T)
    
    if (ncol(ann) == 1) {
          ann = read.csv(paste0(path_setup, set_input, ann_idx), T, sep = ";")
          message("colon separator detected: ", ann_idx)
    }
   
  #### Parameter for stats
  parvec = check_parvector(params)
  
  ### (1) Calculated absolute quantities derived from MTXQCp1
    if ((params$updated == "PeakArea") & (params$updated == "both")) {
      data_q = read.csv(paste0(path_setup, set_output,'quant/CalculationFileData_manVal.csv'), T) 
    } else {
      data_q = read.csv(paste0(path_setup, set_output,'quant/CalculationFileData.csv'), T) 
    }

  
  ### (2) Isotope incorporation (export derived from MAUI or Metmax)
  
    if (idx_exp != "qMTX") {  
      if (params$updated == "Incorporation") {
        data_inc = read.csv(paste0(path_setup, set_input, "inc/Incorporation_manVal.csv"), T)
      } else {
        data_inc = read.csv(paste0(path_setup, set_input, "inc/DataMatrix.csv"), T)
      }
    } else {
      message("Current experiment does not include application of stable isotopes.")
    }

  ### Defined experimental conditions according to annotation file
  knitr::kable(summary(ann),caption = "Statistical summary of annotation file.", format = "pandoc")

```

## Define your color scheme and order of factor levels

This report provides a number of pre-defined colors for an graphical output regarding your personal taste. 

Nevertheless - there is one drawback! If due to missing values for a certain subgroup of your data leads to
an complete absence of a subgroup, the same color might be associated with different levels. 

Therefore I still recommend to optimize the color style in relation to your project in the chunk below. The fixed definition 
of color code and factor level prevents the mis-association of different factor levels with the same color.



```{r color_definitions}
#MOD!
# This are only examples how you could adjust your plots regarding your needs

# Color code (example: )
# Replace names of colors with levels of a factor (entries of one condition)
# Color codes can be either hex-codes or names, e.g., tomato3
# good resource: colorbrewer2.org

# color_cl <- c(
#   `M0` = "#d53e4f",
#   `M0HS` = "#fc8d59",
#   `M1` = "#fee08b",
#   `M1HS` = "#ffffbf",
#   `M2` = "#e8f598",
#   `M2HS` = "#99d594",
#   `Itaconate` = "#3288bd"
# )

# Order of factor levels in plots corresponds to the order of levels
# NA in graphs are a result of an typo in level names

# data_q$Substrate = factor(data_q$Substrate, levels = c("13C15N-Q", "1,5-Cit", "Acetate"))
# ann$Substrate = factor(ann$Substrate, levels = c("13C15N-Q", "1,5-Cit", "Acetate"))


```


# Absolute quantification of metabolites

```{r Quant_parameter, tidy = TRUE, include = TRUE, echo = TRUE}

  #1 - QC cinnamic acid / internal standard == within
  if (params$ca_fac == "QC") {
    data_q_sel = subset(data_q, data_q$IntStd_eval == 'within')
  } else {
    data_q_sel = data_q
  }
     
  
  #2 - Linear range calibration curve
  if (params$analysis == "stringent") {
    data_q_sel = subset(data_q_sel, data_q_sel$islinear == 'linear')
  } else {
    data_q_sel = subset(data_q_sel, data_q_sel$calc_curve == "yes_cal")
  }
  
  write.csv(data_q_sel, paste0(path_setup, set_output, set_pp, 
                               "Input_values_quantification.csv"), row.names = FALSE)
  
```  


# Statistics and data visualization
## Perform the required statistics

```{r statistics, tidy=TRUE}
#MOD!
    comb_par = c("Lettercode", parvec)
    print(comb_par)
    message("Statistics have been performed based on this grouping variables!")

    quant_ss_plain = ddply(data_q_sel, comb_par, summarise, 
                 n_Q = length(get(params$quant)),
                 mean_Q = mean(get(params$quant), na.rm = T),
                 sd_Q = sd(get(params$quant), na.rm = T)) 

  #Export 1
  #plain statistics
   write.csv(quant_ss_plain, paste0(path_setup, set_output, set_pp, "Statistics_absQ_", 
                             params$quant, "_plain.csv"), row.names = F)
   
   if (nrow(quant_ss_plain) == 0) {
     message("WARNING: Empty file detected: Statistics_absQ_par_plain.csv")
   }
   
```


## Advanced plotting of quantities for a complex set of conditions

```{r quantities_plot, echo=FALSE, fig.width=6, fig.height=3}
#MOD!  

if (length(parvec) == 2) {

      for (var in unique(quant_ss_plain$Lettercode)) {
        print( 
          ggplot(subset(quant_ss_plain, quant_ss_plain$Lettercode == var), 
                 aes(factor(get(params$par1)), mean_Q, fill = get(params$par2))) +
                  geom_bar(stat = "identity", size = 0.3, color = "black", 
                            position = position_dodge2(0.9)) +
                  geom_errorbar(aes(ymin = mean_Q - sd_Q, ymax = mean_Q + sd_Q), 
                              width = 0.25, position = position_dodge2(0.9)) +
                  theme_bw() +
                  facet_wrap(~ get(params$par2), scales = "free_y", nrow = 1) +
                  ggtitle(paste0(var)) +
                  xlab("Labeling time in (min)") +
                  ylab("Norm. Quantity in (pmol/1e+6 cells)") +
                  scale_fill_mtxqc(palette = mypal) +
                  #scale_fill_viridis_d() +
                  theme(axis.text.y = element_text(size = 8),
                    axis.text.x = element_text(size = 8, angle = 90, vjust = 0.3),
                    #axis.text.x = element_blank(),
                    strip.background = element_rect(fill = "white"),
                    strip.text.x = element_text(size = 8),
                    strip.text.y = element_text(size = 8)) +
                  theme(legend.position = "bottom", legend.title = element_blank(), 
                          legend.spacing.x = unit(0.25,"cm")) +
                  guides(fill = guide_legend(nrow = 1))
        )
      }
}


if (length(parvec) == 3) {

      for (var in unique(quant_ss_plain$Lettercode)) {
        print( 
          ggplot(subset(quant_ss_plain, quant_ss_plain$Lettercode == var), 
                 aes(factor(get(params$par1)), mean_Q, fill = get(params$par2))) +
                  geom_bar(stat = "identity", size = 0.3, color = "black", 
                            position = position_dodge2(0.9)) +
                  geom_errorbar(aes(ymin = mean_Q - sd_Q, ymax = mean_Q + sd_Q), 
                              width = 0.25, position = position_dodge2(0.9)) +
                  theme_bw() +
                  facet_wrap(get(params$par3) ~ get(params$par2), scales = "free_y") +
                    ggtitle(paste0(var)) +
                  xlab("Labeling time in (min)") +
                  ylab("Norm. Quantity in \n (pmol/1e+6 cells)") +
                  scale_fill_mtxqc(mypal) +
                  theme(axis.text.y = element_text(size = 8),
                    axis.text.x = element_text(size = 8, angle = 90, vjust = 0.3),
                    #axis.text.x = element_blank(),
                    strip.background = element_rect(fill = "white"),
                    strip.text.x = element_text(size = 8),
                    strip.text.y = element_text(size = 8)) +
                   theme(legend.position = "bottom", legend.spacing.x = unit(0.25, "cm")) +
                  guides(fill = guide_legend(nrow = 1))
        )
      }
}

if (length(parvec) == 4) {

      for (var in unique(quant_ss_plain$Lettercode)) {
        print( 
          ggplot(subset(quant_ss_plain, quant_ss_plain$Lettercode == var), 
                 aes(factor(get(params$par1)), mean_Q, fill = get(params$par2))) +
                  geom_bar(stat = "identity", size = 0.3, color = "black", 
                            position = position_dodge2(0.9)) +
                  geom_errorbar(aes(ymin = mean_Q - sd_Q, ymax = mean_Q + sd_Q), 
                              width = 0.25, position = position_dodge2(0.9)) +
                  theme_bw() +
                  facet_wrap(get(params$par4) ~ get(params$par3), scales = "free_y") +
                  ggtitle(paste0(var)) +
                  xlab("Labeling time in (min)") +
                  ylab("Norm. Quantity in \n (pmol/1e+6 cells)") +
                      scale_fill_mtxqc(mypal) +
                  theme(axis.text.y = element_text(size = 8),
                    axis.text.x = element_text(size = 8, angle = 90, vjust = 0.3),
                    #axis.text.x = element_blank(),
                    strip.background = element_rect(fill = "white"),
                    strip.text.x = element_text(size = 8),
                    strip.text.y = element_text(size = 8)) +
                    theme(legend.position = "bottom", legend.title = element_blank(),
                          legend.spacing.x = unit(0.25, "cm")) + 
                guides(fill = guide_legend(nrow = 1))
        )
      }
}

``` 


# Metabolic profile

An automated plot of the metabolic profile has been implemented so far only for statistical
analysis with *two* parameter!

In case of more parameter the metabolic profile is determined and exported, but not
shown with a graphical output. Still under constructions.


```{r qq_profile, results="asis", echo=FALSE, tidy=TRUE}
#MOD! CHOOSE THE CORRECT STATISTICS

  #add annotation
  quant_profile = merge(quant_ss_plain, pathway_profile)

  #rename mean and sd value for stationary pool
  colnames(quant_profile)[grepl("mean_Q", colnames(quant_profile))] <- "Q_ss"
  colnames(quant_profile)[grepl("sd_Q", colnames(quant_profile))] <- "sd_Q_ss"

  #Function call - Metabolic profile
  metabolic_profile = metabolic_profile(quant_profile, parvec, "Q_ss")
  
```  
  


# Relative quantification (no calibration curve)

```{r rel_quantities}
data_nocalc = subset(data_q, data_q$calc_curve == "no_cal")
un_nocalc = length(unique(data_nocalc$Lettercode))

message("Number of metabolites with peak area for relative analysis only: ", 
          un_nocalc)
```


```{r stats_relareas, tidy = TRUE}
data_nocalc_stat = ddply(data_nocalc, comb_par, summarise,
                         n_area = length(PeakArea),
                         mean_area = mean(PeakArea),
                         sd_area = sd(PeakArea))

#Export
write.csv(data_nocalc_stat, paste0(path_setup, set_output, set_pp, 
                                   "Statistics_peakArea_NoCal.csv"), row.names = F)
```

```{r relative_areas, echo=FALSE, fig.height=3, fig.width=6}
#MOD!

  if (length(parvec) == 2 ) {
    for (var in unique(data_nocalc_stat$Lettercode)) {
     print( 
       ggplot(subset(data_nocalc_stat, data_nocalc_stat$Lettercode == var), 
                    aes(factor(get(params$par1)), mean_area, fill = get(params$par2))) +
        geom_bar(stat = "identity", size = 0.3, color = "black", 
                      position = position_dodge2(0.9)) +
        geom_errorbar(aes(ymin = mean_area - sd_area, ymax = mean_area + sd_area), 
                      width = 0.25, position = position_dodge(0.9)) +
        facet_wrap( ~ get(params$par2), scales = "free_y", nrow = 1) +
        ggtitle(paste0(var)) +
        ylab("Peak Area in (-)") +
        xlab("Days") +
        theme_bw() +
            scale_fill_mtxqc(mypal) +
        theme(axis.text.y = element_text(size = 8),
              axis.text.x = element_text(size = 8, angle = 90, vjust = 0.3),
              strip.background = element_rect(fill = "white"),
              strip.text.x = element_text(size = 8),
              strip.text.y = element_text(size = 8)) +
        theme(legend.position = "bottom", legend.spacing.x = unit(0.25, "cm"), legend.title = element_blank()) +
         guides(fill = guide_legend(nrow = 1))
      )
    }
  }

  if (length(parvec) == 3 ) {
    for (var in unique(data_nocalc_stat$Lettercode)) {
     print( 
       ggplot(subset(data_nocalc_stat, data_nocalc_stat$Lettercode == var), 
                    aes(factor(get(params$par1)), mean_area, fill = get(params$par2))) +
        geom_bar(stat = "identity", size = 0.3, color = "black", 
                      position = position_dodge2(0.9)) +
        geom_errorbar(aes(ymin = mean_area - sd_area, ymax = mean_area + sd_area), 
                      width = 0.25, position = position_dodge(0.9)) +
        facet_wrap( ~ get(params$par3), scales = "free_y", nrow = 1) +
        ggtitle(paste0(var)) +
        ylab("Peak Area in (-)") +
        xlab("Days") +
        theme_bw() +
            scale_fill_mtxqc(mypal) +
        theme(axis.text.y = element_text(size = 8),
              axis.text.x = element_text(size = 8, angle = 90, vjust = 0.3),
              strip.background = element_rect(fill = "white"),
              strip.text.x = element_text(size = 8),
              strip.text.y = element_text(size = 8)) +
       theme(legend.position = "bottom", legend.spacing.x = unit(0.25, "cm"), legend.title = element_blank()) +
         guides(fill = guide_legend(nrow = 1))
      )
    }
  }

 if (length(parvec) == 4 ) {
    for (var in unique(data_nocalc_stat$Lettercode)) {
     print( 
       ggplot(subset(data_nocalc_stat, data_nocalc_stat$Lettercode == var), 
                    aes(factor(get(params$par1)), mean_area, fill = get(params$par2))) +
        geom_bar(stat = "identity", size = 0.3, color = "black", 
                      position = position_dodge2(0.9)) +
        geom_errorbar(aes(ymin = mean_area - sd_area, ymax = mean_area + sd_area), 
                      width = 0.25, position = position_dodge(0.9)) +
        facet_wrap(get(params$par4) ~ get(params$par3), scales = "free_y") +
        ggtitle(paste0(var)) +
        ylab("Peak Area in (-)") +
        xlab("Days") +
        theme_bw() +
            scale_fill_mtxqc(mypal) +
        theme(axis.text.y = element_text(size = 8),
              axis.text.x = element_text(size = 8, angle = 90, vjust = 0.3),
              strip.background = element_rect(fill = "white"),
              strip.text.x = element_text(size = 8),
              strip.text.y = element_text(size = 8)) +
       theme(legend.position = "bottom", legend.spacing.x = unit(0.25, "cm"), legend.title = element_blank()) +
         guides(fill = guide_legend(nrow = 1))
      )
    }
  }
    
```


```{r check_qtx}
  inc_idx = as.character(setup_params[which(setup_params$Parameter == "data"), "Value"])

  if (inc_idx == "qMTX") {
    message("The herein project does not include stable isotope labeling data!")
    knitr::knit_exit()
  } 

```


# Stable isotope incorporation
## Statistics and bar plots

```{r inc_statistics}

  if ("QuantMasses" %in% colnames(data_inc)) {
        data_inc_m = melt(data_inc, id.vars = c('Metabolite', 'QuantMasses'), variable.name = c('File'), value.name = c('Inc'))
  } else {
        data_inc_m = melt(data_inc, id.vars = c('Metabolite'), variable.name =  c('File'), value.name = c('Inc'))
  }

  #merge with conversion and annotation file
  data_inc_se = merge(data_inc_m, con_se)
  data_inc_se = merge(data_inc_se, ann)
  
  #add Batch-Id
  data_inc_se$Batch_Id = sapply(strsplit(as.character(data_inc_se$File), "\\_"), "[[", 1)
  
  #replace negative values with NA
  data_inc_se = within(data_inc_se, Inc[Inc < 0.000] <- 0)

```

```{r stats_tr_inc, tidy=TRUE}
#MOD!

  comb_li = c("Lettercode", "Mass_Pos", parvec)
  message("Parameter used for statistics: 13C-Incorporation")  
  print(comb_li)

  #statistics technical repl
  li_stat  =  ddply(data_inc_se, comb_li, summarise,
                    		N_inc  =  length(Inc),
                    		mean_li  =  mean(Inc, na.rm = T),
                    		sd_li  =  sd(Inc, na.rm = T))

  #Export
  write.csv(li_stat, paste0(path_setup, set_output, set_pp ,'13C_statistics.csv'), 
            row.names = F)
```


```{r inc_plot, fig.width=6, fig.height=3}
#MOD!

#Exclude additional fragments for plotting procedure (Cit: m/z 375, 2OG: m/z 288)

li_stat = subset(li_stat, li_stat$Mass_Pos <= 270)


if (length(parvec) == 2) {   
    for (var in unique(li_stat$Lettercode)) {
      print(
        ggplot(subset(li_stat, li_stat$Lettercode == var), 
                      aes(factor(get(params$par1)), mean_li, fill = get(params$par2))) +
            geom_bar(stat = "identity", size = 0.3, color = "black", 
                      position = position_dodge(0.9)) +
            geom_errorbar(aes(ymin = mean_li - sd_li, ymax = mean_li + sd_li),
                      width = 0.25, position = position_dodge(0.9)) +
            facet_wrap( ~ Mass_Pos, scales = "free_y") +
            theme_bw() +
            ggtitle(var) +
                scale_fill_mtxqc(mypal) +
            ylab("13C-Inc. in (-)") +
            xlab(" ") +
            theme(axis.text.y = element_text(size = 8),
                  axis.text.x = element_text(size = 8),
                  strip.background = element_rect(fill = "white"),
                  strip.text.x = element_text(size = 8),
                  strip.text.y = element_text(size = 8)) + 
            theme(legend.position = "bottom", legend.title = element_blank(), legend.spacing.x = unit(0.25, "cm")) +
          guides(fill = guide_legend(nrow = 1))
      )
    }
} 

if (length(parvec) == 3) {   
    for (var in unique(li_stat$Lettercode)) {
      print(
        ggplot(subset(li_stat, li_stat$Lettercode == var), 
                      aes(factor(get(params$par1)), mean_li, fill = get(params$par2))) +
            geom_bar(stat = "identity", size = 0.3, color = "black", 
                      position = position_dodge(0.9)) +
            geom_errorbar(aes(ymin = mean_li - sd_li, ymax = mean_li + sd_li),
                      width = 0.25, position = position_dodge(0.9)) +
            facet_wrap(get(params$par3) ~ Mass_Pos, scales = "free_y") +
            theme_bw() +
            ggtitle(var) +
                scale_fill_mtxqc(mypal) +
            ylab("13C-Inc. in (-)") +
            xlab(" ") +
            theme(axis.text.y = element_text(size = 8),
                  axis.text.x = element_text(size = 8),
                  strip.background = element_rect(fill = "white"),
                  strip.text.x = element_text(size = 8),
                  strip.text.y = element_text(size = 8)) + 
            theme(legend.position = "bottom", legend.title = element_blank(), legend.spacing.x = unit(0.25, "cm")) +
          guides(fill = guide_legend(nrow = 1))
      )
    }
} 




```

# 13C-Quantities

*Previously:* This part mulitplies the levels derived from the metabolic profile with the incorporation rate of a defined time point and visualizes the 13C-labeled quantities of each intermediate.

*Corrected:* Absolute quantities are applied in substrate-specific manner due to the variation of pool sizes.

```{r LabQ_calc}

   #Absolute quantities / Substrate
   stats_q = read.csv(paste0(path_setup, set_output, set_pp,
                       "Statistics_absQ_", params$quant, "_plain.csv"), TRUE)
  
   #13C-Inc statistics
   stats_inc = read.csv(paste0(path_setup, set_output, set_pp,
                       "13C_statistics.csv"), T)
  
   #annotations
   data_combined = merge(stats_q, stats_inc)
  
   #calculate 13C-quantities
   data_combined$mean_labQ = data_combined$mean_li * data_combined$mean_Q
   data_combined$sd_labQ = data_combined$sd_li * data_combined$mean_Q
  
   #Export
   write.csv(data_combined, paste0(path_setup, set_output, set_pp,
             "Statistics_13ClabQuantities_", params$quant, ".csv"), row.names = F)
```


```{r LabQ, echo=FALSE, fig.width=6, fig.height=3}

if (length(parvec) == 2) {
   for (var in unique(data_combined$Lettercode)) {
     print(
       ggplot(subset(data_combined, data_combined$Lettercode == var),
                      aes(factor(get(params$par1)), mean_labQ, fill = get(params$par2))) +
       geom_bar(stat = "identity", size = 0.3, color = "black",
                    position = position_dodge(0.9)) +
       geom_errorbar(aes(ymin = mean_labQ - sd_labQ, ymax = mean_labQ + sd_labQ),
                    width = 0.25, position = position_dodge(0.9)) +
          facet_wrap(~ Mass_Pos, scales = "free_y") +
          theme_bw() +
        	    scale_fill_mtxqc(mypal) +
          ggtitle(paste0(var)) +
          ylab("13C-Quantities in (pmol/mg tissue)") +
         	theme(axis.text.y = element_text(size = 8),
                axis.text.x = element_text(size = 8),
                #axis.text.x = element_blank(),
                strip.background = element_rect(fill = "white"),
                strip.text.x = element_text(size = 8),
                strip.text.y = element_text(size = 8),
                legend.title = element_blank(), legend.spacing.x = unit(0.25, "cm"))
    )
   }
}


if (length(parvec) == 3) {
   for (var in unique(data_combined$Lettercode)) {
     print(
       ggplot(subset(data_combined, data_combined$Lettercode == var),
                      aes(factor(get(params$par1)), mean_labQ, fill = get(params$par2))) +
       geom_bar(stat = "identity", size = 0.3, color = "black",
                    position = position_dodge(0.9)) +
       geom_errorbar(aes(ymin = mean_labQ - sd_labQ, ymax = mean_labQ + sd_labQ),
                    width = 0.25, position = position_dodge(0.9)) +
          facet_wrap(get(params$par3) ~ Mass_Pos, scales = "free_y") +
          theme_bw() +
        	    scale_fill_mtxqc(mypal) +
          ggtitle(paste0(var)) +
          ylab("13C-Quantities in (pmol/mg tissue)") +
         	theme(axis.text.y = element_text(size = 8),
                axis.text.x = element_text(size = 8),
                #axis.text.x = element_blank(),
                strip.background = element_rect(fill = "white"),
                strip.text.x = element_text(size = 8),
                strip.text.y = element_text(size = 8),
                legend.title = element_blank(), legend.spacing.x = unit(0.25, "cm"))
    )
   }
} 

```

\clearpage

# List of abbreviations

```{r abbrev, echo=FALSE}

library(tidyverse)
temp = select(con_se, "Lettercode", "Metabolite_short", "Q1_value" ,"nopsirm") %>%
      arrange(Lettercode)

kable(temp,"latex", longtable = T, booktabs = T) %>%
  add_header_above(c("")) %>%
  kable_styling(latex_options = c("scale_down", "repeat_header"))
```

*End of the document!*