Tutorial 1 - Cancer 9x9

Author

Affiliation

Lieven Clement

Ghent University

1 Load packages

We load the msqrob2 package, along with additional packages for data manipulation and visualisation.

library("QFeatures")  
library("dplyr") 
library("tidyr")
library("ggplot2")
library("msqrob2")    
library("stringr")
library("ExploreModelMatrix")
library("MsCoreUtils") 
library("matrixStats")
library("patchwork")
library("kableExtra")
library("ComplexHeatmap")
library("purrr")
library("tibble")
library("scater")

2 Load QFeatures object

(
  qf <- readRDS(url("https://github.com/statOmics/PDA-DIA/raw/refs/heads/main/data/cancer9x9.rds","rb"))
)

An instance of class QFeatures (type: bulk) with 4 sets:

 [1] peptides_raw: SummarizedExperiment with 26531 rows and 18 columns 
 [2] peptides_log: SummarizedExperiment with 26531 rows and 18 columns 
 [3] peptides_norm: SummarizedExperiment with 26531 rows and 18 columns 
 [4] proteins: SummarizedExperiment with 3506 rows and 18 columns

3 Data exploration

First assess distribution of marginal protein intensities.

qf[, , "proteins"] |>
    longForm(colvars = colnames(colData(qf))) |> 
  data.frame() |> 
  filter(!is.na(value)) |>
  ggplot() + 
  aes(x = value,
      colour = prognosis,
      group = colname) +
  geom_density() +
  theme_minimal() +
   labs(subtitle = "Normalised log2 protein intensities")

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.

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).

qf |> 
  getWithColData("proteins") |> 
  as("SingleCellExperiment") |> 
  scater::runMDS(exprs_values = 1) |> 
  scater::plotMDS(colour_by = "prognosis")

4 Data Modeling (Robust Regression)

4.1 Model estimation

We first define the model. We only have one sources of variability in the experiment that we can model, i.e. the association with prognosis.
We fit the model to each protein in the proteins summarised experiment of the QFeatures object qf using the msqrob function.

model <- ~ prognosis

qf <- msqrob(
  qf,
  i = "proteins",
  formula = model,
  robust = TRUE)

We enabled 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]

$A0AV96
Object of class "StatModel"
The type is "rlm"
There number of elements is 7 

$A0AVT1
Object of class "StatModel"
The type is "rlm"
There number of elements is 7 

$A0FGR8
Object of class "StatModel"
The type is "rlm"
There number of elements is 7

4.2 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 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]]

contrast <- "prognosisPD = 0"
(L <- makeContrast(
   contrast,
   parameterNames = colnames(vd$designmatrix)
))

            prognosisPD
(Intercept)           0
prognosisPD           1

We assess the contrast for each protein.

qf <- hypothesisTest(qf, i = "proteins", contrast = L)

We extract the results table from the proteins summarised experiment in the qf object.

inference <- 
  msqrobCollect(qf[["proteins"]], L)

4.3 Report results

We report the results using a results table, volcano plots and heatmaps.

4.3.1 Results table

We report all results that are significant at the nominal FDR level of 5%.

We define the nominal FDR level
We filter the results
We arrange them according to statistical significance

alpha  <- 0.05
inference |> 
  filter(adjPval < alpha) |> 
  arrange(pval) |> 
  knitr::kable()

	logFC	se	df	t	pval	adjPval	contrast	feature
