Question: limma fold change problem
0
gravatar for Learner
3.1 years ago by
Learner 180
Learner 180 wrote:

Hello,

I know this has been asked many times but I did not figure it out programmatically. So I tried to make a small reproducible example to see if someone knows about it

Here is a data

df <- structure(list(control1 = c(28.2494785829057, 28.1476353572629, 
27.1929113349158, 25.6702823986081, 26.8927288390394, 0, 26.7207761577158, 
26.9122794096471, 31.7603001024038, 29.6522061602385, 28.6278754368064, 
28.6577526656622, 24.9082071983177, 26.7735417008271, 27.0352760202938
), control2 = c(28.2584923577297, 28.1021446453534, 26.2233747158687, 
26.0274297329297, 26.9801099459161, 25.358504222508, 26.7976110787887, 
27.2493429511711, 31.7579581611434, 29.7188834878171, 28.5234274780504, 
28.6846850151125, 0, 27.0461438308603, 27.2389971036243), treatment1 = c(28.8868697499915, 
28.3221674791755, 26.1941520812899, 25.1512405315526, 26.8643148321299, 
25.2739455625142, 26.3634923544076, 27.3338132331068, 31.4077578366191, 
29.6662104574918, 28.3396435112577, 28.5763261654449, 24.9068319368666, 
26.8338833938774, 27.3473948352458), treatment2 = c(28.6182441737351, 
28.3364546698625, 26.5573504729984, 25.4954893749752, 26.9793468747426, 
25.7376090355345, 26.6356105379829, 27.4102977563763, 31.5699592291352, 
29.8796437292758, 28.1989170942089, 28.5062377808235, 24.6275067774911, 
26.8883222436102, 26.8426609128988)), .Names = c("control1", 
"control2", "treatment1", "treatment2"), row.names = c("guk-1", 
"dab-1", "swsn-9", "row-1", "ZC434.7", "rad-50", "sin-3", "hmp-1", 
"dhs-3", "nra-2", "saps-1", "D1081.7", "aph-2", "F15B9.10", "C27A7.5"
), class = "data.frame")

This data has 2 control and 2 treatment conditions. I want to get fold change and p-values based on Limma This is continuous data which is suitable for limma, so please don't ask redundant question like what type of data etc .

Here is what I do

design <- model.matrix(~c(rep(0,2),rep(1,2)))
fit <- lmFit(df, design)
fit2 <- eBayes(fit)
t <- topTable(fit2, coef=2, n=Inf)

which gives me this

               logFC  AveExpr          t    P.Value adj.P.Val         B
dab-1     0.20442107 28.22710  3.9514485 0.03685753 0.2271467 -3.905281
guk-1     0.49857149 28.50327  3.9394325 0.03711921 0.2271467 -3.906530
saps-1   -0.30637115 28.42247 -3.4043286 0.05179096 0.2271467 -3.971773
dhs-3    -0.27027060 31.62399 -3.1719049 0.06057245 0.2271467 -4.007078
row-1    -0.52549111 25.58611 -2.3700569 0.11096701 0.3051381 -4.174433
D1081.7  -0.12993687 28.60625 -2.2558390 0.12205525 0.3051381 -4.205428
sin-3    -0.25964217 26.62937 -1.9634112 0.15745464 0.3298828 -4.294491
hmp-1     0.29124431 27.22643  1.8412188 0.17593747 0.3298828 -4.336009
rad-50   12.82652519 19.09251  1.1581479 0.34133736 0.5256822 -4.609835
aph-2    12.31306576 18.61064  1.1320028 0.35045482 0.5256822 -4.621237
nra-2     0.08738227 29.72924  0.8052242 0.48723341 0.6555833 -4.761318
swsn-9   -0.33239175 26.54195 -0.7312325 0.52446668 0.6555833 -4.791007
F15B9.10 -0.04873995 26.88547 -0.3740694 0.73651132 0.8498207 -4.908211
ZC434.7  -0.01458854 26.92913 -0.1849504 0.86662644 0.8745856 -4.943976
C27A7.5  -0.04210869 27.11608 -0.1737547 0.87458562 0.8745856 -4.945376

if I look at this post, R: How to convert log2FC (Fold Change) obtained by limma's topTable() function to FC I should be able to get the answer by one of the following example but I don't

example1 <- log2(rowMeans(df[,c(3,4)])/(rowMeans(df[,c(1,2)])))

example2 <- log2(rowMeans(df[,c(3,4)])-(rowMeans(df[,c(1,2)])))

example3 <- log(rowMeans(df[,c(3,4)])/(rowMeans(df[,c(1,2)])))

So I have two question,

1- the way I am using the limma is right? (especially the design I make).

2- how really the fold change is calculated. Please avoid words , just show it by R

fold change limma R • 1.3k views
ADD COMMENTlink modified 3.1 years ago by russhh4.7k • written 3.1 years ago by Learner 180
2
gravatar for russhh
3.1 years ago by
russhh4.7k
UK, U. Glasgow
russhh4.7k wrote:

limma will assume your data is already log2 transformed. So, re 2:

example4 <- rowMeans(df[, 3:4]) - rowMeans(df[, 1:2])

gives the log2-fold change. If your data has not been transformed into log-space, you should do this yourself prior to calling limma.

The way you are using limma is ok, it's just a bit ugly: df and t are internal R functions, you shouldn't really be overwriting their defintiions; and you model.matrix will be a bit more self-explanatory if you used a factor to distinguish treatment from control, ie:

treat <- factor(rep(c('control', 'treatment'), each = 2), levels = c('control', 'treatment'))
design <- model.matrix(~ treat)

all the rest is fine

ADD COMMENTlink written 3.1 years ago by russhh4.7k

@russhh thanks for the nice and clear answer right to the point which i liked don't you think that we should get treatment/control for fold change ? also do you know how they calculate p-value in limma?

ADD REPLYlink written 3.1 years ago by Learner 180

re treatment / control values; these aren't provided by default but then, they're superfluous: just exponentiate the log2 FC: 2 ^ logFC

... It's kind of assumed that you have street-fighting levels of maths.

The p-values are a different kettle of fish. You'd be better looking into how p-values are calculated in basic t-tests / ANOVA, before learning how the limma/eBayes-regularised p-values (and subsequently the multiple-comparison adjusted p-values) are determined.

ADD REPLYlink written 3.1 years ago by russhh4.7k
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 2.3.0
Traffic: 1456 users visited in the last hour