Giai Gianetto Quentin
email: quentin.giaigianetto@pasteur.fr

Statistical analysis of post-translational modifications quantified by label-free proteomics across multiple biological conditions with R: illustration from SARS-CoV-2 infected cells

Keywords: Statistics, R, Data quality control, Clustering, Post-translational modifications, Label-free proteomics, Relative quantification

Protein post-translational modifications (PTMs) are essential elements of cellular communication. Their variations in abundance can affect cellular pathways, leading to cellular disorders and diseases. A widely used method for revealing PTM-mediated regulatory networks is their label-free quantitation (LFQ) by high-resolution mass spectrometry. The raw data resulting from such experiments are generally interpreted using specific software, such as MaxQuant, MassChroQ or Proline for instance. They provide data matrices containing quantified intensities for each modified peptide identified. Statistical analyses are then necessary (1) to ensure that the quantified data are of good enough quality and sufficiently reproducible, (2) to highlight the modified peptides that are differentially abundant between the biological conditions under study. The objective of this chapter is therefore to provide a complete data analysis pipeline for analyzing the quantified intensities of modified peptides in presence of two or more biological conditions using the R software. We illustrate our pipeline starting from MaxQuant outputs dealing with the analysis of A549-ACE2 cells infected by SARS-CoV-2 at different time stamps, freely available on PRIDE (PXD020019).

Introduction

Mass spectrometry (MS)-based proteomics allow the identification and quantification of a large number of post-translational protein modifications [START_REF] Witze | Mapping protein post-translational modifications with mass spectrometry[END_REF] [START_REF] Zhao | Modification-specific proteomics: strategies for characterization of post-translational modifications using enrichment techniques[END_REF]. A general workflow consists of (i) a proteolytic digestion of proteins into peptides, (ii) an enrichment step to increase the concentration of modified peptides in the samples, and (iii) the analysis of the samples by liquid chromatography-tandem mass spectrometry (LC-MS/MS) [START_REF] Pieroni | Enrichments of post-translational modifications in proteomic studies[END_REF]. It is often used to study protein phosphorylations [START_REF] Ficarro | Magnetic bead processor for rapid evaluation and optimization of parameters for phosphopeptide enrichment[END_REF] [START_REF] Pinkse | A protocol on the use of titanium dioxide chromatography for phosphoproteomics[END_REF], which is the most commonly studied post-translational modification. This workflow can also be adapted to study other post-translational modifications such as 1 This is a pre-print of the following chapter: Quentin Giai Gianetto, Statistical Analysis of Post-Translational Modifications Quantified by Label-Free Proteomics Across Multiple Biological Conditions with R: Illustration from SARS-CoV-2 Infected Cells published in Statistical analysis of Proteomic Data, edited by Thomas Burger, 2023, Humana, New York, NY reproduced with permission of Springer Nature. The final authenticated version is available online at: https://doi.org/10.1007/978-1-0716-1967-4_12 ubiquitinations [START_REF] Udeshi | Refined preparation and use of anti-diglycine remnant (k-𝜀-gg) antibody enables routine quantification of 10,000 s of ubiquitination sites in single proteomics experiments[END_REF], methylations [START_REF] Carlson | Proteome-wide enrichment of proteins modified by lysine methylation[END_REF], acetylations [START_REF] Kim | Substrate and functional diversity of lysine acetylation revealed by a proteomics survey[END_REF], or several kind of modifications in unified pipelines [START_REF] Mertins | Integrated proteomic analysis of post-translational modifications by serial enrichment[END_REF].

Label-free quantification enables large-scale analyzes and can be applied to experiments composed of many samples and biological conditions, as it can be the case in the field of clinical screening for instance. It allows avoiding drawbacks of labeling methods, whether metabolic (SILAC) or chemical (iTRAQ), which are limited by the cost of labeling reagents, some ineffective labeling, or the limited number of samples that can be analyzed. However, label-free quantification requires careful experimental designs to achieve good reproducibility of quantified values. Indeed, the enrichment step introduces additional variations to the ones traditionally produced by the fluctuation of the liquid chromatography and the ionization conditions in MS-based proteomics. This is the reason why a special care must be taken during this step to carry it out as rigorously as possible in order to maximize the reproducibility of experiments.

