nature of the response factor or how to model it correctly

[Tobias Kockmann]

Dear MSstats support,

can somebody tell me if the signal response factor in label-free LC-MS experiments is peptide-specific or precursor-specific? I recently analyzed a lot of DIA data, where in contrast to targeted experiments like SRM/PRM many peptides are covered by multiple precursors having different charge states. I previously had assumed that the response factor  for a peptide i is equal in all charge states z. Accordingly, one would model the linear relationship between LC-MS signal y and peptide quantity x by:

y~a+bx

over all y, regardless of z (z=+2, +3, ...) instead of fitting a separate model for each charge state (one per precursor). But I see examples in my DIA data where a and b are clearly dependent on z.

Is this due to some nuisance factor, so specific to these hand-picked examples, or do I in general need to fit a model for each precursor, or should I even include z as explanatory variable in the linear model and fit a separate offset for a and b per level of z? How does MSstats cope with this? Looking at the code I can not see any explanatory variable in the lm() or lmer() calls, so I guess the model is fit per precursor? What is the rational behind this?

Greetings,

Tobi

[Olga Vitek]

Dear Tobi,

Thank you for your email! This is a very common concern. Several points:

1) MSstats models the intensities on the log scale, and not on the original scale, as in your plots

2) MSstats models the intensities of the individual replicates within each condition, and not summarized per condition, as your plot seems to indicate.

3) To visualize the difference in log-abundances of the precursors between runs, we use 'profile' plots. They help visualize the extent to which the individual precursors vary in a 'parallel' way between the runs.

The deviations from the 'parallel' shape in DIA can be due to many reasons,  e.g. difficulty of computational deconvolution of the precursor or the fragment peaks, but also to ionization competition, and other experimental interferences. MSstats summarizes all the precursor intensities in a single number per protein per run, and does so with a robust procedure that minimizes the effect of outliers, and of analytes below the limit of detection.

By default MSstats groups all the precursors that belong to a protein. However, you can also give a peptide (or evan a precursor!) name in the protein column, and then the summary will be performed at the peptide or the precursor level instead. This is useful when you have reasons to believe that different peptides contain different quantitative information (e.g., due to differential PTMs), and this better answers your biological question. The price for this is that these higher-granularity summaries are based on  less data, are more likely to overfit, and are more likely to produce multiple comparisons artifacts.

I hope that this helps at least somewhat

Sincerely

Olga 

[Tobias Kockmann]

Hi Olga,

first of all thx for your reply. Here is my point-by-point reply:

1) Yes I know that the default is to model the response on the log-transformed scale. If I understood you correctly, this is primarily done to keep the variance constant, since on the raw scale: var(y) = f(y), right? And lm() and lmer() assume constant variance across covariate levels. Anyway, I still have the impression that many peptides show sig. differences in the estimated slope for different charge states (z = 2, 3, ..., n) when doing a regression using log(y). Is that possible from a physical point of view - in the absence of noise?

Greetings,

Tobi

[Tobias Kockmann]

2) The data I plotted is output by Spectronaut 9. They call it the default peptide report. The report looks like this:

> head(PepQuantDef)

  R.Condition          R.FileName R.Replicate     PG.ProteinAccessions              EG.ModifiedSequence

1           1 20160805_17_DIA_1_2           2 A0A087WSV8;V9HW75;E9PKG6             _LSK[+28]ELDLVSHHVR_

2           1 20160805_17_DIA_1_2           2 A0A087WSV8;V9HW75;E9PKG6 _EK[+28]LQK[+28]ADIEEIK[+28]SGR_

3           1 20160805_17_DIA_1_2           2 A0A087WSV8;V9HW75;E9PKG6           _MREHVMNEVDTNK[+28]DR_

4           1 20160805_17_DIA_1_2           2 A0A087WSV8;V9HW75;E9PKG6 _EK[+28]LQK[+28]ADIEEIK[+28]SGR_

5           1 20160805_17_DIA_1_2           2 A0A087WSV8;V9HW75;E9PKG6             _EHVMNEVDTNK[+28]DR_

6           1 20160805_17_DIA_1_2           2 A0A087WSV8;V9HW75;E9PKG6           _MREHVMNEVDTNK[+28]DR_

                      EG.PrecursorId EG.StrippedSequence EG.Qvalue EG.MinProfileQvalue FG.Quantity

1             _LSK[+28]ELDLVSHHVR_.3       LSKELDLVSHHVR         0                   0   113020.56

2 _EK[+28]LQK[+28]ADIEEIK[+28]SGR_.4     EKLQKADIEEIKSGR         0                   0    57729.28

