Skip to contents

Installation

Run this code below to install MSstatsBioNet from bioconductor

if (!require("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

BiocManager::install("MSstatsBioNet")

Purpose of MSstatsBioNet

The MSstatsBioNet package is a member of the MSstats family of packages. It contains a set of functions for interpretation of mass spectrometry (MS) statistical analysis results in the context of protein-protein interaction networks. The package is designed to be used in conjunction with the MSstats package.

Dataset

We will be taking a subset of the dataset found in this paper.

input = data.table::fread(system.file(
    "extdata/msstats.csv",
    package = "MSstatsBioNet"
))

MSstats Convert from Upstream Dataset

library(MSstatsConvert)
msstats_imported = FragPipetoMSstatsFormat(input, use_log_file = FALSE)
#> INFO  [2025-02-07 00:26:05] ** Raw data from FragPipe imported successfully.
#> INFO  [2025-02-07 00:26:05] ** Using annotation extracted from quantification data.
#> INFO  [2025-02-07 00:26:05] ** Run labels were standardized to remove symbols such as '.' or '%'.
#> INFO  [2025-02-07 00:26:05] ** The following options are used:
#>   - Features will be defined by the columns: PeptideSequence, PrecursorCharge, FragmentIon, ProductCharge
#>   - Shared peptides will be removed.
#>   - Proteins with single feature will not be removed.
#>   - Features with less than 3 measurements across runs will be removed.
#> INFO  [2025-02-07 00:26:05] ** Features with all missing measurements across runs are removed.
#> INFO  [2025-02-07 00:26:05] ** Shared peptides are removed.
#> INFO  [2025-02-07 00:26:05] ** Multiple measurements in a feature and a run are summarized by summaryforMultipleRows: max
#> INFO  [2025-02-07 00:26:05] ** Features with one or two measurements across runs are removed.
#> INFO  [2025-02-07 00:26:05] ** Run annotation merged with quantification data.
#> INFO  [2025-02-07 00:26:05] ** Features with one or two measurements across runs are removed.
#> INFO  [2025-02-07 00:26:05] ** Fractionation handled.
#> INFO  [2025-02-07 00:26:05] ** Updated quantification data to make balanced design. Missing values are marked by NA
#> INFO  [2025-02-07 00:26:05] ** Finished preprocessing. The dataset is ready to be processed by the dataProcess function.
head(msstats_imported)
#>   ProteinName PeptideSequence PrecursorCharge FragmentIon ProductCharge
#> 1      P05023   AVAGDASESALLK               2         b12             1
#> 2      P05023   AVAGDASESALLK               2         b12             1
#> 3      P05023   AVAGDASESALLK               2         b12             1
#> 4      P05023   AVAGDASESALLK               2         b12             1
#> 5      P05023   AVAGDASESALLK               2         b12             1
#> 6      P05023   AVAGDASESALLK               2         b12             1
#>   IsotopeLabelType Condition BioReplicate
#> 1                L       NAT            1
#> 2                L         T            1
#> 3                L       NAT            2
#> 4                L         T            2
#> 5                L       NAT            3
#> 6                L         T            3
#>                                              Run Fraction Intensity
#> 1 CPTAC_CCRCC_W_JHU_20190112_LUMOS_C3L-00418_NAT        1  69915.91
#> 2   CPTAC_CCRCC_W_JHU_20190112_LUMOS_C3L-00418_T        1 101726.53
#> 3 CPTAC_CCRCC_W_JHU_20190112_LUMOS_C3L-00606_NAT        1  36215.18
#> 4   CPTAC_CCRCC_W_JHU_20190112_LUMOS_C3L-00606_T        1  24927.36
#> 5 CPTAC_CCRCC_W_JHU_20190112_LUMOS_C3L-01882_NAT        1        NA
#> 6   CPTAC_CCRCC_W_JHU_20190112_LUMOS_C3L-01882_T        1        NA

We will first convert the input data to a format that can be processed by MSstats. In this example, the MetamorpheusToMSstatsFormat function is used to convert the input data from Metamorpheus to MSstats format. The MSstats format will contain the necessary information for downstream analysis, such as peptide information, abundance values, run ID, and experimental annotation information.

MSstats Process and GroupComparison

library(MSstats)
#> 
#> Attaching package: 'MSstats'
#> The following objects are masked from 'package:MSstatsConvert':
#> 
#>     DIANNtoMSstatsFormat, DIAUmpiretoMSstatsFormat,
#>     FragPipetoMSstatsFormat, MaxQtoMSstatsFormat,
#>     OpenMStoMSstatsFormat, OpenSWATHtoMSstatsFormat, PDtoMSstatsFormat,
#>     ProgenesistoMSstatsFormat, SkylinetoMSstatsFormat,
#>     SpectronauttoMSstatsFormat
#> The following object is masked from 'package:grDevices':
#> 
#>     savePlot
QuantData <- dataProcess(msstats_imported, use_log_file = FALSE)
#> INFO  [2025-02-07 00:26:07] ** Log2 intensities under cutoff = 9.2885  were considered as censored missing values.
#> INFO  [2025-02-07 00:26:07] ** Log2 intensities = NA were considered as censored missing values.
#> INFO  [2025-02-07 00:26:07] ** Use all features that the dataset originally has.
#> INFO  [2025-02-07 00:26:07] 
#>  # proteins: 10
#>  # peptides per protein: 1-8
#>  # features per peptide: 7-12
#> INFO  [2025-02-07 00:26:07] 
#>                     NAT T
#>              # runs   5 5
#>     # bioreplicates   5 5
#>  # tech. replicates   1 1
#> INFO  [2025-02-07 00:26:07] Some features are completely missing in at least one condition:  
#>  TALNHC(UniMod:4)NLC(UniMod:4)R_2_b5_1,
#>  TALNHC(UniMod:4)NLC(UniMod:4)R_2_y2_1,
#>  TALNHC(UniMod:4)NLC(UniMod:4)R_2_y3_1,
#>  TALNHC(UniMod:4)NLC(UniMod:4)R_2_y5_1,
#>  TALNHC(UniMod:4)NLC(UniMod:4)R_2_y6_1 ...
#> INFO  [2025-02-07 00:26:07]  == Start the summarization per subplot...
#>   |                                                                              |                                                                      |   0%  |                                                                              |=======                                                               |  10%  |                                                                              |==============                                                        |  20%  |                                                                              |=====================                                                 |  30%  |                                                                              |============================                                          |  40%  |                                                                              |===================================                                   |  50%  |                                                                              |==========================================                            |  60%  |                                                                              |=================================================                     |  70%  |                                                                              |========================================================              |  80%  |                                                                              |===============================================================       |  90%  |                                                                              |======================================================================| 100%
#> INFO  [2025-02-07 00:26:07]  == Summarization is done.
model <- groupComparison(
    contrast.matrix = "pairwise",
    data = QuantData,
    use_log_file = FALSE
)
#> INFO  [2025-02-07 00:26:07]  == Start to test and get inference in whole plot ...
#>   |                                                                              |                                                                      |   0%  |                                                                              |=======                                                               |  10%  |                                                                              |==============                                                        |  20%  |                                                                              |=====================                                                 |  30%  |                                                                              |============================                                          |  40%  |                                                                              |===================================                                   |  50%  |                                                                              |==========================================                            |  60%  |                                                                              |=================================================                     |  70%  |                                                                              |========================================================              |  80%  |                                                                              |===============================================================       |  90%  |                                                                              |======================================================================| 100%
#> INFO  [2025-02-07 00:26:08]  == Comparisons for all proteins are done.
head(model$ComparisonResult)
#>   Protein    Label     log2FC        SE     Tvalue DF       pvalue  adj.pvalue
#> 1  O00217 NAT vs T  2.0285031 0.4364177   4.648077  4 0.0096753522 0.013821932
#> 2  O00330 NAT vs T  1.3000941 0.1320659   9.844285  3 0.0022285171 0.004457034
#> 3  O60313 NAT vs T  0.9299641 0.2387992   3.894336  4 0.0176257617 0.019584180
#> 4  O60879 NAT vs T -1.9484511 0.1695323 -11.493098  4 0.0003271868 0.001635934
#> 5  O75306 NAT vs T  2.4745040 0.3530857   7.008225  4 0.0021825029 0.004457034
#> 6  P05023 NAT vs T  1.8391155 0.2122058   8.666660  4 0.0009753219 0.003251073
#>   issue MissingPercentage ImputationPercentage
#> 1    NA        0.14166667           0.14166667
#> 2    NA        0.20833333           0.10833333
#> 3    NA        0.29756098           0.29756098
#> 4    NA        0.09583333           0.09583333
#> 5    NA        0.16111111           0.16111111
#> 6    NA        0.11702128           0.11702128

Next, we will preprocess the data using the dataProcess function and perform a statistical analysis using the groupComparison function. The output of groupComparison is a table containing a list of differentially abundant proteins with their associated p-values and log fold changes.

MSstatsBioNet Analysis

ID Conversion

First, we need to convert the group comparison results to a format that can be processed by INDRA. The getSubnetworkFromIndra function requires a table containing HGNC IDs. We can use the annotateProteinInfoFromIndra function to obtain these mappings.

In the below example, we convert uniprot IDs to their corresponding Hgnc IDs. We can also extract other information, such as hgnc gene name and protein function.

library(MSstatsBioNet)
annotated_df = annotateProteinInfoFromIndra(model$ComparisonResult, "Uniprot")
head(annotated_df)
#>   Protein    Label     log2FC        SE     Tvalue DF       pvalue  adj.pvalue
#> 1  O00217 NAT vs T  2.0285031 0.4364177   4.648077  4 0.0096753522 0.013821932
#> 2  O00330 NAT vs T  1.3000941 0.1320659   9.844285  3 0.0022285171 0.004457034
#> 3  O60313 NAT vs T  0.9299641 0.2387992   3.894336  4 0.0176257617 0.019584180
#> 4  O60879 NAT vs T -1.9484511 0.1695323 -11.493098  4 0.0003271868 0.001635934
#> 5  O75306 NAT vs T  2.4745040 0.3530857   7.008225  4 0.0021825029 0.004457034
#> 6  P05023 NAT vs T  1.8391155 0.2122058   8.666660  4 0.0009753219 0.003251073
#>   issue MissingPercentage ImputationPercentage UniprotId HgncId HgncName
#> 1    NA        0.14166667           0.14166667    O00217   7715   NDUFS8
#> 2    NA        0.20833333           0.10833333    O00330  21350     PDHX
#> 3    NA        0.29756098           0.29756098    O60313   8140     OPA1
#> 4    NA        0.09583333           0.09583333    O60879   2877   DIAPH2
#> 5    NA        0.16111111           0.16111111    O75306   7708   NDUFS2
#> 6    NA        0.11702128           0.11702128    P05023    799   ATP1A1
#>   IsTranscriptionFactor IsKinase IsPhosphatase
#> 1                 FALSE    FALSE         FALSE
#> 2                 FALSE    FALSE         FALSE
#> 3                 FALSE    FALSE         FALSE
#> 4                 FALSE    FALSE         FALSE
#> 5                 FALSE    FALSE         FALSE
#> 6                 FALSE    FALSE         FALSE

Subnetwork Query

The package provides a function getSubnetworkFromIndra that retrieves a subnetwork of proteins from the INDRA database based on differential abundance analysis results.

subnetwork <- getSubnetworkFromIndra(annotated_df, pvalueCutoff = 0.05)
#> Warning in getSubnetworkFromIndra(annotated_df, pvalueCutoff = 0.05): NOTICE: This function includes third-party software components
#>         that are licensed under the BSD 2-Clause License. Please ensure to
#>         include the third-party licensing agreements if redistributing this
#>         package or utilizing the results based on this package.
#>         See the LICENSE file for more details.
head(subnetwork$nodes)
#>       id     logFC     pvalue hgncName
#> 3 O60313 0.9299641 0.01958418     OPA1
#> 7 P05067 0.7360012 0.02030666      APP
#> 8 P05090 0.5683951 0.01371505     APOD
head(subnetwork$edges)
#>   source target    interaction evidenceCount paperCount
#> 1 P05067 O60313 DecreaseAmount             5          2
#> 2 P05067 P05090 IncreaseAmount             1          1
#>                                                                                evidenceLink
#> 1 https://db.indra.bio/statements/from_agents?subject=620@HGNC&object=8140@HGNC&format=html
#> 2  https://db.indra.bio/statements/from_agents?subject=620@HGNC&object=612@HGNC&format=html

This package is distributed under the Artistic-2.0 license. However, its dependencies may have different licenses. In this example, getSubnetworkFromIndra depends on INDRA, which is distributed under the BSD 2-Clause license. Furthermore, INDRA’s knowledge sources may have different licenses for commercial applications. Please refer to the INDRA README for more information on its knowledge sources and their associated licenses.

Visualize Networks

The function visualizeNetworks then takes the output of getSubnetworkFromIndra and visualizes the subnetwork. The function requires Cytoscape desktop to be open for the visualization to work.

visualizeNetworks(subnetwork$nodes, subnetwork$edges)
#> Warning in visualizeNetworks(subnetwork$nodes, subnetwork$edges): Visualization
#> is not available in non-interactive mode.

In the network diagram displayed in Cytoscape, you should see two arrows connecting two nodes, P16050 and P84243. These arrows represent the interactions between these two proteins, notably activation and phosphorylation.

Session info

sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] MSstatsBioNet_0.99.9  MSstats_4.14.1        MSstatsConvert_1.16.1
#> [4] BiocStyle_2.34.0     
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.1      viridisLite_0.4.2     IRdisplay_1.1        
#>  [4] dplyr_1.1.4           bitops_1.0-9          RCurl_1.98-1.16      
#>  [7] fastmap_1.2.0         lazyeval_0.2.2        base64url_1.4        
#> [10] XML_3.99-0.18         digest_0.6.37         lifecycle_1.0.4      
#> [13] survival_3.7-0        statmod_1.5.0         r2r_0.1.2            
#> [16] magrittr_2.0.3        compiler_4.4.2        rlang_1.1.5          
#> [19] sass_0.4.9            tools_4.4.2           yaml_2.3.10          
#> [22] data.table_1.16.4     knitr_1.49            htmlwidgets_1.6.4    
#> [25] curl_6.2.0            marray_1.84.0         RColorBrewer_1.1-3   
#> [28] repr_1.1.7            KernSmooth_2.23-24    pbdZMQ_0.3-13        
#> [31] purrr_1.0.4           BiocGenerics_0.52.0   desc_1.4.3           
#> [34] stats4_4.4.2          grid_4.4.2            preprocessCore_1.68.0
#> [37] caTools_1.18.3        colorspace_2.1-1      log4r_0.4.4          
#> [40] ggplot2_3.5.1         scales_1.3.0          gtools_3.9.5         
#> [43] MASS_7.3-61           cli_3.6.3             crayon_1.5.3         
#> [46] rmarkdown_2.29        ragg_1.3.3            reformulas_0.4.0     
#> [49] generics_0.1.3        httr_1.4.7            minqa_1.2.8          
#> [52] cachem_1.1.0          splines_4.4.2         parallel_4.4.2       
#> [55] BiocManager_1.30.25   base64enc_0.1-3       vctrs_0.6.5          
#> [58] boot_1.3-31           Matrix_1.7-1          jsonlite_1.8.9       
#> [61] bookdown_0.42         ggrepel_0.9.6         systemfonts_1.2.1    
#> [64] limma_3.62.2          plotly_4.10.4         jquerylib_0.1.4      
#> [67] tidyr_1.3.1           glue_1.8.0            nloptr_2.1.1         
#> [70] pkgdown_2.1.1         RJSONIO_1.3-1.9       stringi_1.8.4        
#> [73] gtable_0.3.6          lme4_1.1-36           munsell_0.5.1        
#> [76] tibble_3.2.1          pillar_1.10.1         htmltools_0.5.8.1    
#> [79] gplots_3.2.0          RCy3_2.26.0           graph_1.84.1         
#> [82] IRkernel_1.3.2        R6_2.5.1              textshaping_1.0.0    
#> [85] Rdpack_2.6.2          evaluate_1.0.3        lattice_0.22-6       
#> [88] rbibutils_2.3         backports_1.5.0       bslib_0.9.0          
#> [91] Rcpp_1.0.14           uuid_1.2-1            nlme_3.1-166         
#> [94] checkmate_2.3.2       xfun_0.50             fs_1.6.5             
#> [97] pkgconfig_2.0.3