MSProfileR: An Open-Source Software for Quality Control of MALDI-TOF Spectra

2. Details on MSProfileR Tool

2.1. General organization of the “MSProfileR_v1.0” Tool

The MSProfileR tool is a Shiny application assessing the conformity of mass spectra profiles for future analysis. It was developed in R programming language version 4.3.2 and in the RStudio environment version 2023.10.31 (R Core Team, 2020). The general workflow of the application was divided into five major steps: 1) the data loading tool that includes spectra import, 2) the preprocessing tool including successive steps to control the conformity and quality of spectra, 3) the processing tool including peak detection and MS spectra classification, 4) the spectra annotation tool allowing to add complementary metadata to each spectrum and finally 5) the output tool generating several files including a report (Figure 1). Each step of the MSProfileR was composed of several modules performing each one or more tasks. During the pipeline process, the user controls and adjusts the parameters of the MS spectra analysis. However, it is possible to realize the entire workflow without any intervention by using default parameters or loading parameters defined in a previous analysis. To make MSProfileR easier to use, a graphical user interface (GUI) was created using the Shiny R package version 1.7.2. A wrapping of existing R packages, available in Comprehensive R Archive Network (CRAN) repositories was done. Shiny R package allows to build an interactive web application straight from R. The different steps and modules from MSProfileR are described below.

2.2. Data Loading

Data loading, which is the entry point to MSProfileR, contains a module for reading MS spectra, as well as an optional module for selecting and loading the setting parameters of a previous analysis.

2.2.1. Module of the MS Spectra Loading

Functions from the MALDIquant and MALDIquantForeign [16] packages were used for retrieving MALDI-TOF spectra on R system. These packages recognize MALDI-TOF mass spectra file types and load them automatically [16]. Different file formats (txt, tab, csv, fid, ciphergen, mzXML, mzML, imzML, analyze, cdf, msd) can be uploaded, in an automatic detection way, from a folder containing one or many spectra. However, as only MS spectra from Bruker instrument were available, the assessment of MSProfileR tool to these file format. The spectra files are loaded recursively, from which are retrieved, the m/z values and intensity values, in the x- and y-axis, respectively, plus several data related to each spectrum such as the spectrum name, the date of acquisition and the sample name. The sample name must be identical for the spectra of the biological replicates loaded by MSProfileR. Effectively, this parameter will be used during the preprocessing to generate automatically an average spectrum representative for the replicates originated from the same sample.

To download the spectra profiles correctly, the level number of each main folder path must be defined before loading the data. This path is the number of folders existing between the sampleName and the binary files (fid or aqu), which may change between a direct classification of spectra by Bruker Daltonics Software and a manual spectra folder classification by the user if the sample name was not indicated before spectra acquisition. By default, the level of the spectra folders set on the MSProfileR front is four but the user can change it. Moreover, an ID number was attributed to each spectrum which will be useful to follow and to display spectra during the whole pipeline analysis.

2.2.2. Module of Parameter Loading

This module is used to upload all the parameters saved from MSProfileR in a previous analysis. The user can submit distinct datasets to MSProfileR while applying the same parameters, enabling reproducibility.

2.3. Preprocessing

The preprocessing step of the workflow is composed of several modules aiming to evaluate the quality of the spectra and to exclude inconsistent MS profiles (Table 1). All these steps are designed with functions of the MALDIquant package [16], with the exception of the quality control step which is based on the MALDIrppa package [17]. MSProfileR tool allows assessing the MS spectra quality throughout four successive steps: 1) trimming and conformity tests, 2) spectra cleaning, 3) quality control and visualization of potentially outlier spectra, and 4) replicate averaging. All these tasks are using raw mass spectra data. Spectra passing these quality control steps are then selected for processing.

2.3.1. Module of Trimming and Conformity Check

The first task allows fixing the low and high limits of spectra m/z values using the trim function and inspects conformity of the spectra. Without specification, the limits of m/z values of default parameters correspond to a range from 2 to 20kDa. The inspection of data conformity checks whether errors occurred during MS spectra data acquisition. It assesses 3 criteria:

(i)

Completeness: Are there any empty spectra (i.e., no data to load)?

(ii)

Missing values: Are there any spectra with irregular m/z values? Normally, the interval between two successive m/z values should remain equal or increase uniformly (i.e., no missing point or aberrant values).

(iii)

Spectra range: Does the length (ie, m/z values range) of spectra differ?

Empty or irregular spectra, or spectra with too different an m/z range, can compromise the subsequent steps and must therefore be detected. These steps controlled the consistency and conformity of the spectra. The list of non-compliant spectra was removed before continuing the analysis.

2.3.2. Module of Spectra Cleaning

