Question: Limma is doing the right way to calculate the fold change?
0
gravatar for Bioinformatics
2.4 years ago by
Bioinformatics30 wrote:

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 

 

 

proteomics limma • 729 views
ADD COMMENTlink modified 2.4 years ago by Laurent Gatto1.2k • written 2.4 years ago by Bioinformatics30
Answer: Limma is doing the right way to calculate the fold change?
1
gravatar for Laurent Gatto
2.4 years ago by
Laurent Gatto1.2k
Belgium
Laurent Gatto1.2k wrote:

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.

ADD COMMENTlink written 2.4 years ago by Laurent Gatto1.2k
1

@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)))
ADD REPLYlink written 2.4 years ago by SamGG190

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

ADD REPLYlink written 2.4 years ago by Bioinformatics30
1

@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

ADD REPLYlink written 2.4 years ago by SamGG190

@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 

ADD REPLYlink written 2.4 years ago by Bioinformatics30

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.

ADD REPLYlink written 2.4 years ago by Laurent Gatto1.2k

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

ADD REPLYlink written 2.4 years ago by Bioinformatics30
1

Yes: the intensity increases upon treatment.

ADD REPLYlink written 2.4 years ago by Laurent Gatto1.2k

Interesting related answer.

HTH

ADD REPLYlink written 2.4 years ago by SamGG190
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 292 users visited in the last hour