The MS/MS spectra obtained after the LC-MS/MS analysis have to be associated to peptides. The identification of spectra is generally carried out by searching the measured MS/MS spectra in a theoretical database, freely downloadable from UniProt website for instance (https://www.uniprot.org/). The localization of a modification on a specific aminoacid of the identified peptide is assessed using a score or a probability, depending on the method applied. It is usually calculated by comparing the measured spectrum with theoretical spectra where the modification is placed on each possible amino acid of the considered peptide [START_REF] Chalkley | Modification site localization scoring: strategies and performance[END_REF]. Currently, this research is generally done thanks to specific software, such as MaxQuant [START_REF] Tyanova | The maxquant computational platform for mass spectrometry-based shotgun proteomics[END_REF], MassChroQ [START_REF] Valot | Masschroq: a versatile tool for mass spectrometry quantification[END_REF] or Proline [START_REF] Bouyssié | Proline: an efficient and user-friendly software suite for large-scale proteomics[END_REF], but can also be performed from R (see Note 1). Depending on the quantification method, either the MS spectra or the MS/MS spectra are next used to quantify the abundance of each modified peptides identified. In the end, such bioinformatics analyses result in data matrices composed of the intensities measured in each of the samples for all the modified peptides identified.

It is common to identify and quantify several tens of thousands peptides in usual experiments. To extract useful information from these large data sets, statistical analyzes of quantified data matrices are required, to check the reproducibility of the analyses and to determine the modified peptides of interest. The objective of this chapter is therefore to provide a data analysis pipeline for analyzing the quantified intensities of modified peptides, in presence of two or more biological conditions using the R software. This analysis pipeline is provided with examples of R codes adapted to a case study. The informed reader can easily change these codes to adapt them to his/her problematic and to what seems relevant to put forward in his/her analysis. We start from outputs of the MaxQuant software [START_REF] Tyanova | The maxquant computational platform for mass spectrometry-based shotgun proteomics[END_REF], which is one of the most popular software tool for this kind of analysis, but the proposed data analysis pipeline can easily be applied to the outputs of various software, as long as they provide the intensities of the modified peptides and the proteins they belong to.

Material

R and RStudio installation

(https://rstudio.com/products/rstudio/download/) which is a development environment to facilitate the development of R codes and packages. However, all the codes presented herein can be executed from a basic R installation.

Throughout this chapter we assume that the reader has some prior knowledge of R, and we give code examples without dealing with the basics of R, as for example the different types of R objects. For more details on R and its programming language, the reader may refer to the many courses available for free on the R and Bioconductor websites as well as many good books on R and Bioconductor such as [START_REF] Gentleman | R programming for bioinformatics[END_REF], [START_REF] Chambers | Software for data analysis: programming with R[END_REF], [START_REF] Millot | Comprendre et réaliser les tests statistiques à l'aide de r[END_REF].

R packages

In addition to the functions provided in the standard version of R, we will use functions from other R packages in the sample codes provided. These packages are the next: VennDiagram [START_REF] Chen | Venndiagram: a package for the generation of highlycustomizable venn and euler diagrams in r[END_REF], UpSetR [START_REF] Conway | Upsetr: an r package for the visualization of intersecting sets and their properties[END_REF], ggplot2 [START_REF] Wickham | ggplot2: elegant graphics for data analysis[END_REF], grid, ggplotify, jaccard, ggdendro [START_REF] De Vries | ggdendro: Create dendrograms and tree diagrams using 'ggplot2'. R package version[END_REF], ggridges [START_REF] Wilke | ggridges: Ridgeline plots in 'ggplot2'. R package version 053[END_REF], limma [START_REF] Ritchie | limma powers differential expression analyses for rna-sequencing and microarray studies[END_REF], cp4p [START_REF] Gianetto | Calibration plot for proteomics: A graphical tool to visually check the assumptions underlying fdr control in quantitative experiments[END_REF], ssize.fdr [START_REF] Liu | Quick calculation for sample size while controlling false discovery rate with application to microarray analysis[END_REF], imp4p [START_REF] Gianetto | A peptide-level multiple imputation strategy accounting for the different natures of missing values in proteomics data[END_REF], car [START_REF] Fox | car: Companion to applied regression[END_REF], multivariance [START_REF] Böttcher | Copula versions of distance multivariance and dhsic via the distributional transform-a general approach to construct invariant dependence measures[END_REF], factoextra [START_REF] Kassambara | factoextra: Extract and visualize the results of multivariate data analyses[END_REF], cluster, reshape2 [START_REF] Wickham | reshape2: Flexibly reshape data: a reboot of the reshape package[END_REF], ggpubr [START_REF] Kassambara | ggplot2" based publication ready plots[END_REF].

It is up to the reader to download these packages and to check that they are properly installed on his/her machine.

Data type

Throughout this chapter we assume that the data are outputs from software used to identify and quantify the abundance of modified peptides and their proteins, such as MaxQuant [START_REF] Tyanova | The maxquant computational platform for mass spectrometry-based shotgun proteomics[END_REF], MassChroQ [START_REF] Valot | Masschroq: a versatile tool for mass spectrometry quantification[END_REF] or Proline [START_REF] Bouyssié | Proline: an efficient and user-friendly software suite for large-scale proteomics[END_REF]. The quantitative data should fit into two matrices: the first contains quantified intensities for PTMs identified on peptides; the second contains quantified intensities for the non-modified proteins to which the modified peptides of the first matrix are associated (i.e. intensities deduced from non-modified peptides of associated proteins).

For instance, when using MaxQuant [START_REF] Tyanova | The maxquant computational platform for mass spectrometry-based shotgun proteomics[END_REF], the first dataset can be extracted from the "Phospho (STY)Sites.txt" (phosphorylation) or "GlyGly (K)Sites.txt" (ubiquitination), while the second dataset is in the "proteinGroups.txt" file.

The first matrix thus makes it possible to determine the dynamics of the abundance of a PTM between several biological conditions; while the second makes it possible to compare this dynamic with that of the unmodified protein across multiple biological conditions. This comparison is important to ensure that a difference in the abundance of a modified peptide is actually due to the abundance of the modification and not to the dynamic of the associated unmodified protein between compared conditions.

Case study dataset

To illustrate this protocol, we used a dataset that the reader can freely download from PRIDE (PXD020019). On this server, the reader will find MaxQuant outputs as well as the raw data from the mass spectrometers, such as he/she can reanalyze them using their own workflow and software of interest. This dataset corresponds to the analysis of phosphorylated and ubiquitinated peptides in A549-ACE2 cells infected by SARS-CoV-2 at different time stamps [START_REF] Stukalov | Multi-level proteomics reveals host-perturbation strategies of sars-cov-2 and sars-cov[END_REF]. Hereafter, we will focus on the analysis of phosphorylated peptides, but the same analyses can be performed from the ubiquitinated peptides. To reproduce the results presented along this protocol, the careful reader has to download the MaxQuant outputs "proteinGroups.txt" and "Phospho (STY)Sites.txt" from the PRIDE website.

Methods

Data preprocessing steps

1. Importing data in R: First, the software outputs used to identify and quantify the proteins as well as the peptides have to be saved into a format easily loadable in R. Generally, this can be performed by saving the software outputs into a "txt" or "csv" file. Then, these files can be imported in the R session into a data frame object (see Note 2). For instance, using MaxQuant outputs, the read.csv function can be used for this purpose: > data.prot=read.csv("(path_to_your_file)/proteinGroups.txt", header=TRUE,sep="\t",quote="") > data.ptm=read.csv("(path_to_your_file)/Phospho (STY)Sites.txt", header=TRUE,sep="\t",quote="") 2. Filtering peptides which are associated to "Reverse" or "Potential contaminant" peptides:

In MaxQuant, "reverse" peptides are artefactual peptides whose amino acid sequence corresponds to the one of peptides of the input database used for the identification, yet in reverse order (see Note 3). Additionally, "Potential contaminant" peptides are associated with proteins that commonly contaminate samples. Thus, these two kinds of identified peptides can be deleted before subsequent analysis. They are indicated by "+" in the reverse and potential contaminant columns: > data.ptm=data.ptm[which((data.ptm$Reverse!="+")& (data.ptm$Potential.contaminant!="+")),] 3. Filtering peptides which PTM location is not sufficiently well identified: In MaxQuant outputs, the maximum localization probability that has been estimated in one of the analyzed samples is in the "Localization.prob" column. In literature, it is common to keep only the modified peptides for which this probability is greater than 75% (see Experimental design corresponding to the columns of int.prot is defined by: > Cond.prot=factor(c(rep("mock",8),rep("SARS_COV_2",8)), levels=c("mock","SARS_COV_2")) > Time.prot=factor(rep(c(rep("24h",4), rep("6h",4)),2),levels= c("6h","24h")) > CondTime.prot=as.factor(paste(Cond.prot,Time.prot,sep=".")) Both experimental designs can be combined using: > Cond.ptmprot=unlist(list(Cond.ptm,Cond.prot)) > Time.ptmprot=unlist(list(Time.ptm,Time.prot)) Additionally, a vector indicating whether the combined design is for the peptide or protein will be useful in subsequent statistical analysis (see Subheading 3.9, Step 1): #Comp.ptmprot is equal to 1 if the coordinate is related #to the modified peptide and 0 otherwise > Comp.ptmprot=c(rep(1,length(Cond.ptm)), rep(0,length(Cond.prot)))

Checking the reproducibility of identifications

When the samples are highly reproducible, one can expect that the same modified peptides will be identified in each of the samples. The intersections between the sets of modified peptides identified in the samples can be visualized using Venn diagrams or UpSet graphs to highlight any aberrant sample. The following code displays the numbers of modified peptides found in common or specific of each SARS-Cov-2 samples at 6h (as illustrated in Fig. 1A):

.sample){ inter=intersect(peptid[,i],peptid[,j]); n2[i,j]=length(inter[inter!=0]); for (k in j:nb.sample){ inter=intersect(peptid[,i],intersect(peptid[,j],peptid[,k])); n3[i,j,k]=length(inter[inter!=0]); } }}
> require(VennDiagram) > venn.plot = draw.triple.venn( area1=n2[1,1], area2=n2[2,2], area3=n2[3,3], n12=n2[1,2], n13=n2[1,3], n23=n2[2,3], n123=n3[1,2,3]
, category = colnames(peptid), fill = c("dodgerblue", "goldenrod1", "seagreen3"), cat.col = c("dodgerblue", "goldenrod4", "seagreen4"), cat.cex = 1.5,cat.dist=0.1, margin = 0.1, cex = c(1.5,1.5,1.5,1.5,2,1.5,1.5), ind=TRUE); 4. Visualizing the intersections using UpSet graphs: Venn diagrams represent the sets of identified peptides using circles or ellipses. However, the larger the number of samples, the more difficult the Venn diagrams are readable. This is the reason why it may be preferable to use UpSet graphs (see Fig. 1B) when the number of samples becomes large (i.e. superior to > grid.newpage(); > p = grobTree(venn.plot); > venn.plot = as.ggplot(p); > upset.plot = as.ggplot(upset.plot); 6. Clustering samples using Jaccard index-based distances: The reproducibility of the identifications between two samples can be evaluated using Jaccard index. This index can be computed using the jaccard() function of the R package jaccard: > require(jaccard) > mat.jaccard=matrix(NA,ncol(data.upset),ncol(data.upset))

> for (i in 1:ncol(data.upset)){ for (j in 1:ncol(data.upset)){ mat.jaccard[i,j]=jaccard(data.upset[,i],data.upset[,j]); }} > colnames(mat.jaccard)=colnames(data.upset) > rownames(mat.jaccard)=colnames(data.upset) This matrix can be visualized using shades of blue (see Note 10). Additionally, a distance matrix based on this index can be computed with the as.dist() function, and a hierarchical clustering of samples can be performed with the hclust() function: > hv = hclust(as.dist(1-mat.jaccard),method="ward.D2") The clustering can be visualized with the ggdendrogram function of the R package ggdendro: > require(ggdendro) > ggd = ggdendrogram(hv, rotate = TRUE, theme_dendro = FALSE) > ggd = ggd + theme(axis.text.y = element_text(hjust = 1)) > ggd = ggd + xlab("Sample") > ggd = ggd + ylab("Height") > ggd = ggd + ggtitle("Ward's method with a \n Jaccard index-based distance") Additionally, the same approach can be used to evaluate the replication of the nonidentifications between samples (see Note 11).

Checking the reproducibility of quantified values

When the samples are highly reproducible, one can expect the quantified values of a same modified peptides to be similar across the biological replicates.

1. Log2-transforming the intensity values: Before analyzing quantified values, a log2 transformation is usually applied, in particular to decorrelate the intensity levels of peptides from their variances (see [START_REF] Wieczorek | Five simple yet essential steps to correctly estimate the rate of false differentially abundant proteins in mass spectrometry analyses[END_REF] and its supplementary material): > log2.int.ptm=log2(int.ptm.mod); 2. Plotting the distributions of observed values in all samples: The differences between the distributions of the observed values can be used to see whether samples have globally lower or higher quantified values than in the others (see Fig. 1C). > ggd = ggd + ylab("Height") > ggd = ggd + ggtitle("Ward's method with a \nPearson correlation-based distance") Such clustering can be useful for determining which biological conditions are similarly distributed, and which are less so (see Fig. 1E).

Is the number of replicates sufficient in each condition?

A sufficient number of samples to confidently trust the results of a statistical test can be determined by seeking the minimum number of samples reaching a satisfactory statistical power. The calculation of this statistical power depends on several factors. In case of a classical t-test comparing the average intensities between two conditions, it depends on the minimal average difference intensity between conditions that is expected to be detectable, the variance of the intensities, the type 1 error (threshold on the p-values), and the number of samples in each condition (sample size). That is why it is logical to conduct a power analysis from a first MS-based experiment to assess how many samples are sufficient to obtain robust statistical results in a subsequent analysis (see Note 15). This analysis can also be performed after the final data have been acquired to evaluate what is the confidence in the results. The R package ssize.fdr contains several functions that can be used to calculate minimum sample sizes [START_REF] Liu | Quick calculation for sample size while controlling false discovery rate with application to microarray analysis[END_REF]. The purpose of each of these functions is to estimate an average statistical power over all the tests carried out on all the peptides in function of the sample size. In case of more than two conditions, the functions ssize.F and ssize.Fvary can be used by defining a level of false discovery rate, a user-specified proportion of nondifferentially abundant peptides, a design matrix and a standard deviation (or the parameters of inverse gamma distribution followed by variances of peptides for ssize.Fvary). We detailed hereafter how to use the latter one.

1. Creating a design matrix from the vectors defining the design of experiment: Additional explanations on this matrix and the contrast matrix can be found in see Subheading 3.9,

Step 1 and see Subheading 3. #The proportion of non-differentially abundant peptides is #estimated using the distribution of the p-values related #to the contrasted test with default method ("pounds") > prop=estim.pi0(pval_ct) 5. Using the ssize.Fvary function: Additional explanations can be found in [START_REF] Liu | Quick calculation for sample size while controlling false discovery rate with application to microarray analysis[END_REF].

#FDR level required > fdr=0.01 #Minimum statistical power required > pwr=0.75 #Function computing the degrees of freedom related to the #design of experiment: here we have 4 groups (two conditions #at two time points) and two parameters estimated in the model > df=function(n){4*n-2;} #parameters of the inverse gamma distribution of the variances > alph=d0/2 > beta=d0*s02/2 #eps represents values for the model parameters related to #biological conditions and the time. Because we test the #nullity of the parameters, more this value is close to 0 #and the harder it will be on the statistical test to know #if the coefficients are indeed zero or not from data: #the average statistical power will logically be lower. > eps=0.1 > require(ssize.fdr) > ftv=ssize.Fvary(X=design,beta=c(1,eps,-eps),L=ct,dn=df,a=alph, b=beta,fdr=fdr,power=pwr,pi0= prop$pi0.est$pi0.Pounds,maxN=20) Fig. 2C-D represents the average statistical power in function of the sample size for different values of eps. The closer it is to 0, and lower the statistical power will be.

Normalisation step

As explained in the introduction, the enrichment step induces a lower reproducibility among the reported intensities. That is why a normalization step is of prime importance to correct the variations of intensities. In literature, several methods have been proposed [START_REF] Kauko | Label-free quantitative phosphoproteomics with novel pairwise abundance normalization reveals synergistic ras and cip2a signaling[END_REF] [START_REF] Saraei | Phosphonormalizer: an r package for normalization of ms-based label-free phosphoproteomics[END_REF]. Several normalization methods can be used from the wrapper.normalizeD function of the R package DAPAR [START_REF] Wieczorek | Dapar & prostar: software to perform statistical analyses in quantitative discovery proteomics[END_REF] after creating an object of class MSnSet. A classical approach consists to apply a median-centering function in the samples which are a replication of a same biological condition.

Dealing with missing values

In case of multiple biological conditions, modified peptides can have similar detection profiles: they are detected, or not detected, under the same biological conditions. These detection profiles may be of greatest interest. Indeed, non-detection is generally synonymous either with an absence of the modified peptide in the biological condition; or with its low abundance, as missing values mainly result of the limit of detection of the mass spectrometer. The detection profile of the modified peptide have to be compared with those of their belonging protein in order to highlight over-or under-abundance of the modification relatively to the one of the unmodified protein they belong to (see Note 16).

1. Creating a function returning detection profiles (1 if the modified peptide is detected and 0 if not): > detect.prof=function(mat,condition){ lev=levels(condition) detect=matrix(0,nrow(mat),length(lev)) for (k in 1:nrow(mat)){ for (i in 1:length(lev)){ for (j in which(condition ==lev

[i])){ if (!is.na(mat[k,j])){ detect[k,i]=1 } } } } colnames(detect)=levels(condition) return(detect) } 2.
Computing detection profiles for modified peptides and their proteins:

> detect.ptm=detect.prof(log2.int.ptm.norm,CondTime.ptm) > detect.prot=detect.prof(log2.int.prot.norm,CondTime.prot) Missing values can be replaced within a same condition using imputation algorithms. If no value have been measured along replications of a same experiment in a condition, then either a value can be inferred assuming that there is a relationship to the values measured in the other conditions, or no value can be inferred (see Note 17). Two main kinds of missing values arise in MS-based proteomics: either missing not at random (MNAR) values or missing completely at random (MCAR) values [START_REF] Lazar | Accounting for the multiple natures of missing values in label-free quantitative proteomics data sets to compare imputation strategies[END_REF]. The R package imp4p proposes a multiple imputation strategy to deal with these two kinds of missing values and functions to analyze missing value mechanisms [START_REF] Gianetto | A peptide-level multiple imputation strategy accounting for the different natures of missing values in proteomics data[END_REF].

3

Mapping the quantified values of the unmodified protein to the ones of modified peptides.

After the steps of data pre-processing, normalization and imputation, it is useful to create a dataset containing all the modified peptides retained on rows with their values by columns, as well as the values of the potential proteins to which they belong for subsequent statistical analysis. For each modified peptide, these intensities can be retrieved using the identification numbers of their protein groups in MaxQuant: they are in the column "id" of the "proteinGroups.txt" file.

The following R code creates a dataset containing all the information necessary for subsequent analysis:

> int.

Extracting modified peptides with specific detection profiles relatively to their unmodified protein.

Before analyzing the quantified intensities, specific detection profiles can be extracted to highlight modified peptides that are detected while their unmodified protein is not quantified for instance. For instance, this can be performed with the following steps for the "mock" condition of our case study dataset:

1 3.9 Extracting modified peptides evolving significantly differently from their unmodified protein.

A common strategy aims to focus on the evolution of the intensities of modified peptides at different time stamps and under different biological conditions, notably to highlight PTMmediated pathways involved in diseases or cellular disorders. For example, this can happen if one has samples related to patients belonging to different categories, and whose modified proteomes are measured at different time stamps. In this framework, the measured intensity of a modified peptide in a sample can be explained by three factors: the abundance of its membership protein, the patient's category and the time stamps. When no value is missing in the dataset, the following linear model can easily be estimated in each condition using R:

𝑦 𝑗 = 𝛼 + 𝛽1 𝑐𝑜𝑚 𝑝 ( 𝑗)= 𝑝𝑒 𝑝𝑡 + ∑ 𝑡 ∈ [𝑡 1 ,. . . ,𝑡 𝑛 ] 𝜃 𝑡 1 𝑡𝑖𝑚𝑒 ( 𝑗)=𝑡 + 𝛿 𝑡 1 𝑡𝑖𝑚𝑒 ( 𝑗)=𝑡 1 𝑐𝑜𝑚 𝑝 ( 𝑗)= 𝑝𝑒 𝑝𝑡 + 𝜖 𝑗 (1)
where 𝜖 𝑗 is a Gaussian white noise. In this equation, 𝑌 = (𝑦 𝑗 ) 𝑗 ∈ [1, 𝐽 ] is a vector containing the intensities of both the modified peptide and its unmodified protein, and the function 1 𝑐𝑜𝑚 𝑝 ( 𝑗)= 𝑝𝑒 𝑝𝑡 is equal to 1 if the 𝑗 𝑡 ℎ coordinate of the vector 𝑌 is an intensity related to the modified peptide and 0 if it is related to the protein. Similarly, the function 1 𝑡𝑖𝑚𝑒 ( 𝑗)=𝑡 is equal to 1 if the 𝑗 𝑡 ℎ coordinate of the vector y is an intensity related to the time stamp t and 0 otherwise. Here, the measured intensity 𝑦 𝑗 is explained in function of four parameters: the intercept 𝛼 which represents the average intensity level of the protein to which the peptide belongs at time 𝑡 0 ; the 𝛽 parameter which represents the gap between this average at 𝑡 0 and the one of the intensities of the modified peptide; the list of 𝜃 𝑡 which are the effects of each other time point regardless of whether the intensity is that of the peptide or the protein; and the list of 𝛿 𝑡 that correspond to the interaction effects between the peptide intensities and the time stamps.

Testing the nullity of a linear combination of the model coefficients allow to select peptides of interest. For this, a design matrix (defining the linear model) and a constrast matrix (defining the linear combination of model coefficients that have to be tested) has to be specified in R.

1. Defining the design matrix: The model (1) can also be represented as a product between a design matrix X, composed of ones and zeros, and a vector of parameters 𝐵 = (𝛼, 𝛽, 𝜃 To test if the modified peptide has a dynamic similar to that of the protein over time, we have to test that all the interaction parameters 𝛿 𝑡 are equal to 0, i.e. 𝛿 1 = 0 and 𝛿 1 -𝛿 2 = 0, 𝛿 2 -𝛿 3 = 0, etc. For instance, this contrast matrix will test 𝛿 1 = 0 and 𝛿 1 -𝛿 2 = 0:

𝐶 𝑇 𝐵 = 0 ⇐⇒ 0 0 0 • • • 0 1 0 • • • 0 0 0 0 0 • • • 0 1 -1 • • • 0 0 𝛼 𝛽 𝜃 1 . . . 𝜃 𝑇 𝛿 1 . . . 𝛿 𝑇 = 𝛿 1 𝛿 1 -𝛿 2 = 0
(2) In our case study dataset of two time stamps, only one interaction parameter is in the model defined by the design matrix. Using this model, the contrast matrix testing the nullity of this parameter is defined by only one column: > ct=cbind(c(0,0,0,1)) Alternatively, if it is desired to remove all the peptides which are constant over time, it is advisable to test the nullity of all the coefficients in relation to time, i.e. the nullity of the (𝜃 𝑡 ) and (𝛿 𝑡 ) parameters. In our case study, the contrast matrix to be used is: > ct=cbind(c(0,0,1,0),c(0,0,1,-1)) > p.value_LIMMA=fit$F.p.value 5. Selecting peptides with a significantly different dynamic than their parent unmodified protein with a chosen FDR threshold: For this, an adaptive FDR control procedure can be applied to get adjusted pvalues using the cp4p R package [START_REF] Gianetto | Calibration plot for proteomics: A graphical tool to visually check the assumptions underlying fdr control in quantitative experiments[END_REF]: > require(cp4p) > a=adjust.p(p.value_ANOVA,pi0.method = "pounds") > a.pvalue_ANOVA=a$adjp$adjusted.p > ptm.sel.mock=which(a.pvalue_ANOVA<0.01) Here, a 1% FDR threshold has been used and the proportion of true null hypotheses has been estimated using the method of Pounds and Cheng [START_REF] Pounds | Robust estimation of the false discovery rate[END_REF]. Likewise, modified peptides which have a significantly different dynamic from their unmodified protein in the SARS-CoV-2 condition can be recovered by using condition="SARS_COV_2" and listing them in a vector ptm.sel.sarscov2 from previous step (see Subheading 3.9,

Step 4).

6. Extracting peptides of interest: Finally, a dataset containing all the modified peptides evolving significantly differently than their unmodified proteins in at least one condition can be extracted from ptm.

Clustering of modified peptides using their dynamics relatively to their unmodified protein.

When prior information is not available to perform groupings of peptides, unsupervised clustering methods can be applied to highlight clusters of peptides behaving similarly for subsequent biological interpetation.

In label-free proteomics, many physico-chemical properties can lead to higher or lower intensity levels for different peptides with similar abundance in different conditions. Consequently, the cluster analysis of peptides has to be performed without considering the global intensity levels of the peptide, or the one of its protein, but rather their deviation from a reference intensity level. With this purpose, the intensities of the modified peptide and of its unmodified protein can be modeled by a linear model with interaction parameters using the three factors that are time stamps, biological conditions, and the studied component (protein or peptide):

𝑦 𝑗 = 𝛼 + 𝛽 0 1 𝑐𝑜𝑚 𝑝 ( 𝑗)= 𝑝𝑒 𝑝𝑡 + 𝑘 ∈ [1,𝐾 ] 𝛽 𝑘 1 𝑐𝑜𝑛𝑑 ( 𝑗)=𝑘 + 𝛾 𝑘 1 𝑐𝑜𝑛𝑑 ( 𝑗)=𝑘 1 𝑐𝑜𝑚 𝑝 ( 𝑗)= 𝑝𝑒 𝑝𝑡 + 𝑡 ∈ [𝑡 1 ,𝑡 𝑛 ] 𝜃 𝑡 1 𝑡𝑖𝑚𝑒 ( 𝑗)=𝑡 + 𝜗 𝑡 1 𝑡𝑖𝑚𝑒 ( 𝑗)=𝑡 1 𝑐𝑜𝑚 𝑝 ( 𝑗)= 𝑝𝑒 𝑝𝑡 + 𝑘 ∈ [1,𝐾 ] 𝑡 ∈ [𝑡 1 ,𝑡 𝑛 ] 𝛿 𝑡 𝑘 1 𝑡𝑖𝑚𝑒 ( 𝑗)=𝑡 1 𝑐𝑜𝑛𝑑 ( 𝑗)=𝑘 + 𝜑 𝑡 𝑘 1 𝑡𝑖𝑚𝑒 ( 𝑗)=𝑡 1 𝑐𝑜𝑛𝑑 ( 𝑗)=𝑘 1 𝑐𝑜𝑚 𝑝 ( 𝑗)= 𝑝𝑒 𝑝𝑡 +𝜖 𝑗 (3)
where 𝜖 𝑗 is a Gaussian white noise. Thus, the dynamics of the modified peptide compared to the one of its unmodified protein is characterized by all the parameters multiplied by 1 𝑐𝑜𝑚 𝑝 ( 𝑗)= 𝑝𝑒 𝑝𝑡 , with the exception of the 𝛽 0 parameter which represents the overall difference between the intensity level of the modified peptide and that of its unmodified protein. Therefore, the dynamic of the modified peptide relative to that of its protein can be characterized by the parameter vector

(𝛾 1 , • • • , 𝛾 𝐾 , 𝜗 𝑡 1 , • • • , 𝜗 𝑡 𝑛 , 𝜑 𝑡 1 1 , • • • , 𝜑 𝑡 𝑛 𝐾
) 𝑇 , while the dynamic of its protein is characterized by the parameter vector

(𝛽 1 , • • • , 𝛽 𝐾 , 𝜃 𝑡 1 , • • • , 𝜃 𝑡 𝑛 , 𝛿 𝑡 1 1 , • • • , 𝛿 𝑡 𝑛 𝐾 ) 𝑇 .
A clustering of these vectors can be performed to cluster peptides with similar dynamics.

Extracting model parameters using a classical approach:

> param=matrix(NA, nrow(int.ptm.prot.sel),6) > for (i in 1:nrow(int.ptm.prot.sel)){ Int=int.ptm.prot.sel[i,] mod <-lm(Int ~Cond.ptmprot*Time.ptmprot*Comp.ptmprot) param[i,]=mod$coefficients[unique(c(grep("Cond", names(mod$coefficients)), grep("Time",names(mod$coefficients))))] } > colnames(param)=names(mod$coefficients[unique( c(grep("Cond",names(mod$coefficients)), grep("Time",names(mod$coefficients))))]) #Extract coefficients only related to the dynamics of the #modified peptides relatively to their protein > param_pept=param[,grep(":Comp.ptmprot",colnames(param))] 2. Extracting model parameters using limma R package:

> require(limma) > design=model.matrix(~Cond.ptmprot*Time.ptmprot*Comp.ptmprot) > fit=eBayes(lmFit(int.ptm.prot.sel, design),robust=TRUE) > param_LIMMA=fit$coefficients[,unique(c(grep("Cond", names(mod$coefficients)), grep("Time", names(mod$coefficients))))] 3. Choosing the clustering algorithm: Many unsupervised clustering algorithms proposed in the literature could be applied on our dataset. We propose in the next to use the classical "hard" version of the k-means algorithm (see Note 19). This algorithm depends on two main parameters that have to be shortlisted by the user: the distance measure and the number of clusters.

4. Distance measure and dissimilarity matrix: Generally, Euclidean distances or Pearson correlation-based distances between observed values are used (see Note 20). The Euclidean distance aims to measure the absolute deviations between the vector coordinates, while the Pearson correlation-based distance characterizes the linear evolution between the vector coordinates. In our case, we search to cluster parameter vectors and no dynamic is expected between these parameters. Thus, a classical Euclidean distance seems to be an appropriate choice. Once the distance chosen, a dissimilarity matrix can be computed. To this end, the multivariance R package [START_REF] Böttcher | Copula versions of distance multivariance and dhsic via the distributional transform-a general approach to construct invariant dependence measures[END_REF] provides a function fastdist allowing fast computation of Euclidean distance matrix, which is of definite interest with respect to the data set size (several tens of thousands of peptides):

> require(multivariance) > diss.mat=fastdist(param) 5. Determining the optimal number of clusters: Many indices can be used to determine an optimal number of clusters [START_REF] Charrad | Nbclust: An r package for determining the relevant number of clusters in a data set[END_REF]. A classical method consists in using a Gap statistic [START_REF] Tibshirani | Estimating the number of clusters in a data set via the gap statistic[END_REF]. For this, the fviz_nbclust function of the R package factoextra can be used: b. Displaying the two clusters with the highest average norms as in Fig. 3E: > require(reshape2) #for melt > require(ggpubr) #for facet_grid > nb_to_plot=2 > cl_to_plot=as.numeric(names(mc)[1:nb_to_plot]) > dataplot=dataplot[dataplot$cluster%in%cl_to_plot,] > datap=melt(dataplot,id.var=c("cluster","peptID"))

>
> gg= ggplot(datap, aes(x = variable, y = value) ) > gg=gg+geom_boxplot(aes(fill = variable), alpha = 0.5) > gg=gg+facet_grid(cluster~variable,scales="free_x") > gg=gg+scale_x_discrete(labels = "") > gg=gg+theme(legend.position = "none") > gg

Subsequent functional analysis

1. Motif enrichment analysis: Motif analysis can be used to know the main sequence of amino acids surrounding modification sites of interest. For instance, this kind of analysis can be useful to find out which kinases are involved in phosphorylation. In R, motif analysis can be done using the rmotifx R package [START_REF] Wagih | Uncovering phosphorylation-based specificities through functional interaction networks[END_REF] and the ggseqlogo R package can be used to visualize the results [START_REF] Wagih | ggseqlogo: a versatile r package for drawing sequence logos[END_REF] (see Note 21).

2. Enrichment analysis of Gene Ontology terms and Pathways: An analysis of the enrichment of Gene Ontology terms or biological pathways referenced in databases such as KEGG or Reactome is often carried out to highlight biological processes or specific signaling pathways enriched in the lists of modified peptides of interest coming from statistical analysis. Two important points have to be kept in mind when performing these analyses in our framework: i) This enrichment is generally carried out in relation to a "background" using hypergeometric tests. In our context, it is very important to choose the set of modified proteins identified by mass spectrometry as background. Indeed, the enrichment protocols can display a bias in favor of enrichment for peptides with a single modification versus peptides with several modifications, or towards enrichment for hydrophobic versus hydrophilic modified peptides. Such biases will not be taken into account if the whole proteome or genome of the organism studied is used as background.

ii) Most enrichment tools are based on gene-centric databases. However, two modified peptides of the same protein can sometimes evolve differently between two conditions, for example if the abundance of one proteoform of the protein decreases while that of another increases. In such a case, using the annotations of the reference protein without taking the modification into account may lead to falsely highlighting certain biological processes that are associated with a particular proteoform and not with another. Fortunately, recent resources have been developed to take this aspect into account [START_REF] Krug | A curated resource for phosphosite-specific signature analysis[END_REF].

the Cytoscape App Store are particularly interesting for visualizing PTMomics data: stringApp [START_REF] Doncheva | Cytoscape stringapp: network analysis and visualization of proteomics data[END_REF] and Omics Visualizer [START_REF] Legeay | Visualize omics data on networks with omics visualizer, a cytoscape app[END_REF]. stringApp is an application allowing to query the protein-protein interaction database of STRING [46] from Cytoscape, while Omics Visualizer allows to optimize the visualization of the obtained network. The latter is particularly suitable for viewing information on each node of the network, for instance the number of modified peptides regulated for each protein. Specific tutorials are freely available online (https://jensenlab.org/training/).

Notes

1. A number of recently proposed packages allow one to consider pipelines of analyzes performed entirely from R: for example, the rawR package can read RAW files from mass spectrometers [START_REF] Kockmann | rawr-direct access to raw mass spectrometry data in r[END_REF], while the Bioconductor packages rTandem [START_REF] Fournier | rtandem, an r/bioconductor package for ms/ms protein identification[END_REF] or MSGFplus [START_REF] Kim | Ms-gf+ makes progress towards a universal database search tool for proteomics[END_REF] [50] can be used to perform identification of mass spectra from these RAW files. Additionally, the RforProteomics package offer useful functionnalities to manipulate and visualize mass spectra [START_REF] Gatto | Visualization of proteomics data using r and bioconductor[END_REF].

2. The first line of the input file has to contain the column names of the data table. Decimal separator for quantitative values has to be the dot ".". Other R functions, such as read.xlsx() , from the R package openxlsx ; or read_excel() from the R package readxl can also be used if the file is exported from Excel.

3. These peptides are used to calculate a false discovery rate linked to the identification of peptides [START_REF] Couté | Beyond target-decoy competition: Stable validation of peptide and protein identifications in mass spectrometry-based discovery proteomics[END_REF].

4. The higher the probability of localization, the more certain it is that the measured spectrum is associated with the peptide containing the modification.

5. The "proteinGroups.txt" contains intensities for the proteins using either their unmodified peptides or their modified peptides in different columns. In the case study dataset, samples enriched in phosphorylation contain "Phospho", while the ones enriched in ubiquitination contain "Ubi". However, these names can be different with other datasets depending on how they have been defined in MaxQuant or in the used software. Thus, the reader has to adapt these names in such cases. The colnames function can be used to check these names.

6. Similarly to previous note, the names of enriched samples can be different with other datasets depending on how they have been defined in MaxQuant or in the used software: the reader has to adapt these names in such cases.

7. These peptides are generally present because they have been identified. However, no quantified values are available because of the quantification algorithm used in MaxQuant.

8. The same approach can be used with 2, 4 or 5 samples using the draw.pairwise.venn , draw.quad.venn and draw.quintuple.venn functions of the VennDiagram package. 9. It is generally useful to convert a plot into a ggplot object, for example to easily include it in a pdf or powerpoint document created automatically with R Markdown or the officer R package. 10. This can be performed using the qplot function of R package ggplot2 with the following R code: > require(reshape2) > gradientRate=1.2; > text=element_text(colour="black",size=8,face="bold"); > plot.j=qplot(x=Var1,y=Var2,data=melt(mat.jaccard),fill=value, geom="tile") > plot.j=plot.j+theme(axis.text=element_text(size=10), axis.title=element_text(size=10,face="bold"), axis.text.x=element_text(angle=90,vjust=1,hjust=1),legend.text=text) > plot.j=plot.j+theme(legend.title=element_blank()) > plot.j=plot.j+labs(x="",y="") > plot.j=plot.j+scale_fill_gradientn(colours=colorRampPalette( c("white","lightblue","darkblue"))(101),values=c(pexp(seq(0,1,0.01), rate=gradientRate),1),limits=c(0,1)) 11. For this, the same approach has to be applied from vectors equal to 1 is the modified peptide is not identified and 0 if it is identified. Of note, the Jaccard index does not give the same value when evaluating mutual presence or mutual absence, such that the result of this second approach will be different from the first. An alternative to the Jaccard index consists to use the Rand index, this one will lead to a single value measuring at the same time mutual presence and absence. 12. In this way, the differences observed between the various distributions will highlight whether if the intensities measured in one sample are lower or higher than in others, without confounding effects related to different sets of peptides used to plot the distributions. 13. This can be performed similary as in Note 10, except that data=melt(mat.cor) in the qplot function. Also, we advise to choose gradientRate= 3.5 for correlation matrix based on quantified values to see more clearly clusters of samples in the matrix. 14. Because the Principal Component Analysis (PCA) is based on the eigendecomposition of the covariance matrix, it will summarize the same information that appears in the correlation matrix. The prcomp function can be used to perform PCA in R. [START_REF] Chambers | Software for data analysis: programming with R[END_REF]. In general, it is necessary to carry out a first MS-based experiment with a minimum number of replicates per biological condition (at least 3 to be able to estimate the variance of intensities for each modified peptide in each condition). This kind of test experiment is generally intended to check that the enrichment step has worked correctly but can also be used to assess how many samples are sufficient to obtain robust statistical results in a subsequent analysis. [START_REF] Millot | Comprendre et réaliser les tests statistiques à l'aide de r[END_REF]. In each condition, there are therefore two possibilities: either the modified peptide is detected in one of the samples of the condition, or it is not. In the case of 𝑁 conditions, 2 𝑁 -1 detection profiles are possible. These detection profiles should be compared with the ones associated with the corresponding proteins, which further increases the number of possible cases. Hopefully, it can generally be expected that many detection profiles will not be encountered in the experiment. It should be noted that some cases lead to not being able to conclude on the over or under abundance of the modified peptide compared to the unmodified protein. For example, if they are detected in the same condition and not in another, nothing can be concluded by comparing these conditions. 17. For example, with conditions representing different time stamps, a value can be deduced by applying methods developped for time series data [START_REF] Pratama | A review of missing values handling methods on time-series data[END_REF]. In the case where no value can be inferred in the condition from observed ones, then an analysis of the detection profiles can be conducted to highlight interesting modified peptides. 18. Here, the "*" operator is used to specify that interaction parameters has to be included in the model. If one does not want interactor parameters, one uses only the "+" operator. [START_REF] Wickham | ggplot2: elegant graphics for data analysis[END_REF]. Either hierarchical clustering methods or partitional clustering methods (such as the k-means algorithm) can be used [START_REF] Gan | Data clustering: theory, algorithms, and applications[END_REF]. Moreover, either hard clustering versions of these methods, assigning a unique cluster for each peptide, or fuzzy clustering versions, assigning probabilities of belonging to each cluster for each peptide could be used [START_REF] Schwämmle | Vsclust: feature-based variance-sensitive clustering of omics data[END_REF]. The main advantage of hierarchical methods is that they does not need to specify a number of clusters. Their results are generally represented in the form of a dendrogram. However, in our context, it is preferable to cluster the peptides into a finite number of clusters in order to facilitate the subsequent biological interpretation of the results. Partitional clustering methods seem thus more appropriate. The most famous one is the k-means algorithm. The fuzzy version of this algorithm has serious limitations in high dimensional spaces, as it is generally the case in peptidomics datasets [START_REF] Winkler | Fuzzy c-means in high dimensional spaces[END_REF]. 20. However, alternative distances can also be investigated as, in case of time series, Dynamic Time Warping distances [START_REF] Giorgino | Computing and visualizing dynamic time warping alignments in r: the dtw package[END_REF] or Short time series distances [START_REF] Mori | Distance measures for time series in r: the tsdist package[END_REF]. 21. Pre-alignment of the modified peptides around the detected sites of modification has to be performed before to use these packages. Here, param_pept contains coefficients only related to the dynamics of the modified peptides relatively to their protein. The cluster 1 is characterized by peptides with positive coefficients for Time 24h and SARS-Cov-2 conditions, as well as negative coefficients for their interaction.

6 3 .

 3 Using the draw.triple.venn() function of R package VennDiagram to plot a Venn diagram (see Note 8):

  5): > require(UpSetR) > data.upset=as.data.frame(tr.int.ptm) > names(data.upset)=colnames(data.upset) > data.upset[is.na(data.upset)]=0 > upset.plot = upset(data.upset, order.by = "freq", number.angles=315, main.bar.color=4, nsets=ncol(data.upset), sets.x.label="Nb", text.scale = c(1.5,1.5,1.5,1.5,1.5,1)); #To visualize > upset.plot 5. {Optional step} Converting the plots into ggplot objects: It can be useful to convert the obtained UpSet plot or the Venn diagram into ggplot objects (see Note 9). For this, the grid and ggplotify R packages can be used: > require(ggplot2) > require(ggplotify) > require(grid)

Figure 1 :

 1 Figure 1: Graphs to check the reproducibility of experiments. A: Venn diagram checking the reproducibility of modified peptides identified in SARS-CoV-2 samples at 6h. B: UpSet graph checking the reproducibility of modified peptides identified in all samples (3871 are found in common). C: Estimated distributions of log2(intensities) of modified peptides found in all samples. This makes it possible to detect shifts in all the intensity values of one or more samples. D: Pearson correlation matrix between all samples using shades of blue. Darker blues (higher correlation values) should appear between replicates of the same condition. E: Hierarchical clustering of samples from Pearson correlation values. Samples of a same condition should cluster together, and not with the one of another condition if the conditions proteomes are not close. If they are, clustering may have a hard time distinguishing them, but they will appear close (with darker blue) in the correlation matrix plotted in D.

Figure 2 :

 2 Figure 2: A: Histogram of pvalues related to the used test. B: Calibration plot to choose a method to estimate the proportion of true null hypotheses. C: Average power in function of the sample size with eps=0.1 and the proportion of true null hypotheses estimated by the method of Pounds and Cheng [33]. D: Average power in function of the sample size with eps=0.5 and the proportion of true null hypotheses estimated by the method of Pounds and Cheng [33].

Figure 3 :

 3 Figure 3: A: Gap statistic in function of the number of clusters for param vectors, see Subheading 3.10, Step 1. B: Gap statistic in function of the number of clusters for param_pept vectors, see Subheading 3.10, Step 1. C: Visualization using PCA of 20 clusters of param vectors found with kmeans. D: Visualization using PCA of 11 clusters of param_pept vectors found with kmeans. E: Profiles of parameters for the clusters with the highest average norms for param_pept vectors.Here, param_pept contains coefficients only related to the dynamics of the modified peptides relatively to their protein. The cluster 1 is characterized by peptides with positive coefficients for Time 24h and SARS-Cov-2 conditions, as well as negative coefficients for their interaction.

Note 5):

  For the proteins, only the columns relative to the quantification of the unmodified proteins can be kept using (see

	6. Distinguishing between the different versions of a modified peptide: In MaxQuant
	outputs, the quantified intensities derived from one two or three and higher PTMs are in > Cond.ptm=factor(c(rep("mock",6),rep("SARS_COV_2",6)),levels=
	columns whose name ends with ___1, ___2, ___3. Thus three lines must be created for c("mock","SARS_COV_2"))
	each modified peptide of the file. > Time.ptm=factor(rep(c(rep("24h",3),rep("6h",3)),2),levels=
	> int.ptm.mod=rbind(int.ptm[,grep("___1",colnames(int.ptm))], c("6h","24h"))
	int.ptm[,grep("___2",colnames(int.ptm))], > CondTime.ptm=as.factor(paste(Cond.ptm,Time.ptm,sep="."))
	int.ptm[,grep("___3",colnames(int.ptm))])
	> colnames(int.ptm.mod)=unlist(strsplit(colnames(int.ptm)[
	grep("___1",colnames(int.ptm))],split ="___1"))
	Each modified peptide is next characterized by its identification number, the proteins
	to which it may belong (quantified in the "proteinGroups.txt" file) and the number of
	PTMs used to quantify ("multiplicity" column):
	> id.ptm=rep(data.ptm$id,3)
	> id.pg=rep(data.ptm$Protein.group.IDs,3)
	> multiplicity=c(rep("___1",nrow(int.ptm)),rep("___2",
	nrow(int.ptm)),rep("___3",nrow(int.ptm)))
	7. Removing modified peptide with no quantification values: The modified peptides asso-
	ciated with no quantification value in any sample are useless for performing subsequent
	statistical analysis (see Note 7):
	#Compute the number of observed values on each row
	> sum_notna=apply(int.ptm.mod,1,function(x)sum(!is.na(x)))
	#Extract peptides with at least one observed value
	#on each row
	> int.ptm.mod=int.ptm.mod[which(sum_notna>0),]
	> id.ptm=id.ptm[which(sum_notna>0)]
	> id.pg=id.pg[which(sum_notna>0)]
	> multiplicity=multiplicity[which(sum_notna>0)]
	8. Renaming the columns of the dataset: It can generally be useful to rename the columns
	with shorter names for further analyses, in particular to make them appear in graphics:
	> short.name=unlist(strsplit(colnames(int.ptm.mod),
	split="Intensity.Phospho_"))
	> short.name=short.name[short.name!=""]
	> colnames(int.ptm.mod)=short.name
	> data.ptm.mod=data.frame(id.ptm,id.pg,multiplicity,
	int.ptm.mod)
	> int.prot=int.prot[,!colnames(int.prot)%in%grep("Phospho", 9. Defining vectors containing the design of experiments: The experimental design is
	colnames(int.prot), value=TRUE)] specific to each experiment. It is extremely important, especially for the subsequent
	> int.prot=int.prot[,!colnames(int.prot)%in%grep("Ubi", statistical analysis. The vectors containing the experimental design must correspond
	colnames(int.prot), value=TRUE)] to the columns containing the intensity values. It is advisable to do one for modified
	Similarly, the intensities of the modified peptides can be extracted using (see Note 6): peptides, another for proteins, and one containing both. In our case study, we have two
	> int.ptm=as.matrix(data.ptm[,grep("Intensity.Phospho", factors, the time (6h or 24h) and the biological conditions (mock or SARS-Cov-2 infected
	names(data.ptm), value=TRUE)]) samples), 3 replicated samples have been performed in each condition for studying the
	5. Replacing 0 by NA: This step is useful to easily detect missing values in subsequent modifications and 4 replicates have been used for the unmodified proteins. Experimental
	analyses: design corresponding to the columns of int.ptm.mod is defined by:
	> int.prot[int.prot==0]=NA #colnames(int.ptm.mod) can be used to check the column names
	> int.ptm[int.ptm==0]=NA #of int.ptm.mod

Note 4): > data.ptm=data.ptm[which(data.ptm$Localization.prob>0.75),] 4. Separating intensities from metadata: The metadata generally contains multiple information, notably related to the identification of the modified peptides and the proteins they belong to. To separate the quantified intensities from the metadata, columns containing quantified intensities can be searched using the grep function, and the data can be converted into a matrix object in R with the as.matrix function: > int.prot=as.matrix(data.prot[,grep("Intensity.", names(data.prot), value=TRUE)])

  The proportion of missing values is almost 30% in each sample. The proportion of MCAR values is estimated between 0 and 30% depending on the considered sample.

	. Analyzing missing value mechanisms:
	> require(imp4p)
	> log2.int.ptm.impMCAR=impute.mle(log2.int.ptm.norm,
	conditions=CondTime.ptm)
	> resmod=estim.mix(tab=log2.int.ptm.norm,
	tab.imp=log2.int.ptm.impMCAR,conditions= CondTime.ptm)
	> resmod$pi.na
	> resmod$pi.mcar