prognosisPD.P61604	2.6041240	0.3250551	20.64962	8.011330	0.0000001	0.0002967	prognosisPD	P61604
prognosisPD.P06454	3.3789333	0.4510993	20.88319	7.490443	0.0000002	0.0003938	prognosisPD	P06454
prognosisPD.Q96K17	1.6766684	0.2135589	16.66846	7.851081	0.0000005	0.0003938	prognosisPD	Q96K17
prognosisPD.P49006	2.2746124	0.3101449	18.22272	7.334031	0.0000008	0.0003938	prognosisPD	P49006
prognosisPD.Q8NC51	1.6657595	0.2363348	19.80273	7.048305	0.0000008	0.0003938	prognosisPD	Q8NC51
prognosisPD.Q13765	1.4199990	0.2040024	19.68369	6.960697	0.0000010	0.0003938	prognosisPD	Q13765
prognosisPD.O43396	1.9230797	0.2713492	18.75705	7.087104	0.0000010	0.0003938	prognosisPD	O43396
prognosisPD.P54727	1.6173693	0.2316137	19.14945	6.983046	0.0000011	0.0003938	prognosisPD	P54727
prognosisPD.Q9NXG2	1.7143795	0.2540731	20.76058	6.747584	0.0000012	0.0003938	prognosisPD	Q9NXG2
prognosisPD.P07108	4.2542023	0.4749194	11.93903	8.957735	0.0000012	0.0003938	prognosisPD	P07108
prognosisPD.Q9H299	1.9395700	0.2921190	20.03375	6.639657	0.0000018	0.0005377	prognosisPD	Q9H299
prognosisPD.Q9GZT3	2.0157004	0.3200318	20.42718	6.298438	0.0000034	0.0009349	prognosisPD	Q9GZT3
prognosisPD.P63313	3.3989709	0.4818571	14.60186	7.053898	0.0000045	0.0011339	prognosisPD	P63313
prognosisPD.Q9UBE0	1.1997904	0.1964488	20.65961	6.107395	0.0000050	0.0011339	prognosisPD	Q9UBE0
prognosisPD.P27797	1.4613961	0.2399652	20.63039	6.090033	0.0000052	0.0011339	prognosisPD	P27797
prognosisPD.O60635	-1.8704902	0.2634919	12.69920	-7.098852	0.0000092	0.0018736	prognosisPD	O60635
prognosisPD.P09496	1.4512161	0.2534685	20.53051	5.725430	0.0000120	0.0022366	prognosisPD	P09496
prognosisPD.P09012	1.8807719	0.3119256	17.37835	6.029553	0.0000123	0.0022366	prognosisPD	P09012
prognosisPD.O43765	1.0501577	0.1847814	20.45672	5.683244	0.0000134	0.0023090	prognosisPD	O43765
prognosisPD.P30084	1.7956099	0.3160253	19.49759	5.681854	0.0000161	0.0026326	prognosisPD	P30084
prognosisPD.O75368	2.2473966	0.4029769	18.64520	5.576987	0.0000239	0.0036002	prognosisPD	O75368
prognosisPD.P55327	1.7147764	0.3130844	19.58030	5.477042	0.0000250	0.0036002	prognosisPD	P55327
prognosisPD.O75223	1.8772711	0.3465195	20.27451	5.417506	0.0000253	0.0036002	prognosisPD	O75223
prognosisPD.P30085	1.2350160	0.2354281	20.82664	5.245830	0.0000345	0.0046977	prognosisPD	P30085
prognosisPD.P46108	1.6135876	0.2911523	16.93903	5.542075	0.0000363	0.0047460	prognosisPD	P46108
prognosisPD.P22307	1.3504567	0.2562669	19.03567	5.269727	0.0000434	0.0053024	prognosisPD	P22307
prognosisPD.P20042	1.0102350	0.1962977	20.48862	5.146443	0.0000457	0.0053024	prognosisPD	P20042
prognosisPD.P43490	1.4879350	0.2846566	19.24274	5.227123	0.0000461	0.0053024	prognosisPD	P43490
prognosisPD.Q15758	-1.4777359	0.2760829	17.59084	-5.352508	0.0000471	0.0053024	prognosisPD	Q15758
prognosisPD.P06730	2.0005445	0.3485479	13.91310	5.739655	0.0000524	0.0056267	prognosisPD	P06730
prognosisPD.Q15370	1.4689650	0.2842868	19.16415	5.167193	0.0000534	0.0056267	prognosisPD	Q15370
prognosisPD.P51858	1.8539942	0.3597673	19.10748	5.153315	0.0000556	0.0056731	prognosisPD	P51858
prognosisPD.Q9NRX4	1.8243353	0.3541559	17.91090	5.151222	0.0000680	0.0067376	prognosisPD	Q9NRX4
prognosisPD.P19338	1.1055554	0.2215553	18.76278	4.989975	0.0000843	0.0080987	prognosisPD	P19338
prognosisPD.Q16543	1.0422066	0.2115207	19.23078	4.927209	0.0000905	0.0084491	prognosisPD	Q16543
prognosisPD.O95571	1.5581033	0.3002873	15.06486	5.188709	0.0001086	0.0098583	prognosisPD	O95571
prognosisPD.P52434	1.0670443	0.2241178	20.17216	4.761086	0.0001167	0.0103104	prognosisPD	P52434
prognosisPD.Q9UMX5	1.4639249	0.2972751	15.93903	4.924478	0.0001541	0.0131575	prognosisPD	Q9UMX5
prognosisPD.P20962	2.7420169	0.5557565	15.73127	4.933846	0.0001570	0.0131575	prognosisPD	P20962
prognosisPD.P20810	1.4349529	0.3012980	17.59555	4.762570	0.0001651	0.0134862	prognosisPD	P20810
prognosisPD.Q15185	1.0867276	0.2377714	20.73930	4.570472	0.0001709	0.0136255	prognosisPD	Q15185
prognosisPD.Q14980	-1.0853277	0.2375710	19.79668	-4.568436	0.0001910	0.0148256	prognosisPD	Q14980
prognosisPD.Q99497	1.3209141	0.2921426	20.59083	4.521471	0.0001951	0.0148256	prognosisPD	Q99497
prognosisPD.P07858	1.3256342	0.2867481	17.97469	4.622992	0.0002120	0.0157055	prognosisPD	P07858
prognosisPD.O60888	1.1159069	0.2433521	18.45467	4.585565	0.0002163	0.0157055	prognosisPD	O60888
prognosisPD.P02787	1.5014232	0.3372296	20.08783	4.452229	0.0002422	0.0170832	prognosisPD	P02787
prognosisPD.Q9Y2Q5	1.1019938	0.2446334	18.89322	4.504675	0.0002457	0.0170832	prognosisPD	Q9Y2Q5
prognosisPD.P61086	0.9821935	0.2233218	20.64657	4.398108	0.0002597	0.0176822	prognosisPD	P61086
prognosisPD.P12004	1.1987525	0.2729536	19.63179	4.391782	0.0002927	0.0195206	prognosisPD	P12004
prognosisPD.Q15181	1.7855711	0.4021997	18.44766	4.439514	0.0002999	0.0196047	prognosisPD	Q15181
prognosisPD.P49321	1.1268135	0.2598270	20.30155	4.336784	0.0003107	0.0199062	prognosisPD	P49321
prognosisPD.Q92734	1.0659326	0.2461375	19.58764	4.330639	0.0003387	0.0210866	prognosisPD	Q92734
prognosisPD.Q9HB71	1.0586501	0.2454583	19.89652	4.312953	0.0003420	0.0210866	prognosisPD	Q9HB71
prognosisPD.Q9NR45	1.4417491	0.3380157	20.64347	4.265332	0.0003559	0.0215410	prognosisPD	Q9NR45
prognosisPD.Q9Y266	0.9371062	0.2177031	19.40272	4.304514	0.0003667	0.0217177	prognosisPD	Q9Y266
prognosisPD.Q9UHV9	0.8842038	0.2055458	19.22434	4.301736	0.0003759	0.0217177	prognosisPD	Q9UHV9
prognosisPD.P14314	0.8439518	0.1997391	20.66822	4.225271	0.0003906	0.0217177	prognosisPD	P14314
prognosisPD.P62318	0.7907500	0.1873286	20.77508	4.221192	0.0003906	0.0217177	prognosisPD	P62318
prognosisPD.P10412	1.4688990	0.3269031	15.43337	4.493378	0.0004010	0.0217177	prognosisPD	P10412
prognosisPD.O15327	1.3284846	0.2974321	15.73944	4.466514	0.0004047	0.0217177	prognosisPD	O15327
prognosisPD.Q02978	-1.0853069	0.2569914	20.31609	-4.223126	0.0004054	0.0217177	prognosisPD	Q02978
prognosisPD.Q06830	0.7891836	0.1863979	19.77673	4.233865	0.0004160	0.0219255	prognosisPD	Q06830
prognosisPD.P04080	1.4100115	0.3345481	19.93148	4.214675	0.0004286	0.0220237	prognosisPD	P04080
prognosisPD.P52565	1.5729187	0.3731268	19.84467	4.215507	0.0004313	0.0220237	prognosisPD	P52565
prognosisPD.P29966	1.5913581	0.3757393	18.76673	4.235272	0.0004587	0.0230601	prognosisPD	P29966
prognosisPD.Q9NP79	0.8685699	0.2097729	20.52359	4.140525	0.0004836	0.0239456	prognosisPD	Q9NP79
prognosisPD.Q05682	1.3946138	0.3325873	18.70698	4.193226	0.0005076	0.0247587	prognosisPD	Q05682
prognosisPD.P19971	1.9567252	0.4831379	20.47913	4.050035	0.0006013	0.0288964	prognosisPD	P19971
prognosisPD.O43237	0.8678016	0.2144262	19.21564	4.047088	0.0006748	0.0319594	prognosisPD	O43237
prognosisPD.Q16513	1.9939586	0.4289551	10.93903	4.648409	0.0007174	0.0334637	prognosisPD	Q16513
prognosisPD.P61960	1.1897596	0.2949624	18.59383	4.033597	0.0007365	0.0334637	prognosisPD	P61960
prognosisPD.P05387	1.1304980	0.2803130	18.59760	4.032984	0.0007373	0.0334637	prognosisPD	P05387
prognosisPD.Q13451	1.6160995	0.3955798	17.24501	4.085394	0.0007507	0.0336069	prognosisPD	Q13451
prognosisPD.Q5T4S7	-0.8922697	0.2259228	20.28039	-3.949445	0.0007745	0.0340218	prognosisPD	Q5T4S7
prognosisPD.A6NL28	0.9452101	0.2406413	20.64733	3.927880	0.0007923	0.0340218	prognosisPD	A6NL28
prognosisPD.Q9UHQ9	-0.9879144	0.2509873	20.26383	-3.936113	0.0008001	0.0340218	prognosisPD	Q9UHQ9
prognosisPD.Q6IA17	-0.9732841	0.2330255	15.07899	-4.176727	0.0008016	0.0340218	prognosisPD	Q6IA17
prognosisPD.Q9BT22	-0.9597435	0.2368722	17.15594	-4.051735	0.0008154	0.0341632	prognosisPD	Q9BT22
prognosisPD.P22061	1.1696527	0.2972531	19.79422	3.934871	0.0008327	0.0341905	prognosisPD	P22061
prognosisPD.Q15293	1.9064446	0.4853782	19.93903	3.927751	0.0008370	0.0341905	prognosisPD	Q15293
prognosisPD.Q02818	1.0734026	0.2755959	20.77630	3.894842	0.0008488	0.0342438	prognosisPD	Q02818
prognosisPD.Q9NZJ7	-1.0037751	0.2383431	13.93903	-4.211472	0.0008789	0.0348443	prognosisPD	Q9NZJ7
prognosisPD.P07741	1.3601391	0.3478872	19.76498	3.909713	0.0008850	0.0348443	prognosisPD	P07741
prognosisPD.Q9HAB8	0.8434880	0.2183417	20.93903	3.863156	0.0009046	0.0351643	prognosisPD	Q9HAB8
prognosisPD.P21333	1.3454227	0.3444733	19.46718	3.905738	0.0009146	0.0351643	prognosisPD	P21333
prognosisPD.Q15366	0.9534380	0.2466510	19.73583	3.865535	0.0009829	0.0369669	prognosisPD	Q15366
prognosisPD.Q92522	1.5742635	0.3722903	13.00952	4.228591	0.0009841	0.0369669	prognosisPD	Q92522
prognosisPD.P04040	1.0612352	0.2769986	20.64827	3.831194	0.0009961	0.0369909	prognosisPD	P04040
prognosisPD.Q9H2U2	1.1891812	0.3086989	18.10258	3.852237	0.0011570	0.0414649	prognosisPD	Q9H2U2
prognosisPD.Q96KP4	1.0514749	0.2783308	20.00116	3.777789	0.0011821	0.0414649	prognosisPD	Q96KP4
prognosisPD.Q6IAA8	-0.8133843	0.2146964	19.57791	-3.788532	0.0011889	0.0414649	prognosisPD	Q6IAA8
prognosisPD.Q8TE68	1.1278717	0.2777013	13.85091	4.061456	0.0011911	0.0414649	prognosisPD	Q8TE68
prognosisPD.P62310	1.3260664	0.3356226	15.60247	3.951064	0.0011951	0.0414649	prognosisPD	P62310
prognosisPD.Q6I9Y2	2.4577502	0.5661627	10.88284	4.341067	0.0012022	0.0414649	prognosisPD	Q6I9Y2
prognosisPD.P30086	-1.3931548	0.3687186	19.59626	-3.778369	0.0012155	0.0414649	prognosisPD	P30086
prognosisPD.P82650	-0.7594181	0.2024826	20.49157	-3.750536	0.0012181	0.0414649	prognosisPD	P82650
prognosisPD.P12268	0.8651493	0.2306933	20.12625	3.750214	0.0012497	0.0421045	prognosisPD	P12268
prognosisPD.O95471	-1.8078980	0.4573136	14.76823	-3.953300	0.0013106	0.0437060	prognosisPD	O95471
prognosisPD.O43399	1.1930557	0.3087231	16.12812	3.864484	0.0013549	0.0447260	prognosisPD	O43399
prognosisPD.O43583	0.8889550	0.2398751	20.19930	3.705908	0.0013793	0.0450750	prognosisPD	O43583
prognosisPD.P02790	1.3386832	0.3626135	19.66596	3.691763	0.0014777	0.0471062	prognosisPD	P02790
prognosisPD.P16401	1.8920379	0.5048054	17.93903	3.748054	0.0014796	0.0471062	prognosisPD	P16401
prognosisPD.Q15126	1.1932211	0.3214274	18.85252	3.712257	0.0014940	0.0471062	prognosisPD	Q15126
prognosisPD.P67936	1.4828373	0.4027827	19.80493	3.681482	0.0014991	0.0471062	prognosisPD	P67936
prognosisPD.P31937	0.7330570	0.1969090	17.86313	3.722822	0.0015746	0.0487991	prognosisPD	P31937
prognosisPD.P16949	1.7079891	0.4534750	16.68732	3.766446	0.0015828	0.0487991	prognosisPD	P16949

