Run this code below to install MSstatsBioNet from bioconductor
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.
We will be taking a subset of the dataset found in this paper.
library(MSstatsConvert)
msstats_imported = FragPipetoMSstatsFormat(input, use_log_file = FALSE)
#> INFO [2026-10-10 01:09:00] ** Raw data from FragPipe imported successfully.
#> INFO [2026-10-10 01:09:00] ** Using annotation extracted from quantification data.
#> INFO [2026-10-10 01:09:00] ** Run labels were standardized to remove symbols such as '.' or '%'.
#> INFO [2026-10-10 01:09:00] ** 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 [2026-10-10 01:09:00] ** Features with all missing measurements across runs are removed.
#> INFO [2026-10-10 01:09:00] ** Shared peptides are removed.
#> INFO [2026-10-10 01:09:00] ** Multiple measurements in a feature and a run are summarized by summaryforMultipleRows: max
#> INFO [2026-10-10 01:09:00] ** Features with one or two measurements across runs are removed.
#> INFO [2026-10-10 01:09:00] ** Run annotation merged with quantification data.
#> INFO [2026-10-10 01:09:00] ** Features with one or two measurements across runs are removed.
#> INFO [2026-10-10 01:09:00] ** Fractionation handled.
#> INFO [2026-10-10 01:09:00] ** Updated quantification data to make balanced design. Missing values are marked by NA
#> INFO [2026-10-10 01:09:00] ** 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 NAWe 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.
library(MSstats)
QuantData <- dataProcess(msstats_imported, use_log_file = FALSE)
#> INFO [2026-10-10 01:09:00] ** Log2 intensities under cutoff = 9.2885 were considered as censored missing values.
#> INFO [2026-10-10 01:09:00] ** Log2 intensities = NA were considered as censored missing values.
#> INFO [2026-10-10 01:09:00] ** Use top100 features that have highest average of log2(intensity) across runs.
#> INFO [2026-10-10 01:09:00]
#> # proteins: 10
#> # peptides per protein: 1-8
#> # features per peptide: 7-12
#> INFO [2026-10-10 01:09:00]
#> NAT T
#> # runs 5 5
#> # bioreplicates 5 5
#> # tech. replicates 1 1
#> INFO [2026-10-10 01:09:00] Some features are completely missing in at least one condition:
#> ELEAEIQQLR_2_b5_1,
#> ELEAEIQQLR_2_b8_1,
#> NLEAVETLGSTSTIC(UniMod:4)SDK_3_b13_2,
#> NLEAVETLGSTSTIC(UniMod:4)SDK_3_b3_1,
#> NLEAVETLGSTSTIC(UniMod:4)SDK_3_b4_1 ...
#> INFO [2026-10-10 01:09:00] == Start the summarization per subplot...
#> | | | 0% | |======= | 10% | |============== | 20% | |===================== | 30% | |============================ | 40% | |=================================== | 50% | |========================================== | 60% | |================================================= | 70% | |======================================================== | 80% | |=============================================================== | 90% | |======================================================================| 100%
#> INFO [2026-10-10 01:09:00] == Summarization is done.
model <- groupComparison(
contrast.matrix = "pairwise",
data = QuantData,
use_log_file = FALSE
)
#> INFO [2026-10-10 01:09:01] == Start to test and get inference in whole plot ...
#> | | | 0% | |======= | 10% | |============== | 20% | |===================== | 30% | |============================ | 40% | |=================================== | 50% | |========================================== | 60% | |================================================= | 70% | |======================================================== | 80% | |=============================================================== | 90% | |======================================================================| 100%
#> INFO [2026-10-10 01:09:01] == 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.0022285170 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.11702128Next, 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.
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.0022285170 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 GlobalProtein UniprotId
#> 1 NA 0.14166667 0.14166667 O00217 O00217
#> 2 NA 0.20833333 0.10833333 O00330 O00330
#> 3 NA 0.29756098 0.29756098 O60313 O60313
#> 4 NA 0.09583333 0.09583333 O60879 O60879
#> 5 NA 0.16111111 0.16111111 O75306 O75306
#> 6 NA 0.11702128 0.11702128 P05023 P05023
#> EntityNamespace EntityId EntityName IsTranscriptionFactor IsKinase
#> 1 HGNC 7715 NDUFS8 FALSE FALSE
#> 2 HGNC 21350 PDHX FALSE FALSE
#> 3 HGNC 8140 OPA1 FALSE FALSE
#> 4 HGNC 2877 DIAPH2 FALSE FALSE
#> 5 HGNC 7708 NDUFS2 FALSE FALSE
#> 6 HGNC 799 ATP1A1 FALSE FALSE
#> IsPhosphatase
#> 1 FALSE
#> 2 FALSE
#> 3 FALSE
#> 4 FALSE
#> 5 FALSE
#> 6 FALSEThe 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,
statement_types = c("Complex", "IncreaseAmount", "DecreaseAmount", "Inhibition", "Activation", "Phosphorylation")
)
#> INDRA subnetwork: how are 10 selected proteins connected to each other, with no other nodes added?
#> Warning in getSubnetworkFromIndra(annotated_df, pvalueCutoff = 0.05, statement_types = c("Complex", : 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 entity_type entity_name namespace entity_id measured included_in_query
#> 1 O00217 protein NDUFS8 HGNC 7715 TRUE TRUE
#> 2 O60313 protein OPA1 HGNC 8140 TRUE TRUE
#> 3 O75306 protein NDUFS2 HGNC 7708 TRUE TRUE
#> 4 P05023 protein ATP1A1 HGNC 799 TRUE TRUE
#> 5 P05067 protein APP HGNC 620 TRUE TRUE
#> 6 P05090 protein APOD HGNC 612 TRUE TRUE
#> node_role site has_measured_sites logFC adj.pvalue
#> 1 passed_cutoffs <NA> FALSE 2.0285031 0.013821932
#> 2 passed_cutoffs <NA> FALSE 0.9299641 0.019584180
#> 3 passed_cutoffs <NA> FALSE 2.4745040 0.004457034
#> 4 passed_cutoffs <NA> FALSE 1.8391155 0.003251073
#> 5 passed_cutoffs <NA> FALSE 0.7360012 0.020306662
#> 6 passed_cutoffs <NA> FALSE 0.5683951 0.013715050
head(subnetwork$edges)
#> source target interaction directed site confidence evidence_count
#> 1 O75306 O00217 Complex FALSE <NA> 0.7981153 3
#> 2 P08574 O75306 Complex FALSE <NA> 0.8302067 1
#> 3 P05067 O60313 DecreaseAmount TRUE <NA> 0.4850742 3
#> 4 P05067 P05362 IncreaseAmount TRUE <NA> 0.5028536 3
#> 5 P05067 P05090 IncreaseAmount TRUE <NA> 0.4715316 2
#> 6 P08574 P05023 Complex FALSE <NA> 0.8302067 1
#> evidence_url
#> 1 https://db.indra.bio/statements/from_hash/35453498871697150?format=html
#> 2 https://db.indra.bio/statements/from_hash/6349003830434161?format=html
#> 3 https://db.indra.bio/statements/from_hash/-14390698362022978?format=html
#> 4 https://db.indra.bio/statements/from_hash/29503060009374668?format=html
#> 5 https://db.indra.bio/statements/from_hash/13542679335950199?format=html
#> 6 https://db.indra.bio/statements/from_hash/-27854856016676857?format=html
#> statement_id backend_database query_type evidence_sources paperCount
#> 1 35453498871697150 INDRA subnetwork {"biogrid": 3} 1
#> 2 6349003830434161 INDRA subnetwork {"biogrid": 1} 1
#> 3 -14390698362022978 INDRA subnetwork {"reach": 3} 1
#> 4 29503060009374668 INDRA subnetwork {"reach": 3} 1
#> 5 13542679335950199 INDRA subnetwork {"reach": 2} 1
#> 6 -27854856016676857 INDRA subnetwork {"biogrid": 1} 1This 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.
The function previewNetworkInBrowser then takes the
output of getSubnetworkFromIndra and visualizes the
subnetwork. The function requires an internet browser to view the
subnetwork
In the network diagram displayed using CytoscapeJS, you should see three arrows connecting two nodes, APP and APOD. These arrows represent the interactions between these two proteins, notably Activation, IncreaseAmount (e.g. via transcriptional regulation), and forming a protein complex.
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.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.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] MSstatsConvert_1.23.5 rmarkdown_2.32 shiny_1.14.0
#> [4] MSstatsBioNet_1.5.4 MSstats_4.21.1 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 viridisLite_0.4.3 dplyr_1.2.1
#> [4] farver_2.1.2 S7_0.2.2 bitops_1.1-0
#> [7] fastmap_1.2.0 promises_1.5.0 XML_3.99-0.25
#> [10] digest_0.6.39 mime_0.13 lifecycle_1.0.5
#> [13] survival_3.8-12 statmod_1.5.2 magrittr_2.0.5
#> [16] compiler_4.6.1 r2r_0.1.2 rlang_1.3.0
#> [19] sass_0.4.10 tools_4.6.1 yaml_2.3.12
#> [22] data.table_1.18.6.1 knitr_1.52 stopwords_2.3
#> [25] htmlwidgets_1.6.4 curl_8.0.0 marray_1.91.0
#> [28] xml2_1.6.0 RColorBrewer_1.1-3 KernSmooth_2.23-27
#> [31] purrr_1.2.2 sys_3.4.3 grid_4.6.1
#> [34] preprocessCore_1.75.1 caTools_1.18.4 xtable_1.8-8
#> [37] ggplot2_4.0.3 scales_1.4.0 gtools_3.9.5
#> [40] MASS_7.3-66 cli_3.6.6 crayon_1.5.3
#> [43] reformulas_0.4.4 generics_0.1.4 otel_0.2.0
#> [46] httr_1.4.9 minqa_1.2.8 cachem_1.1.0
#> [49] splines_4.6.1 parallel_4.6.1 BiocManager_1.30.27
#> [52] vctrs_0.7.3 boot_1.3-32 Matrix_1.7-6
#> [55] jsonlite_2.0.0 ggrepel_0.9.8 maketools_1.3.2
#> [58] limma_3.99.0 plotly_4.12.1 lgr_0.5.2
#> [61] tidyr_1.3.2 jquerylib_0.1.4 glue_1.8.1
#> [64] nloptr_2.2.1 stringi_1.8.9 gtable_0.3.6
#> [67] later_1.4.8 lme4_2.0-6 mlapi_0.1.1
#> [70] tibble_3.3.1 pillar_1.11.1 htmltools_0.5.9
#> [73] gplots_3.3.0 float_0.3-4 rsparse_0.5.3
#> [76] R6_2.6.1 Rdpack_2.6.6 evaluate_1.0.5
#> [79] lattice_0.23-1 rentrez_1.2.4 rbibutils_2.4.1
#> [82] backports_1.5.1 RhpcBLASctl_0.23-42 memoise_2.0.1
#> [85] httpuv_1.6.17 bslib_0.12.0 text2vec_0.6.6
#> [88] Rcpp_1.1.2 nlme_3.1-171 checkmate_2.3.4
#> [91] xfun_0.61 buildtools_1.0.0 pkgconfig_2.0.3