dataesr/scientific-paragraphs-categorization huggingface dataset:
Comparison True log2 FC without adjustment True log2 FC after adjustment mix2 vs mix1 -1 -1 mix3 vs mix1 0 1 mix4 vs mix1 -1 0 mix3 vs mix2 1 2 mix4 vs mix2 0 1 mix4 vs mix3 -1 -1  As input to our workflow, we used the MSstatsPTM_Summarized. rda object provided on MassIVE. In the FeatureLevel data part of the object, the spiked-in peptides were not annotated and were irretrievable because the heavy peptides can also be present as their light counterparts. However, they were annotated in the ProteinLevel part. Hence, we could not use the low-level data, and had to start from the data that had already been pre-processed and summarised to PTM (for the PTM dataset) and Protein level (for the global profiling dataset) by MSstatsPTM, thus omitting step 4 and 5 from the workflow. We therefore could not assess our entire workflow based on these data, and moreover do not know which preprocessing steps were conducted.
We employed various methods to analyse this dataset. Our primary approach was the msqrob2PTM workflow as described in the workflow section, as well as the normal MSstatsPTM workflow. We also assessed the differential abundance of the PTMs with the standard msqrob2 workflow: DPA-nonNorm, which does no normalisation and hence skips step 3 of the standard workflow entirely and DPA, which only applies median centring in step 3. Because we know the ground truth of this dataset, we can again use the same metrics to assess the performance. Here, we also make ROC (tpr-fpr) curves. Furthermore, the log fold changes estimated by msqrob2PTM and MSstatsPTM were used to generate boxplots showing the observed and expected FCs for each mixture. For MSstatsPTM, the log fold changes were derived from the MSstatsPTM_Model.
Ubiquitination dataset Details of the experimental set-up can be found in reference  (24). The dataset itself is available on MassIVE(25) as MSV000078977. The dataset consists of four conditions: carbonyl cyanide 3-chlorophenylhydrazone (CCCP) treatment, USP30 overexpression (USP30-OE), a combination of both (Combo), and a control group. Per condition, two biological replicates with two technical replicates each were generated. All pairwise comparisons were tested using msqrob2PTM. This dataset has also been used in the MSstatsPTM paper, hence, we can compare our results to theirs for a biological case with unknown ground truth. As input to the msqrob2PTM workflow, we used the usp30_input_data. rda object found in the MassIVE MSstatsPTM analysis container RMSV000000358, which was also used as input to the MSstatsPTM workflow.
This ensures compatibility of the results with those in the MSstatsPTM paper. In this container, the analysis file MSstatsPTM_USP30_Analysis. R can also be found, which was used for the MSstatsPTM results. All steps of the workflow were followed as described above. The normalisation step made use of the available PTM dataset, given the lack of a non-enriched counterpart. Because each condition consists of two biological replicates which in turn consists of two technical replicates, we used the msqrob function with a mixed model as input. The results of both analyses were used to generate line plots with input as well as our normalised PTM level-data and the estimated effects for each condition. The detailed model results in the MSstatsPTM model object allowed us to inspect the model output for each PTM and protein as well as those for PTM upon correction for protein.
Phospho dataset The human phosphorylation datasets consist of 47 samples from condition A and 43 from condition B. Two aliquots were processed for each sample: one dedicated to total proteome analysis, and the other one to the phosphoproteome analysis. The main sample preparation steps were identical for proteomics and phosphoproteomics apart from the additional phosphopeptide enrichment step. Briefly, MeOH precipitation was performed on all samples and protein pellets were resuspended with 0. 1% RapiGest TM surfactant (Waters). Either 20 µg (proteomics) or 100 µg (phosphoproteomics) of samples were subjected to overnight trypsin/lysC (Mass Spec Grade mix, Promega, Madison, USA) digestion at 37°C with an enzyme:protein ratio of 1:25. Peptide samples were then incubated for 45 minutes at 37°C and centrifuged to remove RapiGest.
For total proteome analysis, collected supernatants were loaded on an AssayMAP Bravo (Agilent) for automated peptide clean-up using C18 cartridges. Desalted peptides were injected on a nanoAcquity UltraPerformance LC® (UPLC®) device (Waters Corporation, Milford, MA) coupled to a Q-Exactive Plus mass spectrometer (Thermo Fisher Scientific, Waltham, MA) and analysed using Data Dependent Acquisition (DDA). For phosphoproteomics, collected supernatants were loaded on an AssayMAP Bravo (Agilent) for automated Fe(III)-NTA phosphopeptides enrichment. Enriched samples were then analysed on a nanoAcquity UPLC devise (Waters) coupled to a Q-Exactive HF-X mass spectrometer (Thermo Scientific, Bremen, Germany) using DDA. Generated raw data files were searched against a database containing all human entries extracted from UniProtKB-SwissProt (25/08/2021, 20 339 entries) using MaxQuant (v.
1. 6. 17). The minimal peptide length required was seven amino acids and a maximum of one missed cleavage was allowed. For proteomics data, methionine oxidation and acetylation of proteins' N-termini were set as variable modifications and cysteine carbamidomethylation as a fixed modification. For phosphoproteomics data, serine, threonine and tyrosine phosphorylations were added as variable modifications. For protein quantification, the "match between runs" option was enabled. The maximum false discovery rate was set to 1% at peptide and protein levels with the use of a decoy strategy. Intensities were extracted from the Evidence. txt file to perform the following statistical analysis. All seven steps of the workflow were performed. The dataset can be found on PRIDE (PXD043476).
Mock analyses For the phospho dataset a mock analysis was included, that is an analysis where we only take one treatment arm of the data, so none of the PTMs (peptidoforms) are expected to be differential. We then assign the samples at random to a mock treatment with two levels and assess differential usage between the two conditions (mock vs control). In this way, correct control of the type I error by the statistical method can be assessed. Indeed, every PTM that is called as differentially abundant is a false positive. Hence, we expect the method to return uniform p-values. From the phospho dataset, only the samples from factor 1 condition B, and factor 2 condition y were withheld, i. e. 26 samples. Upon step 4, 13 out of the 26 samples were randomly assigned to condition "mock", the other 13 were assigned as condition "control".
Step 5 was then carried out by testing for a condition effect and the calculated p-values were retained. The randomisation to the mock treatment and step 5 in the analysis was repeated 5 times and a histogram was made for the p-values for each mock simulation. This mock analysis was done for different workflows: we assessed the effect of using robust regression in the modelling step, the use a non-enriched counterpart for normalisation and normalisation based on the enriched dataset, itself. Moreover, we conducted the analysis both on peptidoform as well as PTM-level.
PTM-level We first evaluated our method using the two computer simulations mentioned above. The first simulation consisted of 24 "perfect" datasets with no missing data and ten distinct peptidoforms carrying a specific PTM. Half of the datasets were simulated with a standard deviation of the difference in log-intensities between modified and unmodified peptidoforms of 0. 2, the other half had a standard deviation of 0. 3. The datasets differ in the number of replicates as well as in the number of conditions. Figure  3  shows the true positive rate (tpr, the fraction of the truly differentially abundant PTMs picked up by the method) -false discovery proportion (fdp, the fraction of false positives in the total number of PTMs flagged as differentially abundant) curve for simulation 1 for all 24 datasets.
As expected, both msqrob2PTM and MSstatsPTM perform better in datasets with lower variability and/or a higher number of replicates. Indeed, the true positive rate or sensitivity is higher for the same level of the false discovery proportion when the number of repeats increases while keeping the sd fixed (or when reducing the sd while keeping the number of repeats fixed). msqrob2PTM (solid line) clearly outperforms MSstatsPTM in all datasets (dotted line). Furthermore, MSstatsPTM in particular seems to have issues when the number of replicates is low. Indeed, in four out of six datasets with two replicates, the dotted line immediately veers right instead of up, indicating that non-DU PTMs are returned among the most significant features. This particularly affects datasets with higher variation (sd 0.
3). msqrob2PTM, however, does not suffer from a poor ranking of the PTMs for these four datasets and is still able to report (a few) true positive results at the 5% FDR level. Moreover, the fdp at the 5% FDR level for msqrob2PTM is close to 5% for most datasets, indicating a good control of false positives. Figure  4  shows the tpr -fdp curves for simulation 2 for all 24 datasets. As expected, the higher number of missing values induces a slight drop in performance overall. However, for the larger sample sizes the performance remains very good for msqrob2PTM. Again, msqrob2PTM uniformly outperforms MSstatsPTM and the fdp is close to 5% when adopting a 5% FDR threshold. For two datasets, we see that the far end of the tpr-fdp curve for msqrob2PTM veers straight up (two conditions, two replicates sd 0.
2 and sd 0. 3), which reflects msqrob2's inability to fit the models for a number of PTMs. This happens because these PTMs have too few observations to fit the models due to the missingness introduced in this simulation scenario. msqrob2PTM: differential PTM abundance and usage analysis For further comparison, ROC curves (true positive rate vs false positive rate) are shown in Supplementary figure  1  and 2. These plots give less weight to a few top-ranked false positives. Again, these ROC curves demonstrate superior msqrob2PTM performance. In supplementary table 1 and 2, the performance metrics (false positive rate, sensitivity, specificity, precision and accuracy) that were reported in the MSstatsPTM paper are also given for all datasets for comparison.
Biological spike-in dataset The design of the spike-in dataset (see also figure  2 ) is suboptimal to assess the performance of methods inferring differential PTM usage. This is because the spiked-in peptides and their corresponding protein abundance in the background proteome are not correlated as they would be in real experiments. Indeed, the latter does not contain the actual parent proteins of the spike-in peptides. Moreover, the E. coli proteins in mixes 3 and 4 induce loading differences present across the samples (see also supplementary figure  3 ), which brings additional normalisation issues. We illustrate these issues using ROC curves that compare the performance of different approaches: differential PTM abundance by adopting a conventional msqrob2 workflow directly on the summarized PTM-level intensities without normalisation (DPA-NonNorm), the same workflow upon normalisation with the median peptidoform logintensity (DPA), the default workflow for msqrob2PTM (default msqrob2PTM workflow assessing DPU), and MSstatsPTM (default MSstatsPTM workflow) (figure  7 ).
Every pairwise comparison between mixes is shown. Because all methods report many false positives for this dataset, the tpr-fdp plots quickly became unreadable (see supplementary figure  4 ). When comparing mix 4 to mix 1 (mixmix4), the log2FC after adjustment should be 0, hence, no method should report any differential PTMs. Indeed, this comparison is an internal control, and the ROC curves are expected to lie along the diagonal. Here, DPA-NonNorm and MSstatsPTM show the largest deviations from the diagonal. In the other comparisons, DPA always outperforms the other methods. Note that DPA assesses differential PTM abundance rather than differential usage as it does not normalise for parent protein intensity. This superior performance of the DPA method as compared to the DPA-NonNorm method indicates that it is very important to correct for technical variability resulting from the experimental design, i.
In other words, the design is not suited to benchmark the performance of methods developed to quantify differential peptidoform usage. However, the design is useful for assessing the performance of methods that quantify differential PTM abundance. This can easily be obtained with standard msqrob2 workflows, but is not returned by default by MSstatsPTM. However, because the msqrob2 suite builds upon the QFeatures architecture, the results of a DPA and DPU workflow can both be stored in the same object, thus providing more transparency and reproducibility across the workflows. For completeness, we also plotted the log2 fold changes for all PTMs in supplementary figures 5 and 6, which illustrate that both msqrob2PTM as well as MSstatsPTM provide good estimates for these. Ubiquitination dataset msqrob2(PTM) is capable of handling more complex designs that require mixed model analysis, as well as datasets that lack a non-enriched version of the dataset.
These two aspects apply to the ubiquitination dataset. Note that this is an experimental, biological dataset, and therefore does not come with a known ground truth. Despite the two abovementioned complexities, the standard msqrob2PTM workflow could find differentially abundant ubiquitin sites in most comparisons, except for the USP30_OE vs control comparison. However, table  2  shows that msqrob2PTM generally reports much fewer significant PTMs than MSstatsPTM. Upon closer inspection of the PTMs reported as significant by MSstatsPTM, it was discovered that this large discrepancy can be explained by several reasons.
Contrast First, both methods have a different way of dealing with missing data. Upon inspecting multiple line plots, we observed PTMs that were flagged as significant by MSstatsPTM despite having only one bio-repeat, or even only a single data point available in one of the conditions. In figure  8 , for instance, line plots are shown for two PTMs that are significant in MSstatsPTM when comparing the combination condition (Combo) versus the control condition (Ctrl), but not in msqrob2PTM. Notably, PTM O00154_K205 only presents PTM information for the first biological replicate, while PTM O00159_K0578 contains just one data point within the entire control condition. For these features, msqrob2 therefore did not return a model fit. When examining the results more closely, we noticed that MSstatsPTM uses three different models to fit the data and that the model choice is based on the available data points for each PTM (see UbiquitinationBioData exploration of results file on https://github.
Phospho dataset Two different workflows were employed for this dataset. The first workflow uses the nonenriched counterpart dataset to normalise for differences in protein abundance, while the second workflow only used the enriched dataset, also for the normalisation step. It is important to note that two distinct instrument platforms were used to analyse the total proteome and phosphoproteome samples. The chromatographic conditions were identical as well as the MS instrument geometry but two consecutive generations of Q-Orbitraps were used (Q-Exactive Plus versus Q-Exactive HF-X). This partly explains the observed heterogeneity between enriched and non-enriched datasets. Indeed, we observed a substantial proportion (approximately 25%) of proteins present in the enriched dataset that were absent in the non-enriched one.
This led to some PTMs that could not be normalised, which we opted to exclude from subsequent analysis in workflow 1. Both workflows involved testing multiple contrasts based on two factors: condition (A or B), and subset (x or y). In the first workflow (utilising both datasets), 31 unique differential PTMs were found, of which 25 phosphorylations. Most of these PTMs exhibited significant downregulation in condition A compared to B within subset y. In the second workflow (using only the enriched dataset), fourteen unique significant PTMs were identified, of which eight phosphorylations. The majority of phosphorylations showed significant differential usage between condition A and B within subset y and/or exhibited significant differential usage between condition A and B averaged over subsets x and y.
Supplementary tables S5 and S6 provide detailed results. Interestingly, the results differ between the two workflows. Of the 31 PTMs identified in workflow 1, ten were also found in workflow 2. Instead of solely focusing on significant PTMs, our method is capable of detecting differentially used peptidoforms as well. For this dataset, the first workflow detected twelve peptidoforms as differentially abundant, predominantly showing downregulation in condition A for subset y. In the second workflow, which lacked a global profiling dataset, seven significant peptidoforms were detected across the different comparisons. LPIVNFDYS[Phospho (STY)]M[Oxidation (M)]EEK was picked up as DU by both workflows and is particularly interesting, because both PTMs present on this peptidoform are also returned as significant in the differential PTM usage analysis.
However, the significance of the phospho-PTM might be an artifact triggered by the presence of additional oxidation in one of its underlying peptidoforms. ]EEK) and its corresponding PTMs for the phospho dataset. At the top, the significant peptidoform is depicted in pink. In green is the PTM occurring on that peptidoform, in this case phosphorylation. In grey any other peptidoform carrying that same PTM, and in yellow, the PTM intensity value as estimated by the model. The PTM is represented by two peptidoforms that roughly follow the same pattern, resulting in a PTM that resides in the middle. At the bottom we see the other PTM occurring on that peptidoform, the oxidation. No other peptidoform carries that same modification, resulting in perfect overlap between the line of the significant peptide and that of the PTM.
Here, it is possible that the oxidation is only significant because the phosphorylation is. Indeed, the driving force of the significance of this particular peptidoform could be coming from the phosphorylation (which has two associated peptidoforms). Note that, while these particular line plots were derived using the workflow without a non-enriched dataset, the corresponding plots from workflow 1 are extremely similar. Some PTMs are also significant because they enable aggregating evidence over multiple nonsignificant peptidoforms that all have a similar expression pattern. An example of this can be seen in figure  11  for sp|P10451|OSTP_HUMAN (Phospho (STY)) 280.
Mock analyses As the phospho datasets are biological experiments, the ground truth is unknown. Therefore, we cannot assess the performance of each method. We also do not know if the method provides reliable false positive control. To assess if our workflows provide good type I error control for the case study, we therefore perform a mock analysis. In particular, we introduce a factor for a non-existing effect, implying that all features that are returned significant upon testing for this factor are false positives. Here, we focus on subset y from condition B, so that ample samples remain. When the method provides good false positive control, the p-values upon assessing the mock effect will be uniform. The p-value distribution for the workflow that only uses the enriched dataset is given in Figure  12.
The top panels show the results for the PTM-level analysis and the bottom panels for peptidoform analysis. Both workflows with and without robust regression provide fairly uniform p-values. Supplementary figures 7-10 show similar plots for four other random mock datasets, showing consistency of performance. We did a similar mock analysis for the workflow that uses the non-enriched dataset for usage calculation (Figure  13 ). The workflow on peptidoform level using robust regression showed a slight increase in low p-values, which is also observed in some other random mock datasets (Supplementary Figures  11 12 13 14 ). The remaining workflows generated fairly uniform p-values for all random mock datasets (Figure  13  and Supplementary Figures  11 12 13 14 ). We therefore did not adopt robust regression for the peptidoform analysis.
Therefore, the choice between DPA and DPU depends on the specific research question at hand, or they can both be performed to complement each other. Through analysis of simulated and biological datasets, we have demonstrated that our workflows improve upon the state-of-the-art MSstatsPTM. We showed the advantage of first normalising the peptidoform intensities by the parent protein abundance before conducting the differential analysis. In this way, we can immediately model the usages as opposed to MSstatsPTM that estimates the fold changes for the PTM and protein values separately before differencing these to estimate DPU. Indeed, the peptidoform and protein values from the same sample are correlated, which is explicitly accounted for in our DPU workflow but is ignored by MSstatsPTM. We showed for the latter method that this can lead to artifacts in the estimated fold change for some PTMs upon correction for the fold change in the parent protein.
In general, a standard user is not fully aware of these issues, and the subtleties of interpretation that these require. In contrast, our msqrob2PTM workflow emphasises transparency and reproducibility. While this choice may lead to some PTMs that cannot be estimated using the default workflow, it does ensure that users are fully aware of what was modelled for each PTM. Moreover, we feel that PTMs for which no results are returned due to missingness require the intervention of a skilled data analyst to develop tailored solutions to infer on differential abundance and/or usage; solutions that are moreover supported by the msqrob2 universe. Indeed, we showed that automatic approaches can lead to biased results, and especially in experiments with more complex designs. These differences in normalisation approach and design concept elucidate the variations in performance across the different datasets that were used in our benchmark.
In the simulated datasets, msqrob2PTM capitalises on the within-sample correlation between peptidoforms and proteins that is present in the data, resulting in superior performance compared to MSstatsPTM. However, in the spike-in dataset, where this correlation is absent due to its unrealistic design, the default msqrob2PTM workflow exhibits similar performance to MSstatsPTM. However, for this dataset we show that our workflow for assessing differential PTM abundance analysis uniformly outperforms both the msqrob2PTM and MSstatsPTM workflows assessing differential PTM usage. Indeed, the spike-in study is suited for assessing the performance on differential PTM abundance rather than on differential PTM usage, as the spiked PTMs were not correlated to their corresponding protein in the background.
In the biological ubiquitination dataset, the high amount of missing data, and the absence of a global profiling dataset leads to a high number of PTMs that cannot be fitted with the required model. MSstatsPTM will then resort to other, simpler models that are often suboptimal or even mismatched, while msqrob2PTM will simply not return results for these PTMs, leading to a lower number of reported significant PTMs. These datasets bring to attention a broader issue in the field, specifically the scarcity of suitable datasets for accurately assessing Differential Peptidoform Usage (DPU). When designing such experiments, it is favourable to incorporate a global profiling dataset along with an adequate number of biological replicates. This comprehensive approach not only enables a more thorough evaluation of DPU but also enhances statistical power, yielding more reliable and robust results.