The MS spectra that have passed the previous tests are then subjected to a series of transformations in order to standardize the data: 1) intensity transformation, 2) smoothing, 3) baseline correction, and 4) normalization. The intensity transformation stabilizes the variance and improves the graphical visualization of spectra, by reducing the scale effect flattening the low-intensity peaks. By default, this transformation is performed using the square root method [18], but three log transformations (log, log2 and log10) are also proposed [19,20]. Mass spectra typically contain a mixture of noise and signal. To increase the signal-noise ratio, algorithms are used to improve the measure of peak m/z intensity values and facilitating the peak detection. By default, the Savitzki-Golay-Filter algorithm is used to smooth the spectra with a half window size of 10, but the Moving average is also available [21]. The noise altering the base level of mass spectra is then corrected with the Sensitive Nonlinear Iterative Peak (SNIP) algorithm (number of iterations, n=100) [22] or with the TopHat, Convex hull or Median methods. Subtracting the baseline from the spectra facilitates profile comparison. In order to compare intensities at each m/z value between spectra, normalization is required. Three normalizations are available, the Total-Ion-Current (TIC) [23], the Probabilistic Quotient Normalization (PQN) [24] or the median methods. The TIC is applied by default to each spectrum and corresponds to setting the sum of the relative abundances of all the ions (m/z peaks) to one. All these transformations are based on the mass spectra processing functions (transformIntensity, smoothIntensity, removeBaseline, calibrateIntensity) of the MALDIquant package.

2.3.3. Module of Quality Control

The subsequent steps consist in detecting and filtering MS spectra considered as no conform (i.e., outliers). To this end, the screenSpectra function from the MALDIrppa package is used. It computes for each spectrum an atypicality score (A score). The calculation of this score is based on robust estimators [25]. The A score pointed out MS spectra for which peak intensity profiles diverge from the rest of the dataset. An upper and lower tolerance limits are returned at the same time by the function. MS spectra classified above the upper limit correspond to profiles of low intensity and poor resolution, whereas those below the lower limit have high peak intensity profiles [17]. To avoid the deletion of profiles harbouring high peak intensity, the lower threshold limit is not considered by default. Nevertheless, the user could reactivate it.

Several parameters of the screenSpectra function are available to find the best estimator used in the calculation of the A score, the threshold and the method from which the cut-off is estimated. Two estimators are available, the median absolute deviation (MAD) and the Q. The MAD is defined as the median of the absolute deviations from the median of the intensities. It is a robust estimator for fairly symmetric distributions. The Q estimator is an improved version of the MAD and is more efficient and adequate for non-symmetric distributions [25]. Here, the Q estimator was selected as default parameters. Five methods (i.e., Rousseeuw & Croux (RC), Hampel, extreme studentized deviation (ESD), boxplot and adjusted boxplot) are available to compute the threshold limit for the detection of outliers [26]. By default, the RC method is selected.

The spectra with an A score outside the upper limit are regarded as potentially faulty spectra. By default, the spectra detected as outliers are automatically rejected from further analysis. However, at this step, a graphical interface was added giving to the user the possibility to visualize all MS spectra and notably the outliers. Based on their visualization, the user could refine the decision to exclude spectra near the A score limits or to keep some outliers. Generally, low-quality spectra are characterised by rippled and indented profiles due to low signal to noise intensities which could be absorbed by the background noise. Following this step, all mass spectra for which profiles were considered as atypically and confirmed by the user, are counted and excluded from the next analysis steps.

2.3.4. Module of Spectra Averaging

This module generates and proposes an archetype for MS spectra replicates. When available, the MS spectra replicates that passed the quality control steps were then averaged. An average spectrum is calculated for each sample by applying the mean, the sum or the median methods of its replicates, using the averageMassSpectra function from the MALDIquant R package. Then, one archetype spectrum of the replicates per sample is further analysed. The absence of replicates does not affect the workflow of the MSProfileR application. Spectra without replicates are the archetype profile of the sample. Samples with or without replicates could be downloaded and analysed concomitantly.

2.4. Processing

The processing tool consists in performing, successively, the peak detection, the spectra alignment until the spectra visualization, classification and the creation of an intensity matrix for further analyses.

2.4.1. Module of Peak Detection

This module aims to detect MS peaks excluding background noise for each spectrum taking into account homogeneity of the number of MS peaks detected in the entire dataset. The peak detection is performed on the average spectra using the detect Peaks function from MALDIquant package. It screens each spectrum with a sliding window in order to find the local maximums. A local maximum is detected as a peak whether its intensity is greater than the noise level, estimated using the MAD by default or super Smoother algorithms, multiplied by a signal-to-noise ratio (SNR), with a half window equal to 20, by default. Indeed, the spectra profiles are discretized and uniquely detected peaks are kept. Each spectrum profile can be regarded as a peak list. By default, the SNR value was set at 2 corresponding to a good compromise for detecting as many peaks as possible, without detecting too much parasite peaks which belong to the background noise. The SNR value is a parameter which could be adjusted by the user. To determine the optimal SNR value, the number of detected peaks and their standard deviation were measured for each SNR value tested, ranging from 2 to 7. A box-plot was generated showing the number of peaks detected per SNR value for the entire dataset. The SNR value for which the number of detected peaks was homogeneous and stable among dataset was considered as the optimal trade-off.