We observe that the results only contain spiked UPS proteins.

4.3.2 Volcanoplots

volcanoplot4 <- plot_volcano(inference)
volcanoplot4

Note, that the log fold changes are nicely around the real log2 fold change: \(log2(4/2) = 1\)

4.3.3 Heatmaps

We make the heatmap as follows

We select the names of the proteins that were declared significant between condition A and condition B and extract their quantitative data.
We extract the quantitative data with assay() and scale by rows.
We will create a heatmap using the ComplexHeatmap package, which enables heatmap annotations. We will annotate the heatmap using our model variables condition and lab.
We make the heatmap

sig <- inference |> 
  filter(adjPval < alpha) |> 
  arrange(pval) |> 
  pull(feature) #1.

se <- getWithColData(qf, "proteins")
quants <- t(scale(t(assay(se[sig,])))) #2. 

annotations <- columnAnnotation(
  condition = se$prognosis
) #3.

set.seed(1234) ## annotation colours are randomly generated by default
heatmap <- Heatmap(
                       quants, 
                       name = "log2 intensity",
                       top_annotation = annotations, 
                       column_title = paste0(contrast, " = 0")
                       ) #4.
heatmap

4.3.4 Detail plots

We can explore the data for a protein to validate the statistical inference results. For example, let’s explore the normalised precursor and the summarised protein intensities for the protein with the most significant log2 fold change.

