dataesr/scientific-paragraphs-categorization huggingface dataset:
These simulations were then exploited to generate 3D pharmacophores. Following an optimization with a training dataset, five features were enough to model XIAP-BIR3 synthetic ligands binding to two hydrogen bond donors, one hydrogen bond acceptor and two hydrophobic groups. The correlation between pharmacophoric features and computed MM-PBSA free energy revealed nine residues as crucial for synthetic ligand binding: Thr308, Glu314, Trp323, Leu307, Asp309, Trp310, Gly306, Gln319 and Lys297. Ultimately, and three of them seemed interesting to use to improve XIAP-BR3 versus cIAP-BIR3 selectivity: Lys297, Thr308 and Asp309.
Moreover, it plays an important role in chemoresistance and has therefore become an important target for the malignancy treatment  [4]. XIAP belongs to the inhibitor of apoptosis proteins (IAPs) family, which consists of eight different proteins. They all share the zinc-binding baculovirus-IAP-repeat (BIR) domain comprising around 70 amino acids. Among IAP family members, XIAP is the only antagonist binding affinities, using the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) approach. MM-PBSA is widely employed to predict ligand-receptor free energy thanks to its efficiency and reasonable accuracy. Then, we exploited the molecular dynamics simulations that were producing the best agreement with experimental binding data in order to build a 3D pharmacophore for a synthetic XIAP-BIR3 antagonist.
MM-PBSA Binding Free-Energy Prediction To reach our goal, we started by setting up a method for rapid and reliable predictions of synthetic ligand-binding affinities to the XIAP-BIR3 domain. To assess the predictive performance, we selected four XIAP-BIR3 synthetic ligands for which the X-ray structures of their complexes with XIAP-BIR3 (residues 248-352) were known, and the binding constants were measured experimentally. At the time of our work, 36 3D structures with XIAP-BIR3 were available in the PDB database  [19] ; 4 NMR solution structures and 32 structures were solved using X-ray diffraction. Among them, only 16 were co-crystallized with a non-peptide or non-peptidomimetic ligand. For our study, we selected four complexes from sixteen available ones with the following PDB IDs: 5C7C  [20] , 5M6M  [21] , 5OQW  [22]  and 5M6L  [21].
The choice of these complexes was based on the ligand structural diversity/similarity, as well as on the difference/similarity in their binding activities to the XIAP-BIR3 domain determined by fluorescence polarization assay (see Figure  1 ). simulation strategies. Our strategies included a ranking of relative XIAP-BI3 antagonist binding affinities, using the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) approach. MM-PBSA is widely employed to predict ligand-receptor free energy thanks to its efficiency and reasonable accuracy. Then, we exploited the molecular dynamics simulations that were producing the best agreement with experimental binding data in order to build a 3D pharmacophore for a synthetic XIAP-BIR3 antagonist. The pharmacophore performance was assessed and optimized with a testing chemolibrary.
MM-PBSA Binding Free-Energy Prediction To reach our goal, we started by setting up a method for rapid and reliable predictions of synthetic ligand-binding affinities to the XIAP-BIR3 domain. To assess the predictive performance, we selected four XIAP-BIR3 synthetic ligands for which the X-ray structures of their complexes with XIAP-BIR3 (residues 248-352) were known, and the binding constants were measured experimentally. At the time of our work, 36 3D structures with XIAP-BIR3 were available in the PDB database  [19] ; 4 NMR solution structures and 32 structures were solved using X-ray diffraction. Among them, only 16 were co-crystallized with a non-peptide or non-peptidomimetic ligand. For our study, we selected four complexes from sixteen available ones with the following PDB IDs: 5C7C  [20] , 5M6M  [21]  5OQW  [22]  and 5M6L  [21].
The same trend was also detected for the predictions with the entropic temp calculated using the IE method, except for the simulations without additional parametrization, for which the lengthening of the simulations increased the prediction accuracies this time. Nevertheless, in the case of predictions with the entropic term computed using NM, we did not succeed in deriving rules. Next, we tested the prediction performance protocols on a ligand not belonging to the synthetic ligand category. We chose to use the AVPI tetrapeptide of Smac. To be able to compare the predicted free energies for the synthetic ligands to that of AVPI, the CGENFF force field was also applied to parameterize the AVPI tetrapeptide. As the XIAP-BIR3 3D structure co-crystallized with only AVPI tetrapeptide is not available in the PDB databank, we generated this complex by using the docking (see Supplementary Materials Figure  S4 ) and submitted it to four MD simulation protocols, like for the previous complexes with synthetic ligands.
Construction and Validation of a 3D Pharmacophore of Synthetic XIAP-BIR3 Inhibitors In the final step, we built a representative 3D pharmacophore for the synthetic ligands binding to the XIAP-BIR3 binding groove. Indeed, the pharmacophore approach has become one of the important ligand-based tools in in silico methodologies for drug research and development process. The interest of this approach is to highlight essential motifs for the ligand-target interaction, which can then be applied to guide the design of new ligands. The LigandScout software  [28]  used in this study allowed us to generate the 3D pharmacophore either using only the ligand structures (ligand-based approach) or using the target structure in addition to the ligand structure (structure-based pharmacophore)  [29]. In our study, in order to take into account the essential structural motifs for the interaction with the target, we generated XIAP-BIR3 pharmacophores by using the structure-based approach.
Among the eight pharmacophoric features, Lig-andScout annotated four hydrophobic (H) and one aromatic-hydrophobic features (HAr), one hydrogen bond acceptor (HBA) and two hydrogen bond donors (HBD). The identified features of the pharmacophore reflect the ligand interaction with the XIAP-BIR3 binding site through Lys297, Thr308, Asp309, Glu314, Trp323 and Tyr324 (Figure  4 ). Thr308, Asp309, Glu314 and Trp323 residues already emerged as crucial from our free-energy decomposition per residue (Figure  3 ). In addition, LigandScout software also highlighted Lys297 and Tyr324 as important anchoring residues by positioning the hydrophobic features in their proximity. For Lys297 and Tyr324, we detected via our MM-PBSA analyses their dominant contribution non-polar one, which is in good agreement with the proposed hydrophobic characters of these residues.
Interestingly, comparing XIAP-BIR3 residues to cIAP1-BIR3 and cIAP2-BIR3 residues revealed crucial residues that bind to synthetic ligands (see Figure  5C ). Only three residues are not conserved between XIAP and cIAP1-2: Lys297/Asp, Thr308/Arg and Asp309/Cys. It can therefore be concluded that targeting these three residue side chains could increase the ligand selectivity for XIAP versus cIAP1-2. Among them, the Lys297 side-chain interaction through pharmacomodulations of the G2 group seems to be the most promising way to improve the selectivity. A pharmacophore for XIAP-BIR3 inhibitors has already been published in the literature  [30]. Opo and co-workers also generated their pharmacophore by using LigandScout software, but from the 5OQW X-ray structure, one of the structures used in our MM-PBSA study.
Structure Preparation During this study, X-ray structures of the following complexes with XIAP-BIR3 retrieved from the PDB database were used (Table  4 ): 5C7C  [20] , 5M6M  [21] , 5OQW  [22]  and 5M6L  [21]. Then, to calculate the experimental value of ligand binding free energy ∆G exp , the following equation was applied: ∆G exp = RT × ln IC 50 (1) where the gas constant value was R = 1. 985 8775 × 10 -3 kcal•K -1 •mol -1 and T = 303, 15 K. As the XIAP-BIR3 complex with Smac AVPI tetra peptide was not solved to date, to prepare this complex, we chose to apply the docking strategy. The AVPI tetra peptide, which was retrieved from the PDB structure of XIAP-BIR3 co-crystallized with full SMAC (PDB ID: 1G73  [12] ), was docked into the X-ray crystal structure 5C7C of XIAP-BIR3, without waters and bounded synthetic ligand.
5 kcal/mol/Å 2 force constant for sidechain atoms), followed by 10,000 steps of minimization without any restraints. Next, the minimized systems were heated to 303. 15 K, and the dynamics were temperatureequilibrated during 250 ps via heating reassignment under NVT conditions, with harmonic restraints applied to heavy atoms. Finally, the systems ran freely for 50 ns under NPT conditions, with a time step of 2 fs for the simulations without HMR additional parametrization and 4 fs with HMR parametrization. Langevin dynamics with a damping coefficient of 1 ps -1 was used to maintain the system temperature, and the Nosé-Hover-Langevin piston method was used to control the pressure at 1 atm. Production trajectories were saved every 20 ps, and subsequent analyses were performed using the CHARMM program version c40b2  [41].
G polar was determined by solving the Poisson-Boltzmann equation, whereas the G non-polar term was estimated from a linear relation to the solvent accessible surface area (SASA), using the equation 0. 00542 × SASA + 0. 92. The value of G polar and G non-polar in each trajectory snapshot was calculated as follows: G polar = G PB Complex -G PB Recptor -G PB Ligand (3) G non-polar = G SASA Complex -G SASA Recptor -G SASA Ligand (4) To calculate all contributions to the MM-PBSA CHARMM, homemade scripts were applied. For the last term, entropic contribution, three approaches were tested: firstly, (i) the entropic term was omitted, and then the entropic term was approximated using either (ii) normal mode analysis of harmonic frequencies calculated at the Molecular Mechanics level  [25]  or using the recently published approach (iii) interaction entropy that investigates the entropy change upon binding  [26, 27].
3D Pharmacophore Building The 3D pharmacophore was generated using the 3D structure of human XIAP-BIR3 co-crystallized with a synthetic inhibitor (PDB ID code 5M6L  [21] ) and using the structural data resulting from its MD simulations with additional HMR parametrization. LigandScout software  [28]  was applied for the detection and interpretation of crucial interaction patterns between XIAP-BIR3 and the ligand. LigandScout extracts and interprets ligands and their macromolecular environment from PDB files and automatically creates and visualizes advanced 3D pharmacophores. Besides the features representing the interaction of the ligand with the target (as hydrogen bond donor, hydrogen bond acceptor, hydrophobic group, etc. ), we also generated our pharmacophore exclusion volumes. In order to compare the pharmacophore performance during the optimization rounds, we referred to the evaluation metrics, such as sensitivity, specificity, enrichment factor and area under the ROCs (Receiver Operating Characteristics) curve.
The AUC value between 0. 5 and 1 is mandatory for a model to be validated (AUC = 1 perfect prediction; AUC = 0. 5 random prediction). These tests are defined by the following equations  [45] : where TP is the number of active ligands aligned with the pharmacophore, TN is the number of discarded inactive ligands and FP is the number of inactive ligands aligned with the pharmacophore. Moreover, Htot = (TP + FP) is the total number of aligned ligands, A is the number of active ligands in the initial dataset, I is the number of inactive ligands in the initial dataset and D is the total number of ligands in the initial dataset. In general, an efficient screening should have high values of sensitivity, specificity, enrichment factor and area under the ROC curve. During the screening, we applied an alignment of all features on the ligands and the extraction mode based on the most stable ligand energy conformation.
Testing Chemical Library To test our 3D pharmacophore performance, we created a validation dataset. To build it, we first detected all compounds annotated with XIAP-BIR3 biological activities in CHEMBL database  [46, 47]. From this dataset, the redundant ligands and the peptide ligands or peptide-mimetic ones were deleted. This resulted in a dataset of 487 synthetic XIAP-BIR3 ligands. Then, it was necessary to determine a cutoff (threshold) to classify the ligands in terms of biological activity (active/inactive). As the number of ligands available for XIAP-BIR3 was limited, we applied the approach that is usually used in the literature for targets that are not perfectly annotated  [48]. In this strategy, a pIC 50 (pIC 50 = -log 10 IC 50 ) value of 6 is usually used as a cutoff. All molecules with a pIC 50 value greater than 6 are considered active, and those with a value under 6 are considered inactive.
Moreover, molecules with a pIC 50 value between 4. 5 and 6 are excluded from the dataset (see Table  2 ). This led us to a dataset with a ratio of active to inactive ligands of 1/1. 3 (i. e. , 173 active ligands with respect to 218 inactive ligands). This ratio was not sufficient to obtain statistically significant values during the screening and correctly assess its performance. Thus, an increase in the number of inactive compounds in our dataset was necessary to improve the screening performance, and to do so, we included XIAP-BIR3 decoys from the DUD-E database  [48]  (compounds were assumed to be inactive and built from fragments of active molecules). For our target, there were 5199 decoys available that we added to our starting dataset, resulting in a library of 5686 compounds. This library was then used to evaluate the performance of our 3D pharmacophores (see Table  2 ).
Data Availability Statement: Not applicable. Acknowledgments: This work was supported by a doctoral fellowship (MENRT) held by M. F. The present work was performed thanks to the generous computing resources provided by the CRIANN (Normandy, France). The authors gratefully acknowledge the European Community (ERDF: European Regional Development Fund) for the molecular modelling software. We are also grateful to Thierry Langer and Nathalie Lagarde for helpful discussions during this work.
We argue that this approach is suboptimal as it fails to leverage the inherent correlation between the parent protein and PTMs or peptidoforms, i. e. , a specific peptide with its corresponding modifications. Additionally, the separate modelling and adjustment process in MSstatsPTM can artificially amplify small differences. This phenomenon is demonstrated in figure  1. Here, we can see a PTM for which the PTM intensities closely mimic the protein intensities across the samples. Although no significant differences are observed at the PTM or protein level when comparing the "Combo" and "Ctrl" conditions in the dataset, the adjustment inflates the difference, causing MSstatsPTM to return a significant PTM. Hence, in msqrob2PTM we employ a different normalisation strategy that directly accounts for this correlation between peptidoform and protein.
Additionally, we will not limit ourselves to the analysis of the PTMs. Indeed, our method can manage the analysis of peptidoforms as well. In many studies, each distinct PTM will likely not be characterized by a myriad of peptidoforms. It is therefore possible that a significant PTM effect can be attributed to only one or two strongly significant associated peptidoforms, which may be significant for another reason, i. e. a different PTM occurring on that (those) peptidoform(s). We think it is crucial that potential users thus do not restrict their analysis to the PTM alone, but also assess the individual peptidoforms that carry the specific PTM. We here present a statistical, R-based workflow, based on the msqrob2 R package  (16) , to carry out differential abundance as well as differential usage analysis at the peptidoform and PTM level.
We apply this workflow to simulated datasets, a spike-in study, and to biological datasets, and use these to compare our method to MSstatsPTM. We show that our approach does not suffer from the artifacts that are introduced by uncoupling the within-sample correlation between PTM and parent protein, while maintaining good sensitivity and FDR control. The approach is freely available and can be consulted on https://github. com/statOmics/msqrob2PTMpaper.
Workflow The general workflow for the differential abundance analysis on PTM and peptidoform level was developed in R(  17 ) (version 4. 2) and is mainly based on two R packages: msqrob2 (https://www. bioconductor. org/packages/release/bioc/html/msqrob2. html, version 1. 6. 0) and QFeatures(18) (https://www. bioconductor. org/packages/release/bioc/html/QFeatures. html, version 1. 8. 0). QFeatures provides an infrastructure to store and manage mass spectrometry data across different levels (e. g. peptidoform and protein level) whilst keeping links between the levels where possible. For each preprocessing step a novel, linked assay is constructed. In this way, the original data is not overwritten, and preprocessed data can be traced back to its origin. msqrob2 is a package with updated and modernised versions of the MSqRob  (16)  and MSqRobSum  (19)  tools and builds upon the QFeatures class infrastructure.
e. 𝑦 𝑖,𝑝,𝑃 * = 𝑦 𝑖,𝑝,𝑃 -𝜇 𝑖,𝑃  (1)  With 𝑦 𝑖,𝑝,𝑃 * the normalised log2-transformed intensity for peptidoform p in sample i with parent protein P, 𝑦 𝑖,𝑝,𝑃 the log2-transformed intensity for peptidoform p in sample i with parent protein P before normalisation and 𝜇 𝑖,𝑃 the summarised intensity for protein P in sample i. It is possible to calculate the summarised protein intensity value directly from the PTM dataset itself. However, when the experiment includes both an enriched and non-enriched (global profiling) dataset we recommend using the non-enriched dataset to calculate the summarised protein values. Of note, steps one and two should also be applied to the non-enriched data.
Summarisation of peptidoforms to PTM level For each unique PTM (i. e. unique protein -modification -location combination), we need a summarised intensity value per sample. This is done by taking a subset of the dataset with all peptidoforms containing a specific PTM and summarising all corresponding intensity values into one value per sample. When peptidoforms contain multiple PTMs, these are used multiple times. Here we apply robust regression using the robustSummary function in the MsCoreUtils(20) R package by default to summarise the peptidoform level data at the PTMlevel. In this way, we obtain an intensity assay on the PTM level. This assay can then be added to the existing QFeatures object.
PTM level analysis We use the functionalities of the msqrob2 package for this step. Msqrob2  (16, 19, 21, 22)  provides a robust linear (mixed) model framework for assessing differential abundance in proteomics experiments. To assess differential abundance on the protein level, the workflows can start from raw peptide intensities or summarised protein abundance values. The model parameter estimates can be stabilized by ridge regression, empirical Bayes variance estimation and robust M-estimation. Here we assess differential abundance on the PTM level by first summarising peptidoform expression values (step 5). When one predictor (e. g. condition) is present in the dataset, we perform an msqrob analysis on PTM intensities with the following model: 𝑦 𝑐𝑠 = 𝛽 0 + 𝛽 𝑐 𝑐𝑜𝑛𝑑𝑖𝑡𝑖𝑜𝑛 +𝜀 𝑐𝑠 With 𝑦 𝑐𝑠 the summarised log2-transformed PTM intensity in sample s of condition c, 𝛽 0 the intercept, and 𝛽 𝑐 𝑐𝑜𝑛𝑑𝑖𝑡𝑖𝑜𝑛 , the effect of a condition c.
Data Our novel msqrob2 workflow is tested and benchmarked to MSstatsPTM using two computer simulations developed by the MSstatsPTM team, the spike-in dataset from the MSstatsPTM paper, and data from two real experiments. The computer simulations were specifically developed for testing differential PTM workflows, and also allowed us to directly compare our method to MSstatsPTM. The first simulation produced "perfect" datasets with no missing values and many modified features per PTM, while the second simulation incorporated missing values and limited modified features, producing more lifelike datasets. The spike-in dataset consists of fifty human ubiquitinated peptides that were spiked into four background mixtures in known amounts. Hence, the true log-fold changes and identity of the truly changing PTMs are known.
Computer simulations We used the two computer simulations from the MSstatsPTM team that were found on https://github. com/devonjkohler/MSstatsPTM_simulations/tree/main/data (simulation1_data. rda and simulation2_data. rda). The first simulation consists of data without any missing values, while in the second simulation, missing data is introduced. For each simulation, 24 datasets were created with different experimental designs and intensity variance. In each dataset, 1000 PTMs were simulated. Half of the PTMs were simulated to have a fold change between conditions. However, of the half with differential fold changes on the PTM level, 250 could be confounded with differential fold changes of the parent protein. For further details on the creation of the datasets, we refer to the MSstatsPTM paper  (12)  and to their GitHub page.
Both simulations contain an enriched PTM dataset as well as its non-enriched protein counterpart. From each of the 24 datasets, the FeatureLevelData was extracted from the PTM and the protein dataset. These two datasets were then used as input to the workflow and all seven steps were followed. The protein dataset was used for the normalisation step. Because it is known which PTMs are differentially abundant and/or differentially used, we can readily evaluate the performance of a method in terms of the false positive rate (fpr), sensitivity, specificity, precision and accuracy, and true positive rate (tpr) -false discovery proportion (fdp) plots. Note, that tpr is the fraction of the truly differentially abundant PTMs picked up by the method and fdp is the fraction of false positives in the total number of PTMs flagged as differentially abundant.
On the tpr-fdp plot we also indicate the observed fdp at 5% FDR cut-off, which is expected to be close to 5%. We compared our results with the results obtained with the MSstatsPTM method, on their GitHub page https://github. com/devonjkohler/MSstatsPTM_simulations/tree/main/data (adjusted_models_sim1. rda and adjusted_models_sim2. rda) and included these in the tpr-fdp plots.
Spike-in dataset The MSstatsPTM team also developed a biological spike-in dataset with known ground truth to test their approach. Fifty human ubiquitinated peptides were spiked into four background mixtures consisting of human and E. coli proteins in different amounts. These four mixtures represent four different conditions and for each, two replicates were created. An overview of the experimental design can be seen in figure  2. Because the amount of spiked-in peptides is known, the true log-fold changes between the conditions is known and it is possible to assess whether the method can pick up these fold changes, and if these fold changes differ from the fold change of the corresponding protein in the background. Note, however, that as opposed to real experiments, the ubiquitinated peptides in the spike-in study are not correlated to their corresponding protein in the background.