2.4.2. Module of Spectra Alignment

Prior to align detected peaks from averaged spectra, the determination of reference peaks is required. As no spectrum was selected as reference, the reference Peaks function from the MALDIquant package is used. In such conditions, all peaks from the dataset with an occurrence upper than the threshold (peak frequency) set by default at 0.9, are used as reference peaks for peak alignment, with a default tolerance of 0.002. The strict (default) or relaxed parameters allow to consider all peaks or uniquely the most highest.

The alignment of spectra was performed with alignSpectra function of the MALDIquant R package. All the peaks detected among the average spectra per sample are aligned and calibrated using one function of the correction phase (warping function), lowess, linear, quadratic or cubic. The lowess method was set as default.

Despite this alignment, very close detected peaks, but not identical, are binned into one, in accordance with an applied tolerance. This tolerance corresponds to the maximal relative deviation of a peak position (m/z value) to be considered as identical and is set to 0.002 by default. Thus, the peaks for which the difference of the positions (m/z values) is lower than this threshold get their mass values equalized forming one peak. Two binning functions are available, the strict one, for which bins never contain two or more peaks of the same sample, and the relaxed one, for which multiple peaks of the same sample were combined in one bin. By default, the strict binning function was selected.

The filtering step consists to remove peaks infrequently detected. A minimum frequency per m/z detected peak should be defined. For example, by setting this parameter to 1, uniquely peaks common to all spectra were filtered. By default, this value was set at 0.2, corresponding to an elimination of all peak positions found in less than 20% of the dataset. The number of peaks detected is indicated and then an adjustment of this frequency is possible by the user.

2.4.3. Module of Spectra Clustering

This module consists to visualize detected peaks and to ordinate samples according to their profiles. An intensity matrix is generated, where average spectra per sample are represented and was used for the classification. In this way, a hierarchical clustering algorithm was applied to generate a clustered heatmap using an unsupervised algorithm, the hclust function, and the pheatmap package [27]. This package computes bootstrap values to indicate for each query the positioning relevancy in the dendrogram. The columns correspond to the peak positions in m/z and the rows to the sample spectra. The peak intensity values are represented with a rainbow scale.

2.5. Spectra Annotation

Optionally, it is possible to annotate averaged spectra by the loading of an annotation table. The information introduced in this table from each sample allows to enrich metadata associated to each averaged spectra. The annotation table tool is composed of one module. The module allows to retrieve the list of sample names which passed the preprocessing and the processing tests and are then directly included in the first column headed “sampleName” of a downloadable table (.xlsx file format). At this step, as all spectra replicates were already averaged, a unique “sampleName” was associated to each averaged spectrum. In this way, the folder from each averaged spectrum could be easily paired with data from the annotation table using the “sampleName”. When the user considers the information added sufficient, the annotated table was imported into MSProfileR tool. The module controls whether all averaged spectra were annotated and also whether some annotations were not paired with an averaged spectra. Although the importation of the annotation table is not mandatory, these annotations could be helpful for future comparative results of mass spectra profiles. This module uses the readxl and writexl packages to read and write .xlsx files.

2.6. Output

The output module is created to guarantee traceability for the user about all the process steps, the methods and parameters selected, but also to download and to save data that can be useful for upcoming analyses.

2.6.1. Reports

The report is implemented on R Markdown language [28], offering an interface to generate a downloadable document in pdf format including information about steps performed, methods selected and parameters applied. Moreover, the information about the spectra dataset, such as the name of the directory file, the number of spectra loaded and the results of pipeline analysis, and notably those about the conformity tests, the quality control steps with the list of excluded spectra (non-conform or outliers). In addition, all representations, such as tables, plots, gelviews and clustered heat map are included in the report.

2.6.2. Save Setting Parameters

This module allows, thanks to rjson package, to save in JSON format all methods chosen and associated setting parameters that could be easily downloaded for a future analysis.

2.6.3. Export Intensity Matrix

This module serves to export the data of the matrix intensity into a CSV file format. This matrix was extracted according to two criteria, on the ordinate are indicated the sample names, and on the abscissa the m/z values with respective.

2.6.4. Figure List (Plots, Graph, etc.)

