Proteome-wide PPToP

For PPToP MS raw files were processed using MaxQuant (version 1.6.4.0)⁴¹ using a reference human proteome (uniprot Proteome ID: UP000005640, downloaded 9.6.2020). Data were processed separately for total and phospho-enriched samples, but for each all time points, fractions, and replicates were run together. Default search parameters were used, except as follows: multiplicity: 2; Heavy channel: Arg10, Lys8; variable modifications: Acetyl (Protein N-term), Oxidation (M), and only for the phosphofraction: Phospho (STY); fixed modifications: Carbamidomethyl (C); maximum number of modifications per peptide: 5; maximum missed cleavage sites: 2 (3 for phospho); LFQ: none; re-quantify: unchecked; match between run: checked.

Identified peptides from the MaxQuant evidence file were filtered to remove hits from the reverse database and potential contaminants. All subsequent analysis was done on the modified peptide level (henceforth referred to as “peptide”) including information on Heavy amino acid incorporation, N-terminal acetylation (N-Ac), and phosphorylation, but excluding methionine oxidation. Peptides quantified in only one of the SILAC channels (constituting 45.3% of all identified peptides) were removed. In case a peptide was quantified multiple times, a single entry was chosen by choosing the species with (i) the lower posterior error probability (PEP), and (ii) the highest intensity. Peptides were further filtered for presence in at least 2 replicates and 2 time points.

Cell cycle times for each replicate were estimated from the protein median values in the unmodified fraction. Assuming exponential decay of most proteins, we have the linear relationship¹:

where newold is the SILAC ratio of new and old protein, kdeg is the protein-specific degradation constant, and tcc the cell cycle time. We estimated tcc from the 1% longest-lived proteins where, assuming no active protein degradation (kdeg=0), the slope is defined by ln2tcc, from which we got tcc estimates of 28.0 h, 26.5 h, 27.0 h, and 22.2 h for replicates 1 through 4, respectively. Subsequently, new/old SILAC ratios were transformed into fraction of old remaining (see Supplementary Note 1 for reasons of doing so), and corrected for cell cycle as follows:

Next, peptide entries were filtered for reproducibility between replicates by calculating the distance of the corrected old remaining (φ) value for each peptide for each replicate to the median value of that peptide at that time point over all replicates, and excluding entries with values deviating from the median by more than two standard deviations of the entire distance distribution of all peptides at all time points in that fraction (unmodified or phospho). This removed 2.3% of all measurements.

The clearance of each peptide was compared to the median of all other quantified peptides of that protein from the unmodified fraction by fitting a spline with 3 degrees of freedom to the trace of φ vs time excluding the last time point (28 h). Peptides with data in at least 4 time points, a total of at least 6 data points, and for which the median of the other peptides in that protein included at least 2 unique peptides in the unmodified fraction were included. An F-statistic was calculated for fitting the spline to either both peptide and protein median together (H0 model) or to each separately (H1 model). Due to heteroscedasticity of the data, the resulting F-statistic was calibrated as delineated in ref. ⁴² by estimating the “effective degrees of freedom” (d₁, d₂) from a null dataset, where “peptide” or “median” labels were randomized (thus reducing any differences between the two classes to those occurring by chance). Since our dynamic PTM-SILAC dataset had differing amounts of data points per case, d₁, d₂ were estimated for six separate bins of the data depending on the number of data points available for the comparison. The thusly calibrated F-statistic distribution across both fractions was used to calculate p values for each peptide using the pf() function from the stats package in R and corrected for multiple testing using Benjamini–Hochberg correction. Cases with adjusted p value <= 0.001 were considered as hits.

To compare behavior of peptides in their phosphorylated and unmodified form (e.g., Fig. 3b, ​b,c),c), phosphosites were collated on the site level, i.e., on the specific modified amino acid. Consequently, the specific peptide boundaries were disregarded, when matching unmodified and modified peptides as phosphorylation can cause miscleavages and thus change the tryptic peptide boundaries.

Prediction of protein disorder was taken from the D2P2 database⁴³ (https://d2p2.pro/) using a consensus threshold of 75% across the individual predictor algorithms when determining the disorder status per amino acid. For the enrichment analysis, a peptide was considered intrinsically disordered if it contained at least 40% disordered amino acids.

Peptides were considered to be proximal to a degron sequence, if a stretch of +/−15 amino acids around the center of each peptide (31 amino acids altogether, unless the peptide was close to the N or C terminus) showed any overlap with a degron. Predicted degrons were taken from ref. ²⁵ verified degrons from ref. ²⁶.

Statistical tests were done in R using the following functions: Fisher’s exact test (Figs. 3b, ​b,3d,3d, ​d,4c,4c, ​c,6f):6f): fisher.test(); t test (Figs. 3c, ​c,3i,3i, ​i,6d,6d, ​d,6h,6h, Supplementary Fig. ^2g) and Wilcoxon signed-rank test (Figs. 2h, ​,6g,6g, Supplementary Fig. ^2f): t.test() and stat_compare_means(); linear regression: lm(); Spearman correlation (Fig. 3e, Supplementary Fig. ^6b): cor.test(); Pearson correlation (Fig. 7f): stat_cor(). Sample sizes are shown in each figure. Other covariates tested: on protein level: subcellular location, length; on peptide level: length, presence of signal or propeptides, intensity, non-tryptic cleavage, functional score¹⁴, in interface with other proteins, motifs, secondary structure predictions, predicted phosphorylating kinase.