subject nesting
89 views
Skip to first unread message
Tobias Kockmann
Jul 25, 2016, 7:08:50 AM7/25/16
to MSstats
Dear MSstats support,
I have a question regarding subject nesting in MSstats: My data contains two groups of subjects. (A) patient group (B) control group. From each individual in group A two biopsies were taken, (i) from a diseased tissue (ii) from a matched normal tissue. For individuals from group B a single biopsy was taken. I analysed the data in two ways: I started by ignoring the individual level and simply used the run name as BioReplicate id. This should lead to a statistical analysis treating each run as independent sample.
> str(input.df)
'data.frame': 1715433 obs. of 10 variables:
$ ProteinName : Factor w/ 2516 levels "1433B_HUMAN",..: 2395 2395 2395 2395 2395 2395 2395 2395 2395 2395 ...
$ PeptideSequence : Factor w/ 19585 levels "_[-131]MAALRPLVKPK_",..: 13085 13085 13085 13085 13085 13085 12547 12547 12547 12547 ...
$ PrecursorCharge : Factor w/ 8 levels "2","3","4","5",..: 1 1 1 1 1 1 1 1 1 1 ...
$ FragmentIon : Factor w/ 433 levels "a10H2O","a10H3PO4",..: 433 287 255 153 429 291 417 421 413 399 ...
$ ProductCharge : int 1 1 1 1 1 1 1 1 1 1 ...
$ IsotopeLabelType: Factor w/ 1 level "L": 1 1 1 1 1 1 1 1 1 1 ...
$ Condition : Factor w/ 3 levels "ko","l","nl": 1 1 1 1 1 1 1 1 1 1 ...
$ BioReplicate : Factor w/ 15 levels "20160609_06_classicDIA_37757",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Run : Factor w/ 15 levels "20160609_06_classicDIA_37757",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Intensity : num 6001 3629 24347 742 12315 ...
> table(input.df$BioReplicate, input.df$Condition)
ko l nl
20160609_06_classicDIA_37757 112436 0 0
20160609_07_classicDIA_37750 0 0 112317
20160609_08_classicDIA_37755 0 107410 0
20160609_09_classicDIA_37753 0 114483 0
20160609_10_classicDIA_37749 0 112705 0
20160609_11_classicDIA_37756 0 0 113566
20160609_12_classicDIA_37762 110959 0 0
20160609_14_classicDIA_37748 0 0 117743
20160609_15_classicDIA_37751 0 119712 0
20160609_16_classicDIA_37747 0 115626 0
20160609_17_classicDIA_37766 114848 0 0
20160609_18_classicDIA_37752 0 0 115898
20160609_19_classicDIA_37754 0 0 119652
20160609_20_classicDIA_37759 117787 0 0
20160609_22_classicDIA_37763 110291 0 0
> table(input.df$Run, input.df$Condition)
ko l nl
20160609_06_classicDIA_37757 112436 0 0
20160609_07_classicDIA_37750 0 0 112317
20160609_08_classicDIA_37755 0 107410 0
20160609_09_classicDIA_37753 0 114483 0
20160609_10_classicDIA_37749 0 112705 0
20160609_11_classicDIA_37756 0 0 113566
20160609_12_classicDIA_37762 110959 0 0
20160609_14_classicDIA_37748 0 0 117743
20160609_15_classicDIA_37751 0 119712 0
20160609_16_classicDIA_37747 0 115626 0
20160609_17_classicDIA_37766 114848 0 0
20160609_18_classicDIA_37752 0 0 115898
20160609_19_classicDIA_37754 0 0 119652
20160609_20_classicDIA_37759 117787 0 0
20160609_22_classicDIA_37763 110291 0 0
Next, I used subject ids to populate the BioReplicate column in the input df. This lead to the following:
> str(input.df)
'data.frame': 1715433 obs. of 10 variables:
$ ProteinName : Factor w/ 2516 levels "1433B_HUMAN",..: 2395 2395 2395 2395 2395 2395 2395 2395 2395 2395 ...
$ PeptideSequence : Factor w/ 19585 levels "_[-131]MAALRPLVKPK_",..: 13085 13085 13085 13085 13085 13085 12547 12547 12547 12547 ...
$ PrecursorCharge : Factor w/ 8 levels "2","3","4","5",..: 1 1 1 1 1 1 1 1 1 1 ...
$ FragmentIon : Factor w/ 433 levels "a10H2O","a10H3PO4",..: 433 287 255 153 429 291 417 421 413 399 ...
$ ProductCharge : int 1 1 1 1 1 1 1 1 1 1 ...
$ IsotopeLabelType: Factor w/ 1 level "L": 1 1 1 1 1 1 1 1 1 1 ...
$ Condition : Factor w/ 3 levels "ko","l","nl": 1 1 1 1 1 1 1 1 1 1 ...
$ BioReplicate : Factor w/ 10 levels "ZETH_21","ZETH_22",..: 8 8 8 8 8 8 8 8 8 8 ...
$ Run : Factor w/ 15 levels "20160609_06_classicDIA_37757",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Intensity : num 6001 3629 24347 742 12315 ...
> table(input.df$BioReplicate, input.df$Condition)
ko l nl
ZETH_21 0 112705 112317
ZETH_22 0 119712 115898
ZETH_23 0 114483 119652
ZETH_24 110959 0 0
ZETH_25 117787 0 0
ZETH_26 0 115626 117743
ZETH_27 0 107410 113566
ZETH_28 112436 0 0
ZETH_29 114848 0 0
ZETH_30 110291 0 0
> table(input.df$Run, input.df$Condition)
ko l nl
20160609_06_classicDIA_37757 112436 0 0
20160609_07_classicDIA_37750 0 0 112317
20160609_08_classicDIA_37755 0 107410 0
20160609_09_classicDIA_37753 0 114483 0
20160609_10_classicDIA_37749 0 112705 0
20160609_11_classicDIA_37756 0 0 113566
20160609_12_classicDIA_37762 110959 0 0
20160609_14_classicDIA_37748 0 0 117743
20160609_15_classicDIA_37751 0 119712 0
20160609_16_classicDIA_37747 0 115626 0
20160609_17_classicDIA_37766 114848 0 0
20160609_18_classicDIA_37752 0 0 115898
20160609_19_classicDIA_37754 0 0 119652
20160609_20_classicDIA_37759 117787 0 0
20160609_22_classicDIA_37763 110291 0 0
One can see that individual 21, 22, 23, 26, 27 belong to the patient group, because they provide features coming from condition l and nl. Next, I generated the following comparisons:
> comparisions
ko l nl
l vs ko -1 1 0
nl vs ko -1 0 1
l vs nl 0 1 -1
and compared the results for nested vs. non nested input data. My expectation was that only the comparison l vs nl is affected, but I see that all three comparisons are effected. Can anyone explain me why?
Greetings,
Tobi
Tobias Kockmann
Jul 25, 2016, 10:38:20 AM7/25/16
to MSstats
Here is an example for a protein that is called sig. starting from the non nested data vs. the subject nested data:
> joined.results[index,]
Protein Label log2FC SE Tvalue DF pvalue adj.pvalue issue MissingPercentage ImputationPercentage
2653 H1X_HUMAN l vs ko -4.630406 1.215204 -3.810393 12 0.002482452 0.01717452 <NA> 0.4 0.4
26531 H1X_HUMAN l vs ko -4.630406 1.205021 -3.842593 4 0.018418697 0.12807816 <NA> 0.4 0.4
nesting
2653 run
26531 subject
The FC is the same, but the sig. is completely different due to the DF used to compute the p-value. For the non nested data I am estimating 3 group mean, so 15-3=12, but why do I have 4 DF in the nested analysis?
Meena Choi
Jul 26, 2016, 12:10:21 AM7/26/16
to MSstats
For the second case, it is considered as paired design with 10 different subject and MSstats used the linear model, LogIntensity ~ Group + Subject. Then, degree of freedom 14-2-8=4.
Based on your explanation about your experiment, it should be considered as paired design.
HTH.
- Meena
Tobias Kockmann
Jul 26, 2016, 5:30:07 AM7/26/16
to MSstats
Hi Meena,
I feel we are at the heart of the problem. So you are saying the paired design is modelled by introducing a fixed subject effect into the linear model? I tried to reproduce this using the following R code:
> Runlevel.intensities <- droplevels(QuantData.top10.subjectNesting$RunlevelData[index,])
> str(Runlevel.intensities)
'data.frame': 15 obs. of 12 variables:
$ RUN : Factor w/ 15 levels "1","10","11",..: 1 8 9 10 11 12 13 14 15 2 ...
$ Protein : Factor w/ 1 level "H1X_HUMAN": 1 1 1 1 1 1 1 1 1 1 ...
$ LogIntensities : num 14.08 13.74 13.62 9.83 14.09 ...
$ NumMeasuredFeature: int 10 8 10 4 8 4 4 4 1 7 ...
$ MissingPercentage : num 0 0.2 0 0.6 0.2 0.6 0.6 0.6 0.9 0.3 ...
$ more50missing : Named logi FALSE FALSE FALSE TRUE FALSE TRUE ...
..- attr(*, "names")= chr "1" "2" "3" "4" ...
$ NumImputedFeature : int 0 2 0 6 2 6 6 6 9 3 ...
$ GROUP : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 2 2 2 2 2 ...
$ GROUP_ORIGINAL : Factor w/ 3 levels "ko","l","nl": 1 1 1 1 1 2 2 2 2 2 ...
$ SUBJECT_ORIGINAL : Factor w/ 10 levels "ZETH_28","ZETH_21",..: 1 5 8 9 10 2 3 4 6 7 ...
$ SUBJECT_NESTED : Factor w/ 15 levels "1.1","1.5","1.8",..: 1 2 3 4 5 6 7 8 9 10 ...
$ SUBJECT : Factor w/ 10 levels "1","2","3","4",..: 1 5 8 9 10 2 3 4 6 7 ...
> lm2 <- lm(LogIntensities ~ GROUP_ORIGINAL + SUBJECT_ORIGINAL, data = Runlevel.intensities)
> summary(lm2)
Call:
lm(formula = LogIntensities ~ GROUP_ORIGINAL + SUBJECT_ORIGINAL,
data = Runlevel.intensities)
Residuals:
11999 12000 12001 12002 12003 12004 12005 12006 12007 12008 12009
-1.110e-16 2.752e-17 -3.029e-17 4.181e-17 5.742e-17 -1.382e-01 9.343e-01 -8.800e-01 -1.121e-02 9.514e-02 1.382e-01
12010 12011 12012 12013
-9.343e-01 8.800e-01 1.121e-02 -9.514e-02
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.0754 0.9153 15.378 0.000104 ***
GROUP_ORIGINALl -3.5699 1.1578 -3.083 0.036809 *
GROUP_ORIGINALnl -0.1639 1.1578 -0.142 0.894287
SUBJECT_ORIGINALZETH_21 -1.2735 0.9153 -1.391 0.236520
SUBJECT_ORIGINALZETH_27 -3.2907 0.9153 -3.595 0.022856 *
SUBJECT_ORIGINALZETH_23 -1.1236 0.9153 -1.228 0.286930
SUBJECT_ORIGINALZETH_24 -0.3353 1.2944 -0.259 0.808395
SUBJECT_ORIGINALZETH_26 -4.6256 0.9153 -5.054 0.007213 **
SUBJECT_ORIGINALZETH_22 NA NA NA NA
SUBJECT_ORIGINALZETH_29 -0.4538 1.2944 -0.351 0.743566
SUBJECT_ORIGINALZETH_25 -4.2405 1.2944 -3.276 0.030618 *
SUBJECT_ORIGINALZETH_30 0.0187 1.2944 0.014 0.989167
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.9153 on 4 degrees of freedom
Multiple R-squared: 0.9671, Adjusted R-squared: 0.8849
F-statistic: 11.76 on 10 and 4 DF, p-value: 0.01486
So, as you explained the modelling of the subject effect consumes a lot of DFs. What I am asking myself right now: Does it really make sense to model it this way? I am interested in comparing l vs. ko, nl vs. ko and l vs. nl, for 2 out of 3 comparison (l vs. ko and nl vs. ko) there is no overlap in subjects. I am comparing 5 samples originating from 5 different individuals to 5 samples coming from 5 other individuals. The subject groups are completely disjoint for these comparisons. So why should I start modelling subject effects?
What about the simpler model ignoring subject labels?
> lm1 <- lm(LogIntensities ~ GROUP_ORIGINAL, data = Runlevel.intensities)
> summary(lm1)
Call:
lm(formula = LogIntensities ~ GROUP_ORIGINAL, data = Runlevel.intensities)
Residuals:
Min 1Q Median 3Q Max
-3.238 -1.228 0.651 1.012 2.158
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.0732 0.8593 15.214 3.3e-09 ***
GROUP_ORIGINALl -4.6304 1.2152 -3.810 0.00248 **
GROUP_ORIGINALnl -1.2244 1.2152 -1.008 0.33357
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.921 on 12 degrees of freedom
Multiple R-squared: 0.5651, Adjusted R-squared: 0.4926
F-statistic: 7.797 on 2 and 12 DF, p-value: 0.006765
> BIC(lm1,lm2)
df BIC
lm1 4 69.64494
lm2 12 52.58356
> AIC(lm1, lm2)
df AIC
lm1 4 66.81273
lm2 12 44.08695
Wouldn't this suggest that the simpler model is the better choice to explain my data?
Greetings,
Tobi
Olga Vitek
Jul 26, 2016, 11:40:40 AM7/26/16
to MSstats
Dear Tobias,
Subject should be a random effect, and your design should be expressed with lmer
HTH
Olga
Tobias Kockmann
Jul 26, 2016, 5:13:04 PM7/26/16
to MSstats
Like this (random intercept model)?
> lme1 <- lmer(LogIntensities~GROUP_ORIGINAL+(1|SUBJECT_ORIGINAL), data = Runlevel.intensities)
> summary(lme1)
Linear mixed model fit by REML ['lmerMod']
Formula: LogIntensities ~ GROUP_ORIGINAL + (1 | SUBJECT_ORIGINAL)
Data: Runlevel.intensities
REML criterion at convergence: 50.7
Scaled residuals:
Min 1Q Median 3Q Max
-1.1997 -0.3630 0.1673 0.2598 1.0992
Random effects:
Groups Name Variance Std.Dev.
SUBJECT_ORIGINAL (Intercept) 2.8005 1.6735
Residual 0.8297 0.9109
Number of obs: 15, groups: SUBJECT_ORIGINAL, 10
Fixed effects:
Estimate Std. Error t value
(Intercept) 13.0732 0.8521 15.343
GROUP_ORIGINALl -4.6304 1.2050 -3.843
GROUP_ORIGINALnl -1.2244 1.2050 -1.016
Correlation of Fixed Effects:
(Intr) GROUP_ORIGINALl
GROUP_ORIGINALl -0.707
GROUP_ORIGINALn -0.707 0.886
> coef(lme1)
$SUBJECT_ORIGINAL
(Intercept) GROUP_ORIGINALl GROUP_ORIGINALnl
ZETH_28 13.84638 -4.630406 -1.224353
ZETH_21 13.76061 -4.630406 -1.224353
ZETH_27 12.00363 -4.630406 -1.224353
ZETH_23 13.89119 -4.630406 -1.224353
ZETH_24 13.58771 -4.630406 -1.224353
ZETH_26 10.84098 -4.630406 -1.224353
ZETH_22 14.86978 -4.630406 -1.224353
ZETH_29 13.49626 -4.630406 -1.224353
ZETH_25 10.57507 -4.630406 -1.224353
ZETH_30 13.86080 -4.630406 -1.224353
attr(,"class")
[1] "coef.mer"
How does MSstats calc the p-values for the fixed effects?
Tobias Kockmann
Jul 27, 2016, 11:25:11 AM7/27/16
to MSstats
I also tried it with the alternative LME package:
> lme2 <- lme(LogIntensities~GROUP_ORIGINAL, random = ~1|SUBJECT_ORIGINAL, data = Runlevel.intensities)
> summary(lme2)
Linear mixed-effects model fit by REML
Data: Runlevel.intensities
AIC BIC logLik
60.73744 63.16197 -25.36872
Random effects:
Formula: ~1 | SUBJECT_ORIGINAL
(Intercept) Residual
StdDev: 1.673473 0.9108658
Fixed effects: LogIntensities ~ GROUP_ORIGINAL
Value Std.Error DF t-value p-value
(Intercept) 13.073241 0.8520786 9 15.342763 0.0000
GROUP_ORIGINALl -4.630406 1.2050211 3 -3.842593 0.0311
GROUP_ORIGINALnl -1.224353 1.2050211 3 -1.016043 0.3844
Correlation:
(Intr) GROUP_ORIGINALl
GROUP_ORIGINALl -0.707
GROUP_ORIGINALnl -0.707 0.886
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.1996509 -0.3630210 0.1673296 0.2598298 1.0991450
Number of Observations: 15
Number of Groups: 10
> coef(lme2)
(Intercept) GROUP_ORIGINALl GROUP_ORIGINALnl
ZETH_28 13.84637 -4.630406 -1.224353
ZETH_21 13.76061 -4.630406 -1.224353
ZETH_27 12.00363 -4.630406 -1.224353
ZETH_23 13.89119 -4.630406 -1.224353
ZETH_24 13.58771 -4.630406 -1.224353
ZETH_26 10.84098 -4.630406 -1.224353
ZETH_22 14.86978 -4.630406 -1.224353
ZETH_29 13.49626 -4.630406 -1.224353
ZETH_25 10.57507 -4.630406 -1.224353
ZETH_30 13.86080 -4.630406 -1.224353
Here I get a p-value.
Reply all
Reply to author
Forward
0 new messages