Limma is doing the right way to calculate the fold change?
1
0
Entering edit mode
@bioinformatics-10931
Last seen 2.4 years ago
United States

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

limma proteomics • 2.2k views
1
Entering edit mode
@laurent-gatto-5645
Last seen 6 days ago
Belgium

From the data you show at the top, all proteins are down-regulated because the intensities in your controls are lower than in your treatment, hence negative log2 fold-changes.

If it is the negative values that confuse you, it is because down-regulation (i.e fold-changes between 0 and 1) become negative in log space (-Inf to 0), while up-regulation (i.e. fold-changes between 1 and +Inf) are positive in log space (0 to +Inf).

Hope this helps.

1
Entering edit mode

@Bioinformatics did you try to swap your design? May be it's what you are looking for ;-)

model.matrix(~c(rep(0,2),rep(1,2)))
0
Entering edit mode

@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?

1
Entering edit mode

@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-rA: 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

0
Entering edit mode

@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

0
Entering edit mode

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.

0
Entering edit mode

@Laurent Gatto so it means that they are upregulated ? right?

1
Entering edit mode

Yes: the intensity increases upon treatment.

0
Entering edit mode

HTH