(targetProtein <- inference |> 
   slice(which.min(pval)) |> 
   pull(feature)
)

[1] "P61604"

To obtain the required data, we perform a little data manipulation pipeline:

We use the QFeatures subsetting functionality to retrieve all data related to and focusing on the peptides_log and proteins sets that contains the peptide ion data used for model fitting. We then convert the data with longForm() for plotting. Finally, we plot the log2 normalised intensities for each sample at the protein and at the peptide level. Since multiple peptides are recorded for the protein, we link peptides across samples using a grey line. Samples are colored according to E. coli spike-in condition.

qf[targetProtein, , c("peptides_norm", "proteins")] |> #1
  longForm(colvars = colnames(colData(qf))) |> #2
  data.frame() |> 
  ggplot() +
  aes(x = colname,
      y = value) +
  geom_line(aes(group = rowname), linewidth = 0.1) +
  geom_point(aes(colour = prognosis)) +
  facet_wrap(~ assay, scales = "free") +
  ggtitle(targetProtein) +
  theme_minimal() +
  theme(axis.text.x = element_blank())

qf[targetProtein, , c("peptides_norm", "proteins")] |> #1
  longForm(colvars = colnames(colData(qf))) |> #2
  data.frame() %>% 
  {
  ggplot(.) +
  aes(x = prognosis,
      y = value) +
  geom_boxplot(aes(colour = prognosis)) +
  facet_wrap(~ assay, scales = "free") +
  geom_jitter(aes(shape = rowname)) +
  scale_shape_manual(values = seq_len(dplyr::n_distinct(.$rowname))) +
  ggtitle(targetProtein) +
  theme_minimal() +
  theme(axis.text.x = element_blank()) +
  guides(shape = "none")
  }