---
title: "A Mass Spec Data set for Targeted Proteomics Performance Evaluation"
author:
- name: Tobias Kockmann
  affiliation:
    - &fgcz Functional Genomics Center Zurich - Swiss Federal Institute of Technology in Zurich
- name: Christian Trachsel
  affiliation:
    - *fgcz
- name: Christian Panse
  email: cp@fgcz.ethz.ch
  affiliation:
    - *fgcz
    - &sib Swiss Institute of Bioinformatics, Quartier Sorge - Batiment Amphipole, CH-1015 Lausanne, Switzerland
- name: Jonas Grossmann
  affiliation:
    - *fgcz
    - *sib
- name: Witold E. Wolski
  affiliation:
    - *fgcz
    - *sib
abstract: |
  This document reproduces the Figures presented in [-@msqc1]. For a
  description of the theory behind applications shown here we refer to the
  original manuscript. The results differ slightly due to technical
  changes or bugfixes in `r BiocStyle::Biocpkg('msqc1')` that have
  been implemented after the manuscript was printed.
package: "`r pkg_ver('msqc1')`"
date: "`r Sys.Date()`"
bibliography:
  - msqc1.bib
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Introducing msqc1 - A Mass Spec Data set for Targeted Proteomics Performance Evaluation}
  %\VignetteEncoding{UTF-8}
output:
  BiocStyle::html_document:
    toc_float: true
---

<!--
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{msqc1 - A Mass Spec Data set for Targeted Proteomics Performance Evaluation}
-->


```{r echo = FALSE, message = FALSE}
if (!require(msqc1)){
  ## try http:// if https:// URLs are not supported
  if (!requireNamespace("BiocManager", quietly=TRUE))
      install.packages("BiocManager")
}
```

# Introduction 