4. Imputing missing values of peptides using a multiple imputation strategy combining MCAR and MNAR values with the impute.mi() function: > log2.int.ptm.imp=impute.mi(tab = log2.int.ptm.norm, conditions = CondTime.ptm) > colnames(log2.int.ptm.imp)=colnames(log2.int.ptm) 5. Imputing missing values of proteins using a MCAR-devoted imputation algorithm: > log2.int.prot.impMCAR=impute.mle(log2.int.prot.norm, conditions=CondTime.prot) > colnames(log2.int.prot.impMCAR)=colnames(log2.int.prot.norm)

  .ptm.prot=data.frame(id.ptm.mod,id.pg.mod,id.prot.mod, multiplicity.mod,prof.ptm,prof.prot,int.ptm.prot)

	if (length(data.prot$id%in%id_prot[j])!=0){
	prof.prot=rbind(prof.prot,
	detect.prot[which(data.prot$id%in%id_prot[j]),])
	int.ptm.prot=rbind(int.ptm.prot,c(log2.int.ptm.imp[i,],
	log2.int.prot.impMCAR[which(data.prot$id%in%id_prot[j]),]))
	}else{
	prof.prot=rbind(prof.prot,rep(NA,ncol(detect.prot)))
	int.ptm.prot=rbind(int.ptm.prot,c(log2.int.ptm.imp[i,],
	rep(NA,ncol(log2.int.prot.impMCAR))))
	}
	}}
	> data
	ptm.prot=NULL;id.ptm.mod=NULL;
	> id.pg.mod=NULL;multiplicity.mod=NULL;
	> prof.ptm=NULL;prof.prot=NULL;id.prot.mod=NULL;
	> for (i in 1:nrow(data.ptm.mod)){
	id_prot=unlist(strsplit(data.ptm.mod$id.pg[i],split =";"))
	for (j in 1:length(id_prot)){
	id.ptm.mod=c(id.ptm.mod,data.ptm.mod$id.ptm[i])
	id.prot.mod=c(id.prot.mod,id_prot[j])
	id.pg.mod=c(id.pg.mod,data.ptm.mod$id.pg[i])
	prof.ptm=rbind(prof.ptm,detect.ptm[i,])
	multiplicity.mod=c(multiplicity.mod,
	data.ptm.mod$multiplicity[i])

  1 , . . . , 𝜃 𝑇 , 𝛿 1 , . . . , 𝛿 𝑇 ) 𝑇 , such as 𝑌 = 𝑋 𝐵 + 𝜖, where 𝜖 = (𝜖 𝑗 ) 𝑗 ∈[1,𝐽 ] . In R, the design matrix can be defined by using the model.matrix function (see Note 18): > design = model.matrix(~Comp.ptmprot*Time.ptmprot) 2. Defining the contrast matrix: This matrix (referred to as 𝐶) is used to test if a linear combination of coefficients composing the estimated model is null. In this framework, the null hypothesis of the test is defined by 𝐻 0 : 𝐶 𝑇 𝐵 = 0 and the alternative hypothesis by 𝐻 1 : 𝐶 𝑇 𝐵 ≠ 0.

  [START_REF] Pieroni | Enrichments of post-translational modifications in proteomic studies[END_REF]. For the following steps, we consider a dataset without missing values in order to be able to estimate linear models and perform statistical tests:

	nbna==0),]
	int.ptm.prot=int.ptm.prot[which(nbna==0),]
	4. a. Applying a classical Fisher's F test: Once the design matrix and the contrast matrix
	have been defined, a classical Fisher's F test can be applied in a condition (here, "mock")
	for all the modified peptides of our dataset using the linearHypothesis function of
	the car R package [26]:
	> require(car)
	> condition="mock"
	#condition="SARS_COV_2"
	> p.value_ANOVA=rep(NA, nrow(int.ptm.prot))
	> for (i in 1:nrow(int.ptm.prot)){
	Int=int.ptm.prot[i,Cond.ptmprot==condition]
	Comp=Comp.ptmprot[Cond.ptmprot==condition]
	Time=Time.ptmprot[Cond.ptmprot==condition]
	mod <-lm(Int ~Comp * Time)
	p.value_ANOVA[i]=linearHypothesis(mod,t(ct))$'Pr(>F)'[2]
	}
	b. Applying a Fisher's F test using limma: The limma R package can alternatively be
	used to test the nullity of the linear combination of coefficients with Fisher's F test. It
	uses a regularisation of the variances based on the assumption that they follow an inverse
	Gamma distribution [22]:
	> require(limma)
	> design=model.matrix(~Comp.ptmprot[Cond.ptmprot==condition]
	*Time.ptmprot[Cond.ptmprot== condition])

