Introduction to proteomics data analysis - MaxQuant Data Dependent Acquisition spike-in study
14 April 2026
1 Introduction
Label-Free Quantitative mass spectrometry based workflows for differential expression (DE) analysis of proteins is often challenging due to peptide-specific effects and context-sensitive missingness of peptide intensities.
msqrob2provides peptide-based workflows that can assess for DE directly from peptide intensities and outperform summarisation methods which first aggregate MS1 peptide intensities to protein intensities before DE analysis. However, they are computationally expensive, often hard to understand for the non-specialised end-user, and they do not provide protein summaries, which are important for visualisation or downstream processing.msqrob2therefore also proposes a novel summarisation strategy, which estimates MSqRob’s model parameters in a two-stage procedure circumventing the drawbacks of peptide-based workflows.the summarisation based workflow in
msqrob2maintains MSqRob’s superior performance, while providing useful protein expression summaries for plotting and downstream analysis. Summarising peptide to protein intensities considerably reduces the computational complexity, the memory footprint and the model complexity. Moreover, it renders the analysis framework to become modular, providing users the flexibility to develop workflows tailored towards specific applications.
In this vignette we will demonstrate how to perform msqrob’s summarisation
based workflow starting from a Maxquant search on a subset of the cptac
spike-in study.
Examples on our peptide-based workflows and on the analysis of more complex designs can be found on our companion website msqrob2 book.
Technical details on our methods can be found in (Goeminne, Gevaert, and Clement 2016), (Goeminne et al. 2020), (Sticker et al. 2020) and (Vandenbulcke et al. 2025).
2 Background
This case-study is a subset of the data of the 6th study of the Clinical Proteomic Technology Assessment for Cancer (CPTAC). In this experiment, the authors spiked the Sigma Universal Protein Standard mixture 1 (UPS1) containing 48 different human proteins in a protein background of 60 ng/\(\mu\)L Saccharomyces cerevisiae strain BY4741. Two different spike-in concentrations were used: 6A (0.25 fmol UPS1 proteins/\(\mu\)L) and 6B (0.74 fmol UPS1 proteins/\(\mu\)L) [5]. We limited ourselves to the data of LTQ-Orbitrap W at site 56. The data were searched with MaxQuant version 1.5.2.8, and detailed search settings were described in Goeminne et al. (2016) [1]. Three replicates are available for each concentration.
3 Load packages
We load the msqrob2 package, along with additional packages for
data manipulation and visualisation.
suppressPackageStartupMessages({
library("QFeatures")
library("dplyr")
library("tidyr")
library("ggplot2")
library("msqrob2")
library("stringr")
library("ExploreModelMatrix")
library("kableExtra")
library("ComplexHeatmap")
library("scater")
})4 Data
4.1 Peptide table
We first import the data from the peptide.txt file. This is the file containing
your peptide-level intensities searched by MaxQuant search,
this peptide.txt file can be found by default in the
“path_to_raw_files/combined/txt/” folder from the MaxQuant output,
with “path_to_raw_files” the folder where the raw files were saved.
In this vignette, we use a MaxQuant peptide file which is a subset
of the cptac study.
This data is available in
msqrob2data GitHub repository.
To import the data we use the QFeatures package.
We generate the object peptideFile with the path to the peptide.txt file. This can be a path to the data on the web or a path to data on your local computer.
peptideFile <- "https://github.com/statOmics/msqrob2data/raw/refs/heads/main/dda/cptacAvsB_lab3/peptides.txt"peptides <- data.table::fread(peptideFile)## Registered S3 method overwritten by 'bit64':
## method from
## print.bitstring toolshead(peptides)## Sequence N-term cleavage window C-term cleavage window
## <char> <char> <char>
## 1: AAAAGAGGAGDSGDAVTK EHQHDEQKAAAAGAGG DSGDAVTKIGSEDVKL
## 2: AAAALAGGK QQLSKAAKAAAALAGG AAALAGGKKSKKKWSK
## 3: AAAALAGGKK QQLSKAAKAAAALAGG AALAGGKKSKKKWSKK
## 4: AAADALSDLEIK MPKETPSKAAADALSD ALSDLEIKDSKSNLNK
## 5: AAADALSDLEIKDSK MPKETPSKAAADALSD DLEIKDSKSNLNKELE
## 6: AAAEEFQR RERKEREKAAAEEFQR AAAEEFQRQQELLRQQ
## Amino acid before First amino acid Second amino acid Second last amino acid
## <char> <char> <char> <char>
## 1: K A A T
## 2: K A A G
## 3: K A A K
## 4: K A A I
## 5: K A A S
## 6: K A A Q
## Last amino acid Amino acid after A Count R Count N Count D Count C Count
## <char> <char> <int> <int> <int> <int> <int>
## 1: K I 7 0 0 2 0
## 2: K K 5 0 0 0 0
## 3: K S 5 0 0 0 0
## 4: K D 4 0 0 2 0
## 5: K S 4 0 0 3 0
## 6: R Q 3 1 0 0 0
## Q Count E Count G Count H Count I Count L Count K Count M Count F Count
## <int> <int> <int> <int> <int> <int> <int> <int> <int>
## 1: 0 0 5 0 0 0 1 0 0
## 2: 0 0 2 0 0 1 1 0 0
## 3: 0 0 2 0 0 1 2 0 0
## 4: 0 1 0 0 1 2 1 0 0
## 5: 0 1 0 0 1 2 2 0 0
## 6: 1 2 0 0 0 0 0 0 1
## P Count S Count T Count W Count Y Count V Count U Count Length
## <int> <int> <int> <int> <int> <int> <int> <int>
## 1: 0 1 1 0 0 1 0 18
## 2: 0 0 0 0 0 0 0 9
## 3: 0 0 0 0 0 0 0 10
## 4: 0 1 0 0 0 0 0 12
## 5: 0 2 0 0 0 0 0 15
## 6: 0 0 0 0 0 0 0 8
## Missed cleavages Mass Proteins
## <int> <num> <char>
## 1: 0 1445.6746 sp|P38915|SPT8_YEAST
## 2: 0 728.4181 sp|Q3E792|RS25A_YEAST;sp|P0C0T4|RS25B_YEAST
## 3: 1 856.5131 sp|Q3E792|RS25A_YEAST;sp|P0C0T4|RS25B_YEAST
## 4: 0 1215.6347 sp|P09938|RIR2_YEAST
## 5: 1 1545.7886 sp|P09938|RIR2_YEAST
## 6: 0 920.4352 sp|P53075|SHE10_YEAST
## Leading razor protein Start position End position Unique (Groups)
## <char> <int> <int> <char>
## 1: sp|P38915|SPT8_YEAST 97 114 yes
## 2: sp|Q3E792|RS25A_YEAST 13 21 yes
## 3: sp|Q3E792|RS25A_YEAST 13 22 yes
## 4: sp|P09938|RIR2_YEAST 9 20 yes
## 5: sp|P09938|RIR2_YEAST 9 23 yes
## 6: sp|P53075|SHE10_YEAST 538 545 yes
## Unique (Proteins) Charges PEP Score Identification type 6A_7
## <char> <char> <num> <num> <char>
## 1: yes 2 1.1800e-05 82.942 By MS/MS
## 2: no 2 7.4600e-06 134.810 By MS/MS
## 3: no 2 3.3100e-09 143.730 By MS/MS
## 4: yes 2 9.1600e-23 182.230 By MS/MS
## 5: yes 3 1.5319e-04 73.927 By MS/MS
## 6: yes 2 1.7588e-02 79.474 By matching
## Identification type 6A_8 Identification type 6A_9 Identification type 6B_7
## <char> <char> <char>
## 1: By MS/MS By MS/MS By matching
## 2: By MS/MS By MS/MS By MS/MS
## 3: By MS/MS By MS/MS By MS/MS
## 4: By matching By MS/MS By MS/MS
## 5: By MS/MS By MS/MS By MS/MS
## 6: By matching By matching By matching
## Identification type 6B_8 Identification type 6B_9 Experiment 6A_7
## <char> <char> <int>
## 1: By MS/MS By MS/MS 1
## 2: By MS/MS By MS/MS 2
## 3: By MS/MS By MS/MS 1
## 4: By matching By matching 1
## 5: By MS/MS By MS/MS 1
## 6: By matching By matching NA
## Experiment 6A_8 Experiment 6A_9 Experiment 6B_7 Experiment 6B_8
## <int> <int> <int> <int>
## 1: 1 1 NA 1
## 2: 1 1 2 1
## 3: 1 1 1 1
## 4: 1 1 1 1
## 5: 1 1 1 1
## 6: NA 1 NA NA
## Experiment 6B_9 Intensity Intensity 6A_7 Intensity 6A_8 Intensity 6A_9
## <int> <i64> <num> <num> <num>
## 1: 1 1190800 0 0 66760
## 2: 2 280990000 2441300 1220000 1337600
## 3: 1 33360000 1029200 668040 638990
## 4: 1 54622000 515460 670780 712140
## 5: 1 18910000 331130 420900 365560
## 6: 1 1158600 0 0 51558
## Intensity 6B_7 Intensity 6B_8 Intensity 6B_9 Reverse Potential contaminant
## <num> <num> <num> <char> <char>
## 1: 0 37436 0
## 2: 2850900 935580 1606200
## 3: 777030 641270 562840
## 4: 426580 620510 737780
## 5: 329250 380820 414490
## 6: 0 0 76970
## id Protein group IDs Mod. peptide IDs
## <int> <char> <char>
## 1: 0 859 0
## 2: 1 230 1
## 3: 2 230 2
## 4: 3 229 3
## 5: 4 229 4
## 6: 5 1143 5
## Evidence IDs
## <char>
## 1: 0;1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20;21;22;23
## 2: 24;25;26;27;28;29;30;31;32;33;34;35;36;37;38;39;40;41;42;43;44;45;46;47;48;49;50;51;52;53;54;55;56;57;58;59;60;61;62;63;64;65;66;67;68;69;70;71;72;73
## 3: 74;75;76;77;78;79;80;81;82;83;84;85;86;87;88;89;90;91;92;93;94;95;96;97;98
## 4: 99;100;101;102;103;104;105;106;107;108;109;110;111;112;113;114;115;116;117;118;119;120;121;122;123;124;125;126;127;128;129;130;131;132;133;134;135;136;137;138;139;140;141;142;143
## 5: 144;145;146;147;148;149;150;151;152;153;154;155;156;157;158;159;160;161;162;163;164;165;166;167;168;169;170;171;172;173;174;175;176;177;178;179;180;181;182;183;184
## 6: 185;186;187;188;189;190;191;192;193;194
## MS/MS IDs
## <char>
## 1: 0;1;2;3;4;5;6;7;8;9
## 2: 10;11;12;13;14;15;16;17;18;19;20;21;22;23;24;25;26;27;28;29
## 3: 30;31;32;33;34;35;36;37;38;39;40;41;42;43;44;45;46;47;48;49;50
## 4: 51;52;53;54;55;56;57;58;59;60;61;62;63;64;65;66;67;68;69;70;71;72;73;74;75;76;77;78;79;80;81;82;83;84
## 5: 85;86;87;88;89;90;91;92;93;94;95;96;97;98;99;100;101;102;103;104;105;106;107;108;109;110;111;112;113;114;115;116
## 6: 117;118
## Best MS/MS Oxidation (M) site IDs MS/MS Count
## <int> <char> <int>
## 1: 0 10
## 2: 21 18
## 3: 31 21
## 4: 72 29
## 5: 94 32
## 6: 117 2Each row in the peptide data table contains information about one peptide
(the table above shows the first 6 rows). The columns contains various
descriptors about the peptide, such as its sequence, its charge, the amino acid
composition, etc. Some of these columns (those starting with Intensity)
contain the quantification values for each sample. The table format where the
quantitative values for each sample are contained in a separate column is
depicted as the “wide format”, as opposed to the “long format”
(e.g., the [PSM table]).
quantCols <- grep("Intensity ", names(peptides), value = TRUE)Note, that when you read the data with the default functions in R such as
read.table, the space in Intensity will be replaced with a ..
Unlike for real experiments, we known the ground truth for this dataset: UPS proteins were spiked in at different concentrations and are thus differentially abundant (DA, spiked in). Yeast proteins consist of the background proteome that is the same in all samples and are thus non-DA.
peptides <- peptides |>
mutate(species = grepl(pattern = "UPS",Proteins) |>
as.factor() |>
recode("TRUE"="ups","FALSE" = "yeast"))4.2 Sample annotation
Each row in the annotation table contains information about one sample. The columns contain various descriptors about the sample, such as the name of the sample or the MS run, the treatment (here the spike-in condition), or any other biological or technical information that may impact the data quality or the quantification. Without an annotation table, no analysis can be performed. The sample annotations are generated by the researcher. In this example, the annotations are extracted from the sample names, although reporting a detailed design of experiments in a table is better practice.
- We first generate a variable quantCols, which contain the names of the quantification columns.
- Next we generate variable
sampleId, to pinpoint the different samples in the dataset. Since all quantification columns start with the generic stringIntensity 6we replace the string with an empty string""so we keep the last 3 characters to refer to the sample. - We construct a variable condition, by using the first letter of the sampleId, which refers to the spikein condition.
- We construct a variable rep, by using the third letter of the sampleId.
(annot <- data.frame(quantCols = quantCols) |> #1.
mutate(sampleId = gsub("Intensity 6","", quantCols), #2.
condition = substr(sampleId,1,1), #3.
rep = substr(sampleId, 3, 3)) #4.
)## quantCols sampleId condition rep
## 1 Intensity 6A_7 A_7 A 7
## 2 Intensity 6A_8 A_8 A 8
## 3 Intensity 6A_9 A_9 A 9
## 4 Intensity 6B_7 B_7 B 7
## 5 Intensity 6B_8 B_8 B 8
## 6 Intensity 6B_9 B_9 B 94.3 Convert to QFeatures
The readQFeatures() enables a seamless conversion of tabular data into a
QFeatures object. We provide the peptide table and the sample annotation table.
The function will use the quantCols column in the sample annotation table to
understand which columns in peptides contain the quantitative values, and
automatically link the corresponding sample annotation with the quantitative
values. We also tell the function to use the Sequence column as peptide
identifier, which will be used as rownames. See ?readQFeatures() for more
details.
(qf <- readQFeatures(assayData = peptides,
colData = annot,
name = "peptides",
fnames = "Sequence")
)## Checking arguments.## Loading data as a 'SummarizedExperiment' object.## Formatting sample annotations (colData).## Formatting data as a 'QFeatures' object.## Setting assay rownames.## An instance of class QFeatures (type: bulk) with 1 set:
##
## [1] peptides: SummarizedExperiment with 11466 rows and 6 columns5 Data preprocessing
Since we have a QFeatures object, we can directly make use of QFeatures’
data preprocessing functionality. A major advantage of QFeatures is it provides
the power to build highly modular workflows, where each step is carried out by a
dedicated function with a large choice of available methods and parameters. This
means that users can adapt the workflow to their specific use case and their
specific needs.
5.1 Encoding missing values
The first preprocessing step is to correctly encode the missing values.
It is important that missing values are encoded using NA. For instance,
non-observed values should not be encoded with a zero because true zeros
(the proteomic feature is absent from the sample) cannot be distinguished from
technical zeros (the feature was missed by the instrument or could not be
identified). We therefore replace any zero in the quantitative data with an NA.
qf <- zeroIsNA(qf, "peptides") # convert 0 to NANote that msqrob2 can handle missing data without having to rely on hard-to-verify imputation assumptions, which is our general recommendation. However, msqrob2 does not prevent users from using imputation, which can be performed with impute() from the QFeatures package.
5.2 \(\log_2\) transformation
We perform log2-transformation with logTransform() from the QFeatures
package.
qf <- logTransform(qf, base = 2, i = "peptides", name = "peptides_log")5.3 Peptide-Filtering
Filtering removes low-quality and unreliable peptides that would otherwise introduce noise and artefacts in the data.
5.3.1 Remove failed protein inference
We remove peptides that could not be mapped to a protein. We also remove peptides that cannot be uniquely mapped to a single proteins because its sequence is shared across multiple proteins. This often occurs for homologs or for proteins that share similar functional domains. Shared peptides are often mapped to protein groups, which are the set of proteins that share the given peptides. Protein groups are encoded by appending the constituting protein identifiers, separated by a ";".
qf <- filterFeatures(
qf, ~ Proteins != "" & ## Remove failed protein inference
!grepl(";", Proteins)) ## Remove protein groups## 'Proteins' found in 2 out of 2 assay(s).5.3.2 Remove reverse sequences and contaminants
We now remove the peptides that map to decoy sequences or to
contaminants proteins. In our peptide.txt
file decoys and contaminants are indicated with a “+” character in the
column Reverse and Potential.contaminant, respectively.
Note, that depending on the version of MaxQuant the name of the contaminant variable differs (e.g. Contaminant or Potential contaminant).
Upon converting the peptides table into a QFeatures object spaces in the column names are replaced by a dot.
qf <- filterFeatures(qf, ~ Reverse != "+" & ## Remove decoys
Potential.contaminant != "+") ## Remove contaminants## 'Reverse' found in 2 out of 2 assay(s).'Potential.contaminant' found in 2 out of 2 assay(s).5.3.3 Remove highly missing peptides
The data are characterised by a high proportion of missing values as
reported by nNA() from the QFeatures package. The function
computes the number and the proportion of missing values across
peptides, samples or for the whole data set and return a list of three
tables called nNArows, nNAcols, and nNA, respectively.
nNaRes <- nNA(qf, "peptides")
range(nNaRes$nNAcols$pNA)## [1] 0.4245165 0.4963918The samples contain between 42.5%
and 49.6% missing values. The
missingness within each peptide is more variable, with most peptides
found accross all samples (nNA is 0) and other that could not be
quantified in any sample (nNA is 6), as depicted
by the histogram below.
data.frame(nNaRes$nNArows) |>
ggplot() +
aes(x = nNA) +
geom_histogram() +
ggtitle("Number of missing values for each peptide") +
theme_bw()## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
We keep peptides that were observed at last 3 times out of the \(n = 6\) samples, so that we can estimate the peptide characteristics. We tolerate the following proportion of NAs: \(\text{pNA} = \frac{(n - 3)}{n} = 0.5\), so we keep peptides that are observed in at least 50% of the samples, which corresponds to one treatment condition. This is an arbitrary value that may need to be adjusted depending on the experiment and the data set.
nObs <- 3
n <- ncol(qf[["peptides_log"]])
(qf <- filterNA(qf, i = "peptides_log", pNA = (n - nObs) / n))## An instance of class QFeatures (type: bulk) with 2 sets:
##
## [1] peptides: SummarizedExperiment with 10393 rows and 6 columns
## [2] peptides_log: SummarizedExperiment with 5910 rows and 6 columnsWe keep 5910 peptides upon filtering.
5.4 Normalisation
The most common objective of MS-based proteomics experiments is to understand the biological changes in protein abundance between experimental conditions. However, changes in measurements between groups can be caused due to technical factors. For instance, there are systematic fluctuations from run-to-run that shift the measured intensity distribution. We can this explore as follows:
- We extract the sets containing the log transformed data. This is
performed using
QFeatures’ 3-way subsetting. - We use
longForm()to convert theQFeaturesobject into a long table, including allcolDatafor filtering and colouring. - We visualise the density of the quantitative values within each sample. We colour each sample based on its spike-in condition.
qf[, , "peptides_log"] |> #1.
longForm(colvars = colnames(colData(qf))) |> #2.
data.frame() |>
filter(!is.na(value)) |>
ggplot() + #3.
aes(x = value,
colour = condition,
group = colname) +
geom_density() +
theme_minimal()## Warning: 'experiments' dropped; see 'drops()'## harmonizing input:
## removing 6 sampleMap rows not in names(experiments)
Even in this clean synthetic data set (same yeast background with some spiked UPS proteins), the marginal precursor intensity distributions across samples are not well aligned. Ignoring this effect will increase the noise and reduce the statistical power of the experiment, and may also, in case of unbalanced designs, introduce confounding effects that will bias the results.
There are many ways to perform the normalisation, e.g.
median centering is a popular choice.
If we subtract the sample median at the log scale from each precursor
log2 intensity, this basically boils down to calculating log-ratio’s between
the precursor intensity and its sample median.
Here, we choose to work with offsets, that are chosen in such a way so that
the distributions are centered around the location of a reference sample.
E.g. the median of the sample medians.
Note that users can also use all normalisation functions that are provided in
QFeatures, i.e. by replacing the sweep function by QFeatures::normalize.
qf <- sweep( #Subtract log2 norm factor column-by-column (MARGIN = 2)
qf,
MARGIN = 2,
STATS = nfLogMedian(qf,"peptides_log"),
i = "peptides_log",
name = "peptides_norm"
)## This function aims to calculate norm factors on a log scale,
## the input data are assumed to be on the log-scale!We explore the effect of the global normalisation in the subsequent plot.
Formally, the function applies the following operation on each sample \(i\) across all precursors \(p\):
\[ y_{ip}^{\text{norm}} = y_{ip} - \log_2(nf_i) \]
with \(y_{ip}\) the log2-transformed intensities and \(\log_2(nf_i)\) the log2-transformed norm factor. Upon normalisation, we can see that the distribution of the \(y_{ip}^{\text{norm}}\) nicely overlap (using the same code as above)
qf[, , "peptides_norm"] |>
longForm(colvars = colnames(colData(qf))) |>
data.frame() |>
filter(!is.na(value)) |>
ggplot() +
aes(x = value,
colour = condition,
group = colname) +
geom_density() +
labs(subtitle = "Normalised log2 precursor intensities") +
theme_minimal()## Warning: 'experiments' dropped; see 'drops()'## harmonizing input:
## removing 12 sampleMap rows not in names(experiments)
5.5 Summarization to protein level
Here, we summarise the peptide-level data into protein intensities using
aggregateFeatures(), which streamlines summarisation. It requires the name
of a rowData column to group the peptides into proteins (or
protein groups), here Proteins. We can provide the summarisation
approach through the fun argument.
We use the standard summarisation in aggregateFeatures, which is a
robust summarisation method.
qf <- aggregateFeatures(qf,
i = "peptides_norm",
fcol = "Proteins",
na.rm = TRUE,
name = "proteins"
)##
|
| | 0%##
|
|======================================================================| 100%## The following messages occurred during aggregation:##
## Your quantitative and row data contain missing values. Please read the
## relevant section(s) in the aggregateFeatures manual page regarding the
## effects of missing values on data aggregation.
##
## Occurred during the aggregation of set(s): peptidesNote, that robust summarisation can take a long time for large datasets.
In that case we suggest to use fun = MsCoreUtils::medianPolish, na.rm = TRUE.
Note that the optional argument na.rm = TRUE is needed when using medianPolish
because the data contains missing values as we did not adopt imputation and
medianPolish by default cannot handle missing values.
6 Data exploration and QC
Data exploration aims to highlight the main sources of variation in the data prior to data modelling and can pinpoint to outlying or off-behaving samples.
6.1 Marginal distribution at precursor and protein level
qf[, , "peptides_norm"] |>
longForm(colvars = colnames(colData(qf))) |>
data.frame() |>
filter(!is.na(value)) |>
ggplot() +
aes(x = value,
colour = condition,
group = colname) +
geom_density() +
theme_minimal() +
labs(subtitle = "Normalised log2 precursor intensities")## Warning: 'experiments' dropped; see 'drops()'## harmonizing input:
## removing 18 sampleMap rows not in names(experiments)
qf[, , "peptides_norm"] |>
longForm(colvars = colnames(colData(qf))) |>
data.frame() |>
filter(!is.na(value)) |>
ggplot() +
aes(x = sampleId,
y = value,
colour = condition,
group = colname) +
xlab("sample") +
geom_boxplot() +
theme_minimal() +
labs(subtitle = "Normalised log2 precursor intensities")## Warning: 'experiments' dropped; see 'drops()'## harmonizing input:
## removing 18 sampleMap rows not in names(experiments)
We do the same for the protein-level values.
qf[, , "proteins"] |>
longForm(colvars = colnames(colData(qf))) |>
data.frame() |>
filter(!is.na(value)) |>
ggplot() +
aes(x = value,
colour = condition,
group = colname) +
geom_density() +
theme_minimal() +
labs(subtitle = "Normalised log2 protein intensities")## Warning: 'experiments' dropped; see 'drops()'## harmonizing input:
## removing 18 sampleMap rows not in names(experiments)
qf[, , "proteins"] |>
longForm(colvars = colnames(colData(qf))) |>
data.frame() |>
filter(!is.na(value)) |>
ggplot() +
aes(x = sampleId,
y = value,
colour = condition,
group = colname) +
xlab("sample") +
geom_boxplot() +
theme_minimal() +
labs(subtitle = "Normalised log2 protein intensities")## Warning: 'experiments' dropped; see 'drops()'## harmonizing input:
## removing 18 sampleMap rows not in names(experiments)
6.2 Identifications per sample
qf[,,"peptides_norm"] |>
longForm(colvars = colnames(colData(qf)),
rowvars= c("Sequence",
"Proteins")) |>
data.frame() |>
filter(!is.na(value)) |>
group_by(condition, sampleId) |>
summarise(Precursors = length(unique(Sequence)),
`Protein Groups` = length(unique(Proteins))) |>
pivot_longer(-(1:2),
names_to = "Feature",
values_to = "IDs") |>
ggplot(aes(x = sampleId, y = IDs, fill = condition)) +
geom_col() +
#scale_fill_observable() +
facet_wrap(~Feature,
scales = "free_y") +
labs(title = "Peptide and protein group identificiations per sample",
x = "Sample",
y = "Identifications") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90))## Warning: 'experiments' dropped; see 'drops()'## harmonizing input:
## removing 18 sampleMap rows not in names(experiments)## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by condition and sampleId.
## ℹ Output is grouped by condition.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(condition, sampleId))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.
6.3 Dimensionality reduction plot
A common approach for data exploration is to perform dimension reduction, such as Multi Dimensional Scaling (MDS). We will first extract the set to explore along the sample annotations (used for plot colouring).
getWithColData(qf, "proteins") |>
as("SingleCellExperiment") |>
runMDS(exprs_values = 1) |>
plotMDS(colour_by = "condition", shape_by = "rep", point_size = 3)## Warning: 'experiments' dropped; see 'drops()'
This plot reveals interesting information. First, we see that the
samples are nicely separated according to their spike-in condition.
The also seem to line up according to rep.
6.4 Correlation matrix
corMat <- qf[["proteins"]] |>
assay() |>
cor(method = "spearman", use = "pairwise.complete.obs")
colnames(corMat) <- qf$sampleId
rownames(corMat) <- qf$sampleId
corMat |>
ggcorrplot::ggcorrplot() +
scale_fill_viridis_c() +
labs(title = "Correlation matrix",
fill = "Correlation") +
theme(axis.text.x = element_text(angle = 90))## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the ggcorrplot package.
## Please report the issue at <https://github.com/kassambara/ggcorrplot/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
7 Data Modeling (Robust Regression)
7.1 Sources of variation
For this experiment, we can only model a single source of variation: spike-in concentration. All remaining sources of variability, e.g. pipetting errors, MS run to run variability, etc will be lumped in the residual variability (error term).
Spike-in condition effects: we model the source of variation induced by the experimental treatment of interest as a fixed effect, which we consider non-random, i.e. the treatment effect is assumed to be the same in repeated experiments, but it is unknown and has to be estimated.
When modelling a typical label-free experiment at the protein level, the model boils down to the following linear model:
\[ y_i = \mathbf{x}_i^T\boldsymbol{\beta} +\epsilon_i\]
with the \(\log_2\)-normalised and summarised protein intensities in sample;
\(\mathbf{x}_i\) a vector with the covariate pattern for the sample encoding the
intercept, spike-in condition, or other experimental factors;
\(\boldsymbol{\beta}\) the vector of parameters that model the association between
the covariates and the outcome; and \(\epsilon_i\) the residuals reflecting
variation that is not captured by the fixed effects.
Note that allows for a flexible parameterisation of the treatment
beyond a single covariate, i.e. including a 1 for the intercept, continuous and
categorical variables as well as their interactions. We assume the residuals to
be independent and identically distributed (i.i.d) according to a normal
distribution with zero mean and constant variance,
i.e. \(\epsilon_i \sim N(0, \sigma^2_\epsilon)\), that can differ from protein to
protein.
In R, this model is encoded using the following simple formula:
model <- ~ condition7.2 Model Estimation
We estimate the model with msqrob(). The function takes the QFeatures object,
extracts the quantitative values from the "proteins" set generated during
summarisation, and fits a simple linear model with condition as covariate,
which are automatically retrieved from colData(qf).
Note, that msqrob() can also take a SummarizedExperiment of
precursor, peptide or protein intensities as its input.
qf <- msqrob(
qf,
i = "proteins",
formula = model,
robust = TRUE)We enable M-estimation (robust = TRUE) for improved robustness against outliers.
The fitting results are available in the msqrobModels column of the rowData. More specifically, the modelling output is stored in the rowData as a statModel object, one model per row (protein). We will see in a later section how to perform statistical inference on the estimated parameters.
models <- rowData(qf[["proteins"]])[["msqrobModels"]]
models[1:3]## $`O00762ups|UBE2C_HUMAN_UPS`
## Object of class "StatModel"
## The type is "rlm"
## There number of elements is 5
##
## $`P00167ups|CYB5_HUMAN_UPS`
## Object of class "StatModel"
## The type is "fitError"
## There number of elements is 5
##
## $`P00441ups|SODC_HUMAN_UPS`
## Object of class "StatModel"
## The type is "rlm"
## There number of elements is 57.3 Statistical inference
We can now convert the research question “do the spike-in conditions affect the protein intensities?” into a statistical hypothesis.
In other words, we need to translate this question in a null and alternative hypothesis on a single model parameter or on a linear combination of model parameters, which is also referred to with a contrast.
To aid defining contrasts, we will visualise the experimental design using the ExploreModelMatrix package.
vd <- ExploreModelMatrix::VisualizeDesign(
sampleData = colData(qf),
designFormula = model,
textSizeFitted = 4
)
vd$plotlist## [[1]]
This plot shows that the average protein \(\log_2\) intensity for conditionA is modelled by (Intercept) and that for the conditionB group by (Intercept) + conditionB .
Hence, to assess differential abundance of proteins between the two spike-in conditions, we will have to assess the following null hypothesis: \(\log_2 FC_{B-A} = (Intercept + conditionB) - Intercept = conditionB = 0\).
Based on this null hypothesis we now generate a contrast matrix.
(L <- makeContrast(
"conditionB = 0",
parameterNames = colnames(vd$designmatrix)
))## conditionB
## (Intercept) 0
## conditionB 1We can now test our null hypothesis using hypothesisTest() which
takes the QFeatures object with the fitted model and the contrast we
just built. msqrob2 automatically applies the hypothesis
testing to all features in the assay.
qf <- hypothesisTest(qf, i = "proteins", contrast = L)The results tables for the contrast are stored in the row data of the
proteins assay in qf under the colomnames of the contrast matrix L.
Users can directly access the results tables via de colData. However, msqrob2
also provides the msqrobCollect function. This function will automatically
retrieve the results tables for all contrasts in contrast matrix L and combine
them by default in one table. This table contains two additional columns:
contrast with the name of the contrast and feature with the name of the
modelled feature, here protein.
inferences <-
msqrobCollect(qf[["proteins"]], L)
head(inferences)## logFC se df t
## conditionB.O00762ups|UBE2C_HUMAN_UPS 1.5026502 0.1418246 6.066186 10.595131
## conditionB.P00167ups|CYB5_HUMAN_UPS NA NA NA NA
## conditionB.P00441ups|SODC_HUMAN_UPS -0.4291633 0.1782764 6.117213 -2.407292
## conditionB.P00709ups|LALBA_HUMAN_UPS 1.3257554 0.4567491 6.117213 2.902590
## conditionB.P00915ups|CAH1_HUMAN_UPS NA NA NA NA
## conditionB.P00918ups|CAH2_HUMAN_UPS 1.4592637 0.1961286 7.117213 7.440342
## pval adjPval contrast
## conditionB.O00762ups|UBE2C_HUMAN_UPS 3.873132e-05 0.004522272 conditionB
## conditionB.P00167ups|CYB5_HUMAN_UPS NA NA conditionB
## conditionB.P00441ups|SODC_HUMAN_UPS 5.198259e-02 0.893152586 conditionB
## conditionB.P00709ups|LALBA_HUMAN_UPS 2.664631e-02 0.668425430 conditionB
## conditionB.P00915ups|CAH1_HUMAN_UPS NA NA conditionB
## conditionB.P00918ups|CAH2_HUMAN_UPS 1.329479e-04 0.010449706 conditionB
## feature
## conditionB.O00762ups|UBE2C_HUMAN_UPS O00762ups|UBE2C_HUMAN_UPS
## conditionB.P00167ups|CYB5_HUMAN_UPS P00167ups|CYB5_HUMAN_UPS
## conditionB.P00441ups|SODC_HUMAN_UPS P00441ups|SODC_HUMAN_UPS
## conditionB.P00709ups|LALBA_HUMAN_UPS P00709ups|LALBA_HUMAN_UPS
## conditionB.P00915ups|CAH1_HUMAN_UPS P00915ups|CAH1_HUMAN_UPS
## conditionB.P00918ups|CAH2_HUMAN_UPS P00918ups|CAH2_HUMAN_UPS7.3.1 Results table for significant proteins
We will return a table of proteins for which the contrast is significant at the 5% FDR level.
- We set the significance level
- We filter the significant results for the contrast from the table
- We print the table
alpha <- 0.05 #1.
sigList <- inferences |>
filter(adjPval < alpha) #3.
kable(sigList) |>
kable_styling(full_width = FALSE) |>
scroll_box(height = "250px") #4.| logFC | se | df | t | pval | adjPval | contrast | feature | |
|---|---|---|---|---|---|---|---|---|
| conditionB.O00762ups|UBE2C_HUMAN_UPS | 1.502650 | 0.1418246 | 6.066186 | 10.595131 | 0.0000387 | 0.0045223 | conditionB | O00762ups|UBE2C_HUMAN_UPS |
| conditionB.P00918ups|CAH2_HUMAN_UPS | 1.459264 | 0.1961286 | 7.117213 | 7.440342 | 0.0001329 | 0.0104497 | conditionB | P00918ups|CAH2_HUMAN_UPS |
| conditionB.P01031ups|CO5_HUMAN_UPS | 1.744301 | 0.1856229 | 5.528016 | 9.397013 | 0.0001323 | 0.0104497 | conditionB | P01031ups|CO5_HUMAN_UPS |
| conditionB.P01127ups|PDGFB_HUMAN_UPS | 1.475470 | 0.1622265 | 6.836031 | 9.095124 | 0.0000460 | 0.0045223 | conditionB | P01127ups|PDGFB_HUMAN_UPS |
| conditionB.P01344ups|IGF2_HUMAN_UPS | 1.481863 | 0.1289776 | 6.117213 | 11.489306 | 0.0000228 | 0.0033600 | conditionB | P01344ups|IGF2_HUMAN_UPS |
| conditionB.P01375ups|TNFA_HUMAN_UPS | 1.858881 | 0.2319031 | 6.117213 | 8.015766 | 0.0001824 | 0.0134425 | conditionB | P01375ups|TNFA_HUMAN_UPS |
| conditionB.P02144ups|MYG_HUMAN_UPS | 1.367831 | 0.2173194 | 7.117213 | 6.294107 | 0.0003799 | 0.0248859 | conditionB | P02144ups|MYG_HUMAN_UPS |
| conditionB.P02787ups|TRFE_HUMAN_UPS | 1.462212 | 0.1518081 | 6.604679 | 9.631979 | 0.0000397 | 0.0045223 | conditionB | P02787ups|TRFE_HUMAN_UPS |
| conditionB.P06732ups|KCRM_HUMAN_UPS | 1.615841 | 0.1315682 | 7.088294 | 12.281393 | 0.0000049 | 0.0014503 | conditionB | P06732ups|KCRM_HUMAN_UPS |
| conditionB.P07339ups|CATD_HUMAN_UPS | 1.901800 | 0.1838925 | 6.117213 | 10.341911 | 0.0000422 | 0.0045223 | conditionB | P07339ups|CATD_HUMAN_UPS |
| conditionB.P08311ups|CATG_HUMAN_UPS | 1.543161 | 0.1117913 | 7.117213 | 13.803946 | 0.0000021 | 0.0008392 | conditionB | P08311ups|CATG_HUMAN_UPS |
| conditionB.P09211ups|GSTP1_HUMAN_UPS | 1.171263 | 0.1876560 | 6.006373 | 6.241542 | 0.0007802 | 0.0459900 | conditionB | P09211ups|GSTP1_HUMAN_UPS |
| conditionB.P10636-8ups|TAU_HUMAN_UPS | 1.264263 | 0.1138747 | 7.117213 | 11.102235 | 0.0000095 | 0.0016011 | conditionB | P10636-8ups|TAU_HUMAN_UPS |
| conditionB.P12081ups|SYHC_HUMAN_UPS | 1.411765 | 0.1810004 | 6.117213 | 7.799792 | 0.0002128 | 0.0147581 | conditionB | P12081ups|SYHC_HUMAN_UPS |
| conditionB.P51965ups|UB2E1_HUMAN_UPS | 1.646747 | 0.1484305 | 7.117213 | 11.094395 | 0.0000095 | 0.0016011 | conditionB | P51965ups|UB2E1_HUMAN_UPS |
| conditionB.P61769ups|B2MG_HUMAN_UPS | 2.007500 | 0.3295815 | 7.029565 | 6.091059 | 0.0004874 | 0.0302429 | conditionB | P61769ups|B2MG_HUMAN_UPS |
| conditionB.P63165ups|SUMO1_HUMAN_UPS | 1.044336 | 0.1317043 | 7.117213 | 7.929393 | 0.0000883 | 0.0080056 | conditionB | P63165ups|SUMO1_HUMAN_UPS |
| conditionB.P63279ups|UBC9_HUMAN_UPS | 1.804380 | 0.1170894 | 7.117213 | 15.410281 | 0.0000010 | 0.0008392 | conditionB | P63279ups|UBC9_HUMAN_UPS |
| conditionB.P99999ups|CYC_HUMAN_UPS | 2.043480 | 0.1462042 | 7.117213 | 13.976892 | 0.0000020 | 0.0008392 | conditionB | P99999ups|CYC_HUMAN_UPS |
| conditionB.Q06830ups|PRDX1_HUMAN_UPS | 1.564310 | 0.1362011 | 7.114248 | 11.485293 | 0.0000075 | 0.0016011 | conditionB | Q06830ups|PRDX1_HUMAN_UPS |
7.3.2 Volcanoplots
A volcano plot is a common visualisation that provides an overview of
the hypothesis testing results, plotting the \(-\log_{10}\) p-value1 Note,
that upon this transformation a value of 1 represents a p-value of 0.1, 2
a p-value of 0.01, 3 a p-value of 0.001, etc. as a function of
the estimated log fold change. Volcano plots can be used to highlight
the most interesting proteins that have large fold changes and/or are
highly significant. We can use the table above directly to build a
volcano plot using ggplot2 functionality. We also highlight which
proteins are UPS standards, known to be differentially abundant by
experimental design.
Experienced users can make the plot themselves. However, in msqrob2 we also
provide the plotVolcano function that generates volcanoplots based on msqrob2
inference tables generated with the hypothesisTest function.
inferences |>
plotVolcano()
7.3.3 Heatmaps
We construct heatmaps for the significant proteins for each contrast.
- We get the names of the significant features.
- We extract the quant data for the significant features and scale it.
- We extract annotation
- We make the heatmap using the quants data. We do not cluster columns (samples) to keep them together according to the design.
sig <- sigList |> pull(feature) #1.
quants <- assay(qf,"proteins")[sig,] |>
t() |>
scale() |>
t() #2.
colnames(quants) <- qf$sampleId #3.
annotations <- columnAnnotation(condition = qf$condition) #3.
set.seed(1234) ## annotation colours are randomly generated by default
Heatmap(
quants,
name = "log2 intensity",
top_annotation = annotations
) #4.
Note, that 20 DA proteins are returned and that they are all UPS proteins!
A typical workflow for the analysis of a proteomics experiment will end here.
8 Assess performance
Note, that the evaluations in this section are typically not possible on real experimental data. Indeed, we are using a spike-in dataset so we know the ground truth: all human UPS proteins are differentially abundant (DA) between the conditions and the yeast proteins are non-DA.
We use a nominal FDR level of 0.05
alpha <- 0.058.1 Real Fold changes
As this is a spike-in study with known ground truth, we can also plot the log2 fold change distributions against the expected values, in this case 0 for the yeast proteins and the difference of the log concentration for the spiked-in UPS standards, i.e. \(logFC = \log_2 (0.74 fmol/0.25 fmol) = \log_2 2.96 = 1.57\).
realLogFC <- log2(0.74/0.25)We collect the inference table again and now add the species information.
inferences <-
msqrobCollect(qf[["proteins"]], L) |>
mutate(species = rowData(qf[["proteins"]])$species)We now evaluate if the estimated fold changes correspond to the real fold changes.
- We plot the log fold changes in function of the spikein condition and color according to spikein condition.
- We add a boxplot layer.
- We use custom colors.
- We add a horizontal line at 0, which corresponds to the known log2 fold change for yeast proteins
- We add a horizontal lines for known log2 fold changes of the spiked UPS proteins.
logFC <- inferences |>
filter(!is.na(logFC)) |>
ggplot(aes(x = species, y = logFC, color = species)) + #1.
geom_boxplot() + #2.
theme_bw() +
scale_color_manual(
values = c("grey20", "firebrick"), #3.
name = "",
labels = c("yeast", "ups")
) +
geom_hline(yintercept = 0, color="grey20") + # 4.
geom_hline(yintercept = realLogFC, color="firebrick") #5.
logFC
Note, that the msqrob2 estimates are nicely distributed around the true fold changes.
8.2 True and false positives
All human UPS proteins are differentially abundant (DA) between the conditions and the yeast proteins are non-DA.
inferences |>
filter(adjPval < alpha) |>
summarise("TP" = sum(species == "ups"),
"FP" = sum(species != "ups"),
"FDP" = mean(species != "ups"))## TP FP FDP
## 1 20 0 0Note, that at the 5% FDR level 20 true positives (TP, UPS proteins) are returned as DA, 0 false positives (FP, yeast proteins) are returned DA, resulting in a false discovery proportion (FDP) of 0, confirming that msqrob2 is providing strong FDR control for this dataset.
8.3 TPR - FDP curves
We generate the TPR-FDP curves to assess the performance to prioritise differentially abundant proteins. Again, these curves are built using the ground truth information about which proteins are differentially abundant (spiked in) and which proteins are constant across samples. We create two functions to compute the TPR and the FDP.
computeFDP <- function(pval, tp) {
ord <- order(pval)
fdp <- cumsum(!tp[ord]) / 1:length(tp)
fdp[order(ord)]
}
computeTPR <- function(pval, tp, nTP = NULL) {
if (is.null(nTP)) nTP <- sum(tp)
ord <- order(pval)
tpr <- cumsum(tp[ord]) / nTP
tpr[order(ord)]
}We apply these functions and compute the corresponding metric using the statistical inference results and the ground truth information.
performance <- inferences |> group_by(contrast) |>
na.exclude() |>
mutate(tpr = computeTPR(pval, species == "ups"),
fdp = computeFDP(pval, species == "ups")) |>
arrange(fdp)We also highlight the observed FDP at a \(alpha = 5\%\) FDR threshold.
workPoints <- performance |>
group_by(contrast) |>
filter(adjPval < 0.05) |>
slice_max(pval)ggplot(performance) +
aes(
y = fdp,
x = tpr,
) +
geom_line() +
geom_point(data = workPoints, size = 3) +
geom_hline(yintercept = 0.05, linetype = 2) +
coord_flip(ylim = c(0, 0.2)) +
theme_minimal()
9 References
Goeminne, L. J. E., A. Sticker, L. Martens, K. Gevaert, and L. Clement. 2020. “MSqRob Takes the Missing Hurdle: Uniting Intensity- and Count-Based Proteomics.” Anal Chem 92 (9): 6278–87.
Goeminne, L. J., K. Gevaert, and L. Clement. 2016. “Peptide-level Robust Ridge Regression Improves Estimation, Sensitivity, and Specificity in Data-dependent Quantitative Label-free Shotgun Proteomics.” Mol Cell Proteomics 15 (2): 657–68.
Sticker, A., L. Goeminne, L. Martens, and L. Clement. 2020. “Robust Summarization and Inference in Proteome-wide Label-free Quantification.” Mol Cell Proteomics 19 (7): 1209–19.
Vandenbulcke, Stijn, Christophe Vanderaa, Oliver Crook, Lennart Martens, and Lieven Clement. 2025. “Msqrob2TMT: Robust Linear Mixed Models for Inferring Differential Abundant Proteins in Labeled Experiments with Arbitrarily Complex Design.” Mol. Cell. Proteomics 24 (7): 101002.
Appendix
Note, that more examples can be found in our msqrob2 book.
sessionInfo()## R version 4.6.0 alpha (2026-04-05 r89794)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.23-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [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: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] grid stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] scater_1.39.4 scuttle_1.21.6
## [3] SingleCellExperiment_1.33.2 ComplexHeatmap_2.27.1
## [5] kableExtra_1.4.0 ExploreModelMatrix_1.23.1
## [7] stringr_1.6.0 msqrob2_1.19.1
## [9] ggplot2_4.0.2 tidyr_1.3.2
## [11] dplyr_1.2.1 QFeatures_1.21.2
## [13] MultiAssayExperiment_1.37.4 SummarizedExperiment_1.41.1
## [15] Biobase_2.71.0 GenomicRanges_1.63.2
## [17] Seqinfo_1.1.0 IRanges_2.45.0
## [19] S4Vectors_0.49.1 BiocGenerics_0.57.0
## [21] generics_0.1.4 MatrixGenerics_1.23.0
## [23] matrixStats_1.5.0 BiocStyle_2.39.0
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 rstudioapi_0.18.0 jsonlite_2.0.0
## [4] shape_1.4.6.1 magrittr_2.0.5 magick_2.9.1
## [7] ggbeeswarm_0.7.3 farver_2.1.2 nloptr_2.2.1
## [10] rmarkdown_2.31 GlobalOptions_0.1.4 vctrs_0.7.3
## [13] Cairo_1.7-0 minqa_1.2.8 tinytex_0.59
## [16] BiocBaseUtils_1.13.0 htmltools_0.5.9 S4Arrays_1.11.1
## [19] BiocNeighbors_2.5.4 SparseArray_1.11.13 sass_0.4.10
## [22] bslib_0.10.0 htmlwidgets_1.6.4 plyr_1.8.9
## [25] cachem_1.1.0 igraph_2.2.3 mime_0.13
## [28] lifecycle_1.0.5 iterators_1.0.14 pkgconfig_2.0.3
## [31] rsvd_1.0.5 Matrix_1.7-5 R6_2.6.1
## [34] fastmap_1.2.0 rbibutils_2.4.1 shiny_1.13.0
## [37] clue_0.3-68 digest_0.6.39 colorspace_2.1-2
## [40] irlba_2.3.7 textshaping_1.0.5 beachmat_2.27.5
## [43] labeling_0.4.3 abind_1.4-8 compiler_4.6.0
## [46] bit64_4.6.0-1 withr_3.0.2 doParallel_1.0.17
## [49] S7_0.2.1 ggcorrplot_0.1.4.1 BiocParallel_1.45.0
## [52] viridis_0.6.5 MASS_7.3-65 DelayedArray_0.37.1
## [55] rjson_0.2.23 tools_4.6.0 vipor_0.4.7
## [58] otel_0.2.0 beeswarm_0.4.0 httpuv_1.6.17
## [61] glue_1.8.0 nlme_3.1-169 promises_1.5.0
## [64] cluster_2.1.8.2 reshape2_1.4.5 gtable_0.3.6
## [67] data.table_1.18.2.1 BiocSingular_1.27.1 ScaledMatrix_1.19.0
## [70] xml2_1.5.2 XVector_0.51.0 ggrepel_0.9.8
## [73] foreach_1.5.2 pillar_1.11.1 limma_3.67.1
## [76] later_1.4.8 rintrojs_0.3.4 circlize_0.4.18
## [79] splines_4.6.0 lattice_0.22-9 bit_4.6.0
## [82] tidyselect_1.2.1 knitr_1.51 gridExtra_2.3
## [85] reformulas_0.4.4 bookdown_0.46 ProtGenerics_1.43.0
## [88] svglite_2.2.2 xfun_0.57 shinydashboard_0.7.3
## [91] statmod_1.5.1 DT_0.34.0 stringi_1.8.7
## [94] lazyeval_0.2.3 yaml_2.3.12 boot_1.3-32
## [97] evaluate_1.0.5 codetools_0.2-20 MsCoreUtils_1.23.9
## [100] tibble_3.3.1 BiocManager_1.30.27 cli_3.6.6
## [103] xtable_1.8-8 systemfonts_1.3.2 Rdpack_2.6.6
## [106] jquerylib_0.1.4 dichromat_2.0-0.1 Rcpp_1.1.1
## [109] png_0.1-9 parallel_4.6.0 assertthat_0.2.1
## [112] AnnotationFilter_1.35.0 lme4_2.0-1 viridisLite_0.4.3
## [115] scales_1.4.0 purrr_1.2.2 crayon_1.5.3
## [118] GetoptLong_1.1.1 rlang_1.2.0 cowplot_1.2.0
## [121] shinyjs_2.1.1