Actually i am trying to find upregulated proteins using limma package which is very straight forward
The problem is that I cannot find a logic behind of it. As an example
This is the data obtained from Mass spectrometery
|
Control_1 |
Control_2 |
Treated_1 |
Treated_2 |
P1 |
27.2 |
34.9 |
1871 |
1710.7 |
P2 |
62.2 |
74.5 |
1258.1 |
1485.8 |
P3 |
4.5 |
3.3 |
437.6 |
207.7 |
P4 |
33.6 |
58.3 |
2854.6 |
2580.8 |
P5 |
1202.9 |
795.3 |
922.7 |
1489.1 |
P6 |
1012.3 |
1307.3 |
1545.4 |
2236 |
Q7 |
36 |
46.1 |
1564.9 |
1633.2 |
P8 |
4.6 |
8.4 |
174.3 |
206.1 |
In R format
df<- structure(list(X = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 8L, 7L ), .Label = c("P1", "P2", "P3", "P4", "P5", "P6", "P8", "Q7"), class = "factor"), Control_1 = c(27.2, 62.2, 4.5, 33.6, 1202.9, 1012.3, 36, 4.6), Control_2 = c(34.9, 74.5, 3.3, 58.3, 795.3, 1307.3, 46.1, 8.4), Treated_1 = c(1871, 1258.1, 437.6, 2854.6, 922.7, 1545.4, 1564.9, 174.3), Treated_2 = c(1710.7, 1485.8, 207.7, 2580.8, 1489.1, 2236, 1633.2, 206.1)), .Names = c("X", "Control_1", "Control_2", "Treated_1", "Treated_2"), class = "data.frame", row.names = c(NA, -8L))
When I do the Limma, with a design like below
design <- model.matrix(~c(rep(1,2),rep(0,2)))
I get the following
|
Fold Change |
AveExper |
t |
p_value |
adj.p_value |
B |
P1 |
-1759.8 |
910.95 |
-25.5988 |
2.47E-04 |
0.005736 |
-4.13352 |
P2 |
-1303.6 |
720.15 |
-13.3532 |
0.001444 |
0.012069 |
-4.14176 |
P3 |
-318.75 |
163.275 |
-3.23851 |
0.05468 |
0.10845 |
-4.28044 |
P4 |
-2671.75 |
1381.825 |
-22.70456 |
0.000342 |
0.005982 |
-4.134362 |
P5 |
-206.8 |
1102.5 |
-0.692565 |
0.542838 |
0.57722 |
-4.683112 |
P6 |
-730.9 |
1525.25 |
-2.274524 |
0.116155 |
0.183387 |
-4.374939 |
Q7 |
-1558 |
820.05 |
-52.42197 |
3.50E-05 |
0.002334 |
-4.131131 |
P8 |
-183.7 |
98.35 |
-13.04355 |
0.001538 |
0.012383 |
-4.142291 |
in R format df2<- structure(list(X = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 8L, 7L ), .Label = c("P1", "P2", "P3", "P4", "P5", "P6", "P8", "Q7"), class = "factor"), Fold.Change = c(-1759.8, -1303.6, -318.75, -2671.75, -206.8, -730.9, -1558, -183.7), AveExper = c(910.95, 720.15, 163.275, 1381.825, 1102.5, 1525.25, 820.05, 98.35), t = c(-25.59881, -13.353164, -3.238507, -22.70456, -0.6925649, -2.274524, -52.42197, -13.04355), p_value = c(0.000247059, 0.001443824, 0.05467982, 0.000342436, 0.5428383, 0.1161553, 3.50056e-05, 0.001538154), adj.p_value = c(0.005736091, 0.012069176, 0.1084498, 0.005981734, 0.5772197, 0.183387, 0.002334204, 0.01238259 ), B = c(-4.133516, -4.141755, -4.280436, -4.134362, -4.683112, -4.374939, -4.131131, -4.142291)), .Names = c("X", "Fold.Change", "AveExper", "t", "p_value", "adj.p_value", "B"), class = "data.frame", row.names = c(NA, -8L))
I know that I should log2 the value at the beginning but it is strange to me that a fold change be minus so it will all be down regulated proteins ?
Any comment is appreciated
@Bioinformatics did you try to swap your design? May be it's what you are looking for ;-)
@SamGG can you please explain it to me, why you think I should swap the design? I think the fold change should be treated/control , right?
@Bioinformatics Although limma sounds easy and convenient, you should encode your design carefully. What I like is to use named values as "Ctrl" and "Treat" instead of 1 and 0, and to define the difference I want to look at using makeContrasts() instead of relying on the default contrast convention.
Read the doc of limma of course, but also https://stats.stackexchange.com/questions/64249/creating-contrast-matrix-for-linear-regression-in-r, A: Not able to reproduce limma results from affy chip with or without CEL file from and other posts concerning contrast coding with limma. NB: the name of the intercept should always be "(Intercept)' as limma search and removes it.
HTH
@Laurent Gatto thanks for your comment which I appreciate it a lot. I thought the up regulation means when the abundance in the treated is significantly increased in treated sample compare to that of control! do I miss something? if yes, how do you define up-regulation is such a case ? Thanks again
See @SamGG 's comment. Also, the log2 fold-changes are symetrical, so -1759.8 when comparing control vs. treatment becomes 1759.8 when comparing treatment vs control.
@Laurent Gatto so it means that they are upregulated ? right?
Yes: the intensity increases upon treatment.
Interesting related answer.
HTH