3           _MREHVMNEVDTNK[+28]DR_.4     MREHVMNEVDTNKDR         0                   0    46998.97

4 _EK[+28]LQK[+28]ADIEEIK[+28]SGR_.3     EKLQKADIEEIKSGR         0                   0    13105.55

5             _EHVMNEVDTNK[+28]DR_.3       EHVMNEVDTNKDR         0                   0    10646.84

6           _MREHVMNEVDTNK[+28]DR_.3     MREHVMNEVDTNKDR         0                   0     7002.79

I did not summarise per condition:

> xtabs(~R.Condition+R.Replicate, data = PepQuantDef)

           R.Replicate

R.Condition     1     2     3

          0 27133 27013 26989

          1 27045 26696 26594

          2 26097 27244 26616

          3 25866 25625 26768

          4 26389 27008 27057

          5 26141 26770 25987

The unequal numbers are because low quality information is not reported explicitly by Spectronaut (EG.Q-value filtered). How Spectronaut aggregates from the fragment ion level up to the precursor-level is unknown to me.

3) I prefer to plot y ~ condition in combinations with grouping:

xyplot(FG.Quantity~R.Condition | EG.PrecursorId, ...) instead of having one x value per run. But it should give the same picture. Precursors should behave in a parallel fashion.

[Tobias Kockmann]

You said that deviations from the parallel shape can have many reasons...but you describe all of them as experimental interferences, noise, etc... so factors that smear the underlying true signal. Do I read you correctly: In an ideal experiment, they have to be parallel? What I picked up is due to noise and non-ideal inference methods?

[Olga Vitek]

Dear Tobias,

To some of your points:

* Log transform: The main reason for the log is that it transforms the multiplicative signal on the original scale to the additive signal on the log scale. Additive signals are much easier to model than multiplicative signals, especially in complex designs. In particular, linear (mixed) models all work with additive signals.

* 'Parallel' profiles: there may be many biological reasons for deviations from 'parallel' profiles on the log scale, such as post-translational modifications. Changes in site occupancy between conditions (as opposed to changes in abundance) will definitely result in systematic difference of profiles of various peptides. If you have reasons to believe that this happens in your situation, then it is not appropriate to perform protein-level inference (or, at least, it is not appropriate to include peptides with these modification sites). Protein-level inference focuses on 'overall' changes in protein abundance, and views peptides as 'replicates'. If the peptides are not genuine replicates, then peptide-level inference is preferred. However, as I mentioned previously, peptide-level inference is more likely  to lead to overfitting, and is more susceptible to artifacts of multiple testing.

* Plotting 'condition' vs 'run' on the x axis: we prefer plotting 'run' on x axis, as it allows us to distinguish the extent of random interferences within conditions and between conditions for each peptide ion. Plotting 'condition' as x axis confounds these two types of variation.

* Constant variance: Since MSstats 3.0, the 'variance' refers to the biological variation of the *entire protein* between runs, after summarization. Generally the log transformation does not ensure that this variance is constant, in particular it can depends on the protein abundance. Linear models can handle this type of unequal variance.

HTH

Olga

[Tobias Kockmann]

Dear Olga,

first of all a happy new year and thx for your answer. I am indeed interested in performing peptide-level inference, since my experimental conditions affect protein processing, so it is a kind of PTM analysis. But I really would like to make sure we are on the same page: Different precursors of the same peptide - in the absence of noise or factors like unexpected modifications - should behave alike in a linear regression? Thought experiment: I would buy 10 synthetic peptides, mix them, pipette a linear dilution series and measure all available precursors by DIA. Will the slope be the same for peptide i in z = 2, 3 , 4, ... or can peptide i | z=2 be different from peptide i | z=3 given the signals allow a high quality regression for both charges. Does the response factor depend on z? 

Greetings,

Tobi

[Olga Vitek]

Dear Tobias,

I am not sure which linear regression you have in mind. See, e.g. the figure on p.15 of the user manual http://msstats.org/wp-content/uploads/2016/11/MSstats_v3.6_manual.pdf

All the features that represent a protein (in your case, a peptide) are expected to have same *systematic* patterns of variation between the runs. I.e., they are expected to systematically go up and down between the runs in the same manner, as in the left part of the figure.  All the deviations from this systematic up and down variation is assumed to be non-systematic, irreproducible, noise. The summarization step finds a consensus pattern of the systematic variation, while accounting for censored values, missing values, and outliers, as shown in the right part of the figure.

HTH

Olga