This module exports all the figures generated during the analysis of the current dataset by the MSProfileR tool. All the figures are in SVG format and zipped into a unique file. This list includes ‘mass_bounds.svg’ for the defined extremities spectra in the trimming module, ‘screening.svg’ for the plot excluding the outliers in the quality control module, ‘peak_count_by_SNR.svg’ for the boxplot showing the number of peaks depending on the different SNR values in the peak detection module, ‘reference_peaks.svg’ for the generated peak list, ‘spectra_aligned.svg’ for the gelview displaying the alignment, ‘peaks_binned.svg’ for the gelview of the binned peaks, ‘peaks_filtering.svg’ for the gelview of the filtered peaks, the last four figures are illustrated in the spectra alignment module, ‘dendrogram.svg’ represents a clustering of sample spectra, ‘intensity_matrix_8x4.svg’,‘intensity_matrix_16x8.svg’ and ‘intensity_matrix_32x16.svg’, the figures of the intensity matrix can be loaded with three different dimension to improve the classification visualisation, which are illustrated in the clustering module.

2.6.5. Hierarchical Data Format Version 5 (HDF5) File

This module allows the storage in HDF5 file of the four main outputs of this analysis corresponding to the parameters selected, the averaged spectra, the intensity matrix and the annotation table. This HDF5 file is an open source format supporting a large, heterogeneous and complex data. The creation of this file is done using the hdf5r library, which generates a file format that can store and manage large, complex and heterogeneous data [29].

2.6.6. Module Downloading All Files

This last module was created in order that the user could download all files listed in the five preceding modules. It is then possible to download either some particular data or all the data generated by the analysis.

2.7. The User Interface (UI)