nbna=apply(int.ptm.prot,1,function(x)sum(is.na(x))) data.ptm.prot=data.ptm.prot[which(> fit=eBayes(contrasts.fit(lmFit(int.ptm.prot[,Cond.ptmprot== condition], design),ct),robust=TRUE)

  There is two main ways to visualize cluster results. One way is to use dimension reduction methods as Principal Component Analysis or Multidimensional Scaling Methods. The graph=TRUE option of the eclust function or the fviz_cluster function can be used with this purpose. A second way is to plot the profiles of the model parameters used to cluster. The clusters associated to strong deviations from 0 of these parameters are synonym of profiles strongly discriminant, and are thus of main interest. These clusters can be highlighted by computing the average norms of parameter vectors.

	require(factoextra)
	> gap.stat=fviz_nbclust(param,diss = diss.mat, kmeans,
	nstart = 25, method = "gap_stat", nboot = 10, k.max = 30)
	> gap.stat=gap.stat+labs(subtitle = "Gap statistic method")
	> gap.stat
	Fig.3A-B display the Gap statistics in function of the number of clusters for param and
	param_pept vectors see Subheading 3.10, Step 1.
	6. Performing the clustering: The eclust function of factoextra R package can be
	used for enhancing the workflow of clustering analyses and ggplot2-based elegant data
	visualization.
	#Using a k-means clustering with 20 clusters.
	> res.km=eclust(param, hc_metric = "euclidean", "kmeans",
	k=20, nstart = 10, graph=TRUE)
	Alternatively, the k-medioids algorithm can be used using the pam R function from
	cluster R package, while a visualization of this clustering using Principal Component
	Analysis is obtained with the fviz_cluster function:
	> require(cluster)
	> res.pam=pam(x=diss.mat,k=20,diss=TRUE)
	> res.pam$data=param
	> fviz_cluster(res.pam,main = "K MEDIOIDS Clustering")
	7. Visualization of results: a. Estimating the average norms of parameter vectors by cluster:
	> dataplot=data.frame(param_pept,cluster=res.pam$clustering,
	peptID=1:nrow(param))
	#compute the norms of parameter vectors for param_pept
	> normv=apply(param_pept,1,function(x)sqrt(sum(x^2)))
	#average norm for each cluster
	> mc=NULL
	> for (i in 1:length(levels(as.factor(dataplot$cluster)))){
	mc=c(mc,mean(normv[dataplot$cluster==i]))
	}
	> names(mc)=levels(as.factor(dataplot$cluster))
	#rank clusters by their average norms
	> mc=mc[order(-mc)]

R can be freely installed on most operating systems (Windows, Linux, Mac OS) from https://cran.r-project.org/. The basic R installation includes the language itself as well as packages, most of which are geared towards statistical analysis. One of the strengths of R is the large number of packages that can be downloaded, mainly from the CRAN repository (https://cran.r-project.org/) and, as far as bioinformatics is concerned, from Bioconductor (http://bioconductor.org/). In addition to R, we recommend that the reader install RStudio

Visualizing networks of PTMomics data: This step can be insightful to find modified proteins with known interactions. To do so, the reader can start from the statistical analyzes carried out from R, export them and then use the Cytoscape software[START_REF] Shannon | Cytoscape: a software environment for integrated models of biomolecular interaction networks[END_REF] to visualize interaction networks. Two Cytoscape applications freely downloadable from

Acknowledgements:

The author wants to acknowledge Mariette Matondo, Thibaut Douché, Thibault Chaze and Magalie Duchateau of the Proteomics platform of the Institut Pasteur for fruitful discussions about proteomics and phosphoproteomics.