Here we provide an experiment data package containing all the spread sheets used for figures and supplemental figures published in @msqc1.
The raw instrument files are accessible for registered users through 
[https://fgcz-bfabric.uzh.ch](https://fgcz-bfabric.uzh.ch) [-@panse2022bridging]
as project 1959. Raw LC-MS data from all five platforms were imported into
[Skyline 3.1](https://skyline.gs.washington.edu), see @pmid20147306, processed, 
and are available though 
[https://panoramaweb.org/labkey/MSQC1.url](https://panoramaweb.org/labkey/MSQC1.url) when published.
This package contains data.frame objects exported by Skyline making the data available for the R world. 

# MSQC1 Peptide properties

## Figure - reference L:H ratio versus the on-column amount

```{r "Supplementary Figure1", echo=FALSE, fig.width=6, fig.height=6, fig.retina=3, warning=FALSE, message = FALSE}
plot(10 * msqc1::peptides$SIL.peptide.per.vial, msqc1::peptides$actual.LH.ratio,
     log='xy',
     ylab='reference L:H ratio',
     xlab='on column amount [fmol]',
     axes=FALSE,
     xlim=c(0.8, 2000));

axis(1, 10 * peptides$SIL.peptide.per.vial, 10 * peptides$SIL.peptide.per.vial )
axis(2, peptides$actual.LH.ratio, peptides$actual.LH.ratio)
text(10 * peptides$SIL.peptide.per.vial, peptides$actual.LH.ratio, peptides$Peptide.Sequence, cex=0.5 ,pos=4, srt=11)
box()
```

The scatter-plot displays the reference Light:Heavy ratio versus the on-column amount of heavy peptide of the MSQC1 peptide mix. Note, x and y axis are drawn in log scale.

# Data

show the dilution series
```{r}
table(msqc1_dil$relative.amount)
```

show the peptide frequency
```{r}
table(msqc1_dil$Peptide.Sequence)
```

show ion types
```{r}
table(msqc1_dil$Fragment.Ion)
```


show instruments
```{r}
table(msqc1_dil$instrument)
```

# Illustrations

## Figure - Sigma mix peptide level signals
```{r "Figure1", echo=FALSE, fig.width=14, fig.height=9, fig.retina=3, warning=FALSE, message = FALSE}
msqc1:::.figure1(data=msqc1_8rep, peptides=peptides)
```

**Sigma mix peptide level signals** - The graph displays the log2 L:H area ratios of eight technical replicates over 14 peptides from five MS platforms. The 14 panels were ordered by the on column amount of heavy peptide per vial (0.8, 4, 20, 40, 80, 200, 500, 1000 fmol). The black line indicates the theoretical L:H ratio as reported in the QC certificate by Sigma-Aldrich. In a perfect setting all data points would be located close to the black line indicating a perfect match between theoretical and measured log2 L:H ratios.
The dark gray boxes correspond to a 0.5 and the light grey boxes to a deviation of 1 from the expected value (black line).


## Figure - Volcano Plot


```{r, echo=FALSE}
# taken from WEW
.compute_volcano_tuple <- function(x, y, adjust=TRUE, alternative="two.sided"){
  stopifnot(nrow(x) == nrow(y))

  res <- lapply(1:nrow(x), function(i){
    r <- t.test(x[i,], y[i,], paired =FALSE, alternative=alternative)
    data.frame(FC=(r$estimate[1] - r$estimate[2]), p.value=r$p.value)
  })
  
  res <- do.call('rbind', res)
  
  if(adjust){
    res$pval <- p.adjust(res$p.value, method="BH")
  }
  return(res)
} 

.shape_for_volcano <- function(S, peptides){
  
  S <- S[grep("[by]", S$Fragment.Ion), ]
  S <- S[S$Peptide.Sequence %in% peptides$Peptide.Sequence, ]
  S <- aggregate(Area ~ Peptide.Sequence + Isotope.Label.Type + instrument + File.Name.Id,
                data=S, FUN=sum)
  
  
  S <- data.frame(
    Peptide.Sequence=paste(S$Peptide.Sequence, S$Isotope.Label.Type, sep=','),
    Area=S$Area,
    instrument=paste(S$instrument), rep=S$File.Name.Id)
  
  
  S <- reshape(S,
               v.name="Area",
               idvar=c("Peptide.Sequence", 'instrument'),
               timevar="rep",
               direction = "wide")
  
  p <- gsub(".light", "", gsub(".heavy", "", S$Peptide.Sequence))

  row.names(S) <- 1:nrow(S)
  S <- (droplevels(S[which(p %in% names(table(p)[table(p) == 10])), ]))
  S <- S[order(paste(S$instrument, S$Peptide.Sequence )), ]
  
  idx.l <- grep('light', S$Peptide.Sequence )
  idx.h <- grep('heavy', S$Peptide.Sequence )
  
  L<-S[idx.l, 3:10]
  H<-S[idx.h, 3:10]
  v <- .compute_volcano_tuple(log2(L), log2(H))
  
  Peptide.Sequence <- gsub(".light", "", S$Peptide.Sequence[idx.l])
  instrument <- S$instrument[idx.l]
  
  S <- data.frame(Peptide.Sequence, instrument, log2FC=v$FC, p.value=v$p.value)
}

```

```{r fig.width=5, fig.height=10, fig.retina=3, echo=TRUE}
S <- .shape_for_volcano(S=msqc1_8rep, peptides)

msqc1:::.figure_setup()

xyplot(-log(p.value, 10) ~ log2FC | instrument, data=S, group=Peptide.Sequence,
       panel = function(...){
         panel.abline(v=c(-1,1), col='lightgray')
         panel.abline(h=-log(0.05,10), col='lightgray')
         panel.xyplot(...)
       },
       ylab='-log10 of p-value',
       xlab='log2 fold change',
       layout=c(1,5),
       auto.key = list(space = "right", points = TRUE, lines = FALSE, cex=1))

```


## Figure - Ratio stability upon analyte dilution

```{r "Figure3a", echo=FALSE, fig.width=14, fig.height=9, fig.retina=3, warning=FALSE, message = FALSE}
msqc1:::.figure4(data=msqc1_dil, peptides=peptides)
```

**Ratio stability upon analyte dilution** - Each scatter-plot panel displays the experimental derived log2 L:H ratios versus the relative amount. The panels are ordered by the SIL on column amount (lower left to upper right).
Color grouping was done by instrument. The loess fit curves were added for visualizing the trend.
The SIL value given in each panel legend is valid for the relative amount of 1.
The horizontal black line indicates the theoretical log2 L:H ratio.

## Figure - Accuracy

```{r "Figure3b", echo = FALSE, fig.width = 14, fig.height = 5, fig.retina = 3, warning = FALSE, message = FALSE}
msqc1:::.figure5(data=msqc1_dil, data_reference = msqc1_8rep)
```

**Accuracy** - The graph displays in each panel a sensitivity curves for a given relative amount.



# Session information
```{r}
sessionInfo()
```

# References