For unfamiliar users with the R language, a user-friendly interactive web interface was created making possible to exchange information following a user-machine interaction. The UI or graphical interface was based on the Shiny dashboard template (https://rstudio.github.io/shinydashboard/), using Shiny, a “package that makes it easy to build interactive web apps straight from R and Python”.(https://rdrr.io/cran/shiny/). The UI was compartmentalized in five tabs corresponding to the five steps of workflow, including 1) Data loading, 2) Preprocessing, 3) Processing 4) Annotation and 5) Output. The architecture of MSProfileR tool was organized in a modular way (Figure 2).

2.7.1. Web Interface Development Architecture

All R packages dependencies used for the creation of the MSProfileR tool are available in Bioconductor (https://www.bioconductor.org/) or Comprehensive R Archive Network (CRAN) repositories (https://cran.r-project.org/). The execution of MSProfileR tool modules is in cascade, with a specific modularity allowing the insertion of new functionalities and tasks without rewriting all R code. MSProfileR was made with the motivation to be used by everyone especially for individuals who have no programming skills. The installation of the R packages used for the tool building and the opening of the front of the interface, are obtained by installing the “MSProfileR” package on the R environment (version ≥ 4.3.0), which load all the needed dependencies for the launch of the application. The application launches into web browser.

2.7.2. UI of Data Loading Tab

The interface of the data loading tab was divided into two panels corresponding to the different modules from this tab (Supplementary Figure S1, Table 1). The MSProfileR tool starts with the loading of spectra. The level number of the binary files path sets at four by default should not be modified before spectra loading. The user selects the directory containing one or many spectra to analyse. The data related to each spectrum are gathered into a summarizing table, including the id spectra, a number generated by a programmatic way for each spectrum, the sample name, the replicate name and the acquisition date. The sample name will be used to detect sample replicates when available. The data from ten spectra are presented per table. Spectra are classified per loading order and the others are available by the page selection at the bottom of the table. The number of loaded spectra as well as the path are indicated at the top of the table. Thanks to the table, the user can control whether the sampleName and replicateName columns were correctly filled. If a shift was noticed, this problem could be rectified by adjusting the level number of the binary files path until the correct classification was obtained.

The second panel is optional. The user has the possibility at this step to upload setting parameters used in a previous analysis. This option allows to analyse the current spectrum dataset with the same parameters applied previously for another dataset.

2.7.3. UI of Preprocessing Tab

The preprocessing tool is represented on the interface by five panels corresponding to the four modules of this tab (Supplementary Figure S2, Table 1). It consisted of successive steps to evaluate the quality of the spectra to exclude inconsistent MS profiles. The quality control module was divided into two panels to distinct the automatic selection process from the optional manual validation of conform spectra. All these steps are interactive, the user could select a method and associated parameters when available, for dataset analysis.

The first panel consists to trim and to test the conformity of spectra (Supplementary Figure S2). The lower and upper limits of spectra range could be determined by the user and are plotted in a graduated scale of m/z (kiloDaltons, kDa). Following this step, conformity tests were automatically done and the results were presented in a table. The table rows are coloured in green if all spectra are compliant and in grey if some spectrum did not pass the test(s) successfully. Moreover, the number of spectra considered as non-conform per test was indicated and these spectra are automatically excluded from the rest of the analysis. Nevertheless, the user can keep these non-conform spectra by deselecting the checkbox untitled “Exclusion of non-conform MS spectra”.

The cleaning panel realizes successive steps to adjust spectra profiles including transformation, smoothing and normalization of intensity plus a baseline correction. The user can select different methods by radio button for each step and visualizes the result of treatment to each spectrum by a plot through a reactive window.

The panel of spectra quality control allows to visualize the results of the spectra screening according to the method selected, into a plot representing the atypical score (A score) of each spectrum indexed by a number corresponding to their loading order (i.e., id spectra). The dotted line represented the upper limit of the A score threshold above which spectra are considered as outliers. Numbers corresponding to atypical spectra are then coloured in red whereas those kept for the rest of the analysis are indicated in blue. Interestingly, it is also possible to detect spectra considered as outliers due to their low A score by deselecting the checkbox “Include spectra below the lower threshold”. In such condition, all spectra obtaining an A score outside the lower or upper limits will be considered as atypical. Nevertheless, based on our background, notably for MS spectra from entomological origin, protein profiles which obtained very low A score correspond generally to MS spectra with higher peak intensities. To avoid their classification as outliers, uniquely the upper limit is applied by default.

The validation of spectra classification as conform or atypical remains possible through the selection panel. This panel was composed of two boxes separating spectra according to their compliant status. One box lists the outliers (i.e., “Atypical spectra”), whereas the second box lists spectra kept for the next steps (i.e., “Selected spectra”). Spectra in each box are identified by their id spectra, sample name and A score. Using these box lists, the user can choose a spectrum which is plotted below. By default, atypical spectra are automatically removed from the rest of the analysis. However, after visual spectra checking, the user can decide to keep some “atypical spectra” or to remove some “selected spectra” by moving them from one box to another by using arrow buttons.

The last panel of this tab consists to average spectra replicates which succeeded in passing all the preprocessing steps. The user can select one of the three methods to perform the spectra averaging. Averaged spectra are presented in a table containing one representative spectrum per replicated sample, classified by the sample name. The number of averaged spectra was indicated at the top of the table. This averaging reduces the number of spectra to analyse in the processing part.

2.7.4. UI of Processing Tab

The interface of the processing tab is separated into five panels (Supplementary Figure S3). The first panel concerns the peak detection. The user can select one method of noise estimator and visualize the results of peak detection on a boxplot graph showing the variance of peak number detected per spectra in the dataset for each SNR value from 2 to 7. Based on this result, the user chooses the optimal SNR value for the current analysis using a graduated slider. A visualization of the peaks detected for each averaged spectrum is available. On the graph, the detected peaks are indexed in ascending order from highest to lowest intensities.

The next three panels correspond to the spectrum alignment modules, which were split in the interface to visualize and to control the outcomes of each step. Firstly, to realize spectrum alignment the requirement of reference peaks is compulsory. In this way, the user should select a method to obtain reference peaks among the spectra dataset based on two parameters, the minimum of peak occurrence (i.e., frequency) and the alignment tolerance. The number of peaks used as reference and their distribution in the range of m/z values are presented on a plot. Before to apply spectra alignment, an adjustment of parameters remains possible by the user to increase the number of the reference peaks. Then, the user launches the alignment process, by choosing one of the four methods of warping functions. The visualization of the alignment is presented by a gelview where abscissa and ordinate correspond respectively to peak position (m/z) and averaged spectra. Peak intensity is represented by a grey scale. Once the peaks are aligned, the next two panels concern the peak binning and peak filtering. For these two steps, different methods and parameters are offered to the user to reduce the number of peaks which will be retained for the creation of the matrix of peak intensities. The results of the binning and the filtering steps are illustrated by gelviews in the UI, and the total number of peaks (m/z values) from the dataset was also indicated at each step. Based on these information the user could decide either to adjust parameters or to continue the analysis.

The last panel consists to classify averaged spectra according to their profiles. In this way, a hierarchical clustering was performed and represented by a heatmap linked to a dendrogram. At this stage of the application, an intensity matrix for all the detected peaks is generated for future analysis. Each peak does not necessarily exist in all the spectra. By default, when peaks are missing, intensity values of the spectra are used to fill the intensity matrix. A checkbox allows this filling to be deactivated. In this case, when peaks are missing, NA values are inserted into the intensity matrix. In the heatmap, missing values appear in grey whereas the detected peaks were coloured in a rainbow scale according to the intensity values.

2.7.5. UI of Annotation Tab

This tab divided into three panels, aims to add sample information (metadata) to each averaged spectrum (Supplementary Figure S4). An annotation table (.xlsx) generated automatically by MSProfileR is downloadable in the first panel. Uniquely the first column “sampleName” of the annotation table is filled with the names of the samples which passed all previous tests. In the subsequent columns, the information of interest could be listed. The body part used, the arthropod family, its geographical origin, the method and mode of storing are some examples of added information which could be useful for the next steps of the analysis. Once the annotation table is enriched with metadata, it could be uploaded using the second panel of this Table Finally, in the annotation processing panel, the success of the annotation and pairing with averaged MS spectra could be controlled in a summary table generated automatically.

2.7.6. UI of Output Tab

The output tab was constituted of six panels (Supplementary Figure S5). The first one serves to download the report generated during the analysis. The next four panels allow users to download, respectively, the parameters applied, the intensity matrix, the figures and the HDF5 files, including parameters, averaged spectra, intensity matrix and the annotated table. All these outputs can be downloaded as a single archive in the last panel named “All files”.

3. Assessment of MSProfileR Tool

Although MSProfileR can be applied beyond the entomological domain, its performance was assessed here with two datasets from arthropods. The two datasets of MS spectra were very different from each other, both in terms of sample diversity, but also in terms of the number of spectra tested. The first dataset was intentionally heterogeneous and included several distinct arthropod families. The second consisted in a single mosquito genera but included 13 species with a larger dataset.

3.1. Use Case N°1: Arthropod Families

This first dataset is composed of several arthropod species all laboratory reared. It comprises three arthropod families, mosquitoes (Culicidae), ticks (Ixodidae) and fleas (Pulicidae), including six species at two developmental stages for some of them. According to family and developmental stages, different body parts were submitted to MS analysis. Several recent papers established standardized operational protocols for sample preparation of these arthropod families [7,8,11]. Details about the arthropod species used for MS spectra analysis are summarized in Table 2. The breeding conditions and sample preparations are also indicated in this table, with respective references. This first dataset consisted in a total of 192 MS spectra coming from 48 samples, analysed in quadruplicate. These spectra come from our home-made reference spectra DB [30] and are then considered as high quality (Supplementary file S1).

The 192 MS spectra were uploaded in MSProfileR. The number of files uploaded and the name of samples can be controlled on the “data loading” tab (Supplementary Figure S1). As it is the first dataset of spectra analysed with this tool, no file containing MSProfileR parameters set in a previous analysis was available. Then, all the parameters were selected in the present analysis.

In the preprocessing tab (Supplementary Figure S2A), MS spectra were trimmed in the range of 2-20 kDa and were submitted to conformity tests (Supplementary Figure S2B and S2C). All spectra passed these tests. The parameters selected in the spectra cleaning panel were: square root method for transforming intensity, Savitzky-Golay method with a half-window size of 10 for smoothing intensity, SNIP method with 100 iterations for removing baseline and TIC method for normalizing the intensity of MS spectra. The result of the cleaning steps can be visualized on a plot for each spectrum (Supplementary Figure S2D). The Q estimator and the RC method with a threshold of 3 were computed for detecting of outliers in the quality control panel (Supplementary Figure S2E). In the present dataset, an upper atypical score (A score) limit of 0.59 revealed six spectra among 192 as outliers (3.1%). Interestingly, unticking the “Include spectra below the lower threshold” button did not exclude additional spectra. This underlines that none of the spectra with a low A score were considered outlier. All spectra classified as outliers are automatically removed from the next steps of the analysis. However, it remains possible to visualize all outliers and keep them (Supplementary Figure S2F). Using MSProfileR, all outlier spectra were compared with their respective replicates (Supplementary Figure S6). These spectra with higher A score (i.e., upper the limit) were clearly confirmed to be of lower quality, particularly because they presented a higher background noise. These spectra came from 5 distinct samples, two replicates from one Ae. albopictus at larval stage, two spectra from legs of two Ae. aegypti specimens and two spectra from cephalothoraxes of two Ct. felis specimens. As no sample had four replicates considered outliers, removing the six outlier spectra did not suppress any sample, which was confirmed in the averaging panel (Supplementary Figure S2G). The averaging of selected replicates was done with the mean methods revealing that 48 samples were available for further analysis.

In the processing tab (Supplementary Figure S3A), the MAD method with a half-window of 20 was applied for peak detection. Peak detection is directly linked to the signal-to-noise (SNR) value selected. In this dataset, when the SNR increases from 2 to 7, the mean number of peaks detected per spectrum decreases from about 144 to 39, respectively (Supplementary Figure S3B). The choice of the optimal SNR value could be decisive for the next steps of the analysis. To determine the most appropriate SNR value, a comparison of the peak list detected per arthropod family and per body part was carried out based on the SNR value from two to five (Supplementary Figure S7). It was noticed that at SNR equal to two, numerous peaks were detected between 2 to 8 kDa, among which several are unspecific and correspond to background noise. At SNR equal to three, the majority of these peaks of very low intensity were deleted. Moving to SNR equal to four, regardless of the sample types, some peaks which do not correspond to background noise were no longer detected. In such condition, some information about the protein profile of the sample was lost. To avoid this phenomenon, a SNR value of three was selected for this dataset. At this SNR value, the number of peaks detected per spectrum was 103.3 ± 11.5 (mean ± standard deviation (SD)).

After peak detection, the alignment of spectra was done (Supplementary Figure S3C). The strict method with a minimum frequency of 50% and a tolerance of 0.002 were set for selecting reference peaks. A total of 33 peaks met these criteria and were used to align spectra from the dataset by applying lowess method for spectra warping. For the binning step, the strict method with a tolerance of 0.002 was set and the parameter for peak filtering was a prevalence of at least 20% for a m/z value to be conserved (Supplementary Figure S3D and S3E). Among the 516 peak positions obtained after the binning step, the filtering reduced this number to 202. Finally, these 202 peaks were used to classify each sample based on their averaged MS profiles using a hierarchical clustering (Supplementary Figure S3F). The clustering confirmed that all averaged spectra from the same sample type were grouped per species and body parts. Interestingly, the first criterion of classification is the body part followed by the species (Supplementary Figure S8).

In the annotation tab (Supplementary Figure S4A), the list of the sample name of the 48 averaged spectra can be downloaded and was completed with several metadata from each sample prior to be uploaded (Supplementary Figure S4B and S4C). The annotation processing table allows to verify the correspondence of the averaged spectra with the metadata (Supplementary Figure S4D). Here, as the number of averaged spectra and annotated samples are identical (n=48) and as none of them was not paired, then all averaged spectra from the dataset N°1 were annotated (Supplementary fileS1).

In the output tab (Supplementary Figure S5A and S5(a)), several kinds of files could be downloaded. The report file detailed modules, methods and parameters applied, plus information about the dataset N°1 and the results (tables, plots, boxplot, gelviews and clustered heatmap) of the preprocessing and processing parts (Supplementary Figure S5B and S5(b)). The parameters could be uploaded as a JSON file and used for analysing a new dataset (Supplementary Figure S5C and S5(c)). This also ensures the traceability and the reproducibility of the analysis. The plots generated by the application during dataset N°1 analysis can also be exported (Supplementary Figure S5E and S5(e)), as well as the intensity matrix (Supplementary Figure S5D and S5(d)) and the HDF5 files (Supplementary Figure S5F and S5(f)). All the outputs from this dataset are downloadable in one zipped file which is available in the Supplementary file S2.

3.2. Use Case N°2: Culex Genus

This second dataset included MS spectra from Neotropical Culex mosquitoes (Diptera: Culicidae) coming from a work recently published [31]. The methods and protocols used for mosquito collection, for their morphological and/or molecular identification, for sample preparation and MS submission were detailed in this previous study [31]. Briefly, this dataset was composed of MS spectra from legs and thoraxes of 13 distinct Culex species collected in the field, in French Guiana. The number of specimens per species varied according to the availability of the sample, ranging from 1 to 34 (Supplementary file S3).

The dataset N°2 has interesting characteristics for the evaluation of MSProfileR. First, it includes cryptic species, morphologically undistinguishable. Secondly, specimens were collected in the field and then higher inter-sample variations could occurred among specimens from the same species compared to laboratory reared specimens. Thirdly, it encompasses a large number of samples, 169 mosquito specimens submitted to MALDI-TOF MS on two distinct body parts (legs and thoraxes) and loaded in quadruplicates, corresponding to a total of 1352 spectra (169 specimens × two body parts × four replicates). Finally, in the previous work using the same dataset [31], some spectra were excluded by the authors based on their visual inspection which were considered as low quality. The MSProfileR tool offers now the opportunity to detect and visualize all spectra considered as outliers and can then automatically exclude them from the dataset. The comparison of the spectra list excluded by the experimenters on the same dataset in the previous work will allow to assess the relevance of MSProfileR tool.

For analyzing the dataset N°2, the parameters applied for the dataset N°1 were uploaded. Uniquely the modified parameters compared to the dataset N°1 were specified. In the preprocessing part, all the 1352 spectra passed the conformity tests. However, the quality control steps revealed that 65 spectra (4.8%) exceeded the upper limit of the atypical threshold (A score = 0.72) and were considered as outliers (Figure 3D). Interestingly, all spectra classified as outliers had leg origin (Figure 3A, Supplementary file S4).

Several studies reported a higher diversity of spectra between body parts than between species [9,32]. Moreover, spectra from thoraxes generally had more numerous and higher intensity peaks than the paired leg spectra [11,31]. These two factors could likely explain why numerous spectra originating from legs were classified as outliers. To avoid this bias of exclusion, dataset N°2 was analyzed by splitting the list of spectra per body part, legs or thoraxes.

By applying the same parameters, “A” score thresholds of 0.93 and 0.59 were obtained for leg and thorax spectra, respectively (See report in Supplementary file S5 and Supplementary file S6). Nine and 20 spectra from leg and thorax, respectively, were classified as outliers. The inspection of these spectra confirmed the low quality of these outlier profiles compared to typical spectra. However, some spectra with a low-quality profile now scored below the threshold. (Figure 3E). Then, the threshold parameter of the quality control step was adjusted at 1.5 in order to exclude these spectra of low quality. “A” score thresholds of 0.73 and 0.52 were obtained for leg (Figure 3B) and thorax (Figure 3C) spectra, respectively. Details of the spectra listed as outliers are available in the reports of leg (Supplementary file S7) and thorax datasets (Supplementary file S8). Although 62 (9.2%) and 70 (10.4%) spectra from leg and thorax samples were classified as outliers (Figure 3D), the number of samples excluded remains modest, less than 5%. Effectively, the four spectra replicates from leg and thorax of eight and six specimens were excluded from the analysis. It is interesting to note that the eight specimens for which leg spectra were classified as outliers, encompassed those from the two Culex idottus which were excluded by the authors due to the low intensity and inter-sample heterogeneity of MS profiles in the previous work [31]. The inspection of the other spectra for which the four replicates were classified as outliers from legs or thoraxes confirmed the lower quality of these protein profiles. These samples were then excluded from the analysis. A total of 161 and 163 averaged spectra from legs or thoraxes, respectively, passed the preprocessing steps.

To check whether MSProfileR tool succeeded in classifying spectra according to species for each body part, averaged spectra were submitted to the peak detection steps by applying the parameters used in the dataset N°1. The filtering step revealed that 204 and 194 peaks for legs (Supplementary file S7) and thoraxes (Supplementary file S8), respectively, were retained for averaged spectra classification. Hierarchical clustering revealed that thorax averaged spectra were grouped per species. Solely two Cx. usquatus (#TH_24 and #TH_82) were not clustered with the other spectra from the same species; and the thorax averaged spectra from the unique Cx. adamesi (#TH_2) was classified inside the Cx. dunni group (Supplementary figure S10A). The clustering of leg averaged spectra appeared less efficient, with several samples intertwining with other species (Supplementary figure S10B). The lower classification of the leg averaged spectra compared to thorax was attributed, in part, to their lower spectra intensity. Effectively, as shown in the present study and in concordance with previous works, leg spectra were generally less intense than those of thoraxes from the same specimens [11,31].

For thoraxes, among the seven samples, from the previous analysis [31], which did not reach the threshold value (LSV>1.8) for relevant identification, for five of them, two or more of their spectra replicates were considered as outliers, leading to the exclusion of three samples, for which all replicates overtaken the “A” score threshold. The two other thorax samples (#TH_233 and #TH_234) from Cx. dunni conserved by the quality control step, obtained nearly relevant identification scores, 1.78 and 1.74, respectively, in the previous study [31]. The clustering of these last two thorax averaged spectra (#TH_233 and #TH_234) with the other samples from the same species support the conservation of their respective spectra by quality control steps.

Interestingly, the previous work [31] reported that the selection of the top ten of mass peak list per Culex species and per body part appeared sufficient to discriminate these 13 Culex species with a correct classification upper than 90%. Here, during the filtering step, all peaks with a frequency lower than 20 % (i.e., 0.2) across the dataset were excluded from the analysis. In the dataset N°2, for five species, one to five specimens were available which represents less than 1 to 4% of the number of averaged spectra of the dataset N°2 (161 for legs and 163 for thoraxes). It is then possible that some peaks specific to these species were not included in the intensity matrix due to their low representation (i.e., too few specimens from the same species to reach filtering threshold), which could explain the imperfect clustering of the samples per Culex species, notably for legs.

The decrease of the peak filtering from 20% to 0.5% (i.e., 0.005) led to the inclusion of 878 peaks in the intensity matrix without improving drastically the classification of leg averaged spectra. Among the peaks added in the intensity matrix, about 81.7% (n=717), some of them are heterogeneous between averaged spectra from the same species which should perturb clustering. MSProfileR then appears well adapted for the detection of atypical spectra, but also for their classification at condition that a subgroup is not too under-representative, which could alter the classification. Moreover, spectra of high intensity remain essential for relevant classification.

The duration of computational analyse is another parameter to take into account. Generally, the quickness of the analysis is directly linked to the power of the computer used. In the present work, the analysis was done on a laptop with classical configuration (Processor: Intel ^(R) Core^(TM) i7-10510U CPU @ 1.80 GHz, Hard disk: RAM: 16 GiB, Graphics card: 00:02.0 VGA compatible controller: Intel Corporation UHD Graphics (rev 02), Operating System: Ubuntu 20.04.). The complete pipeline process for analysing the 1352 spectra of dataset N°2 took less than four minutes, by applying default parameters. This processing speed allows to the user to compare/adjust methods and parameters throughout the workflow without a loss of time.