Limma, Log2FC and pvalue changing depending on the order of the samples.
1
0
Entering edit mode
eiudkann • 0
@bdf25980
Last seen 16 months ago
United States

I have two sets of samples: ORG and FT. I am using limma do do DE analysis. In the first pass, I put ORG in the second group. In the second pass, I put ORG in the first group. I get a different Log2FC and intercept at the end. What am I doing wrong? It took me many hours to get to this point and searching the site didn't give me answers, so I'd be sincerely grateful if someone could explain this to me.

First pass (ORG is second group):

> intercept <- c(1,1,1,1,1,1)
> ORG.FT.logfc <- c(0,0,0,1,1,1)
> design= <- cbind(intercept, ORG.FT.logfc)
> 
> data.test.contrast <- data %>%
+   select(all_of("Name"), contains("Log.FT_Control"), contains("Log.ORG_Control")) %>%
+   remove_rownames() %>%
+   column_to_rownames(var = "Name")
> data.test.contrast
            Log.FT_Control1 Log.FT_Control2 Log.FT_Control3 Log.ORG_Control1 Log.ORG_Control2 Log.ORG_Control3
ABI1 - P23  -0.057560361      -0.118932196        0.14073431     -0.212971126      -0.34317744    -2.939038e-01
ABI1 - P306 -0.721475735      -0.067965386       -0.72147574     -0.270382213       0.46787382    -2.713435e-01
ABI2 - P393 0.290695541       0.017626361        0.38626237      0.515744621       0.83827463     4.650932e-01

> rownames(design.contrast) <- data.test.contrast %>%
+   select(contains("Log.")) %>%
+   colnames()
> design.contrast
                  intercept ORG.FT.logfc
Log.FT_Control1         1              0
Log.FT_Control2         1              0
Log.FT_Control3         1              0
Log.ORG_Control1          1              1
Log.ORG_Control2          1              1
Log.ORG_Control3          1              1
> 
> fit.mine.type2 <- lmFit(data.test.contrast, design.contrast)
> head(fit.mine.type2$coefficients,3)
           intercept ORG.FT.logfc
ABI1 - P23 -0.01191942   -0.271431375
ABI1 - P306 -0.50363895    0.479021650
ABI2 - P393 0.23152809    0.374842719

> fit.mine.type2eb <- eBayes(fit.mine.type2)
> fit.mine.type2tt <- as.data.frame(topTable(fit.mine.type2eb, coef = 1, number = Inf,sort.by = "none"))

> fit.mine.type2tt
            Log2.FC       AveExpr            t        p.val       adj.p            B
ABI1 - P23 -0.0119194171 -0.1476351043 -0.122179081 9.062931e-01 0.939942639 -7.035021201
ABI1 - P306 -0.5036389522 -0.2641281272 -2.546622916 3.926962e-02 0.111267727 -4.459653621
ABI2 - P393 0.2315280894  0.4189494488  1.896802680 1.009532e-01 0.206808000 -5.402760959

Second pass (ORG first)

> intercept <- c(1,1,1,1,1,1)
> FT.ORG.logfc  <- c(0,0,0,1,1,1)
> design.contrast <- cbind(intercept, FT.ORG.logfc)
> data.test.contrast <- data %>%
+   select(all_of("Name"), contains("Log.ORG_Control"), contains("Log.FT_Control")) %>%
+   remove_rownames() %>%
+   column_to_rownames(var = "Name")
> data.test.contrast
            Log.ORG_Control1 Log.ORG_Control2 Log.ORG_Control3 Log.FT_Control1 Log.FT_Control2 Log.FT_Control3
ABI1 - P23  -0.212971126      -0.34317744    -2.939038e-01      -0.057560361      -0.118932196        0.14073431
ABI1 - P306  -0.270382213       0.46787382    -2.713435e-01      -0.721475735      -0.067965386       -0.72147574
ABI2 - P393 0.515744621       0.83827463     4.650932e-01       0.290695541       0.017626361        0.38626237
> design.contrast
                  intercept FT.ORG.logfc
Log.ORG_Control1          1              0
Log.ORG_Control2          1              0
Log.ORG_Control3          1              0
Log.FT_Control1         1              1
Log.FT_Control2         1              1
Log.FT_Control3         1              1

> fit.mine.type2 <- lmFit(data.test.contrast, design.contrast)
> head(fit.mine.type2$coefficients,10)
              intercept FT.ORG.logfc
ABI1 - P23  -0.28335079    0.271431375
ABI1 - P306  -0.02461730   -0.479021650
ABI2 - P393  0.60637081   -0.374842719

> fit.mine.type2eb <- eBayes(fit.mine.type2)
> fit.mine.type2tt <- as.data.frame(topTable(fit.mine.type2eb, coef = 1, number = Inf,sort.by = "none"))
> colnames(fit.mine.type2tt) <- c("Log2.FC","AveExpr","t", "p.val","adj.p","B")
> fit.mine.type2tt
                   Log2.FC       AveExpr             t        p.val       adj.p          B
ABI1 - P23  -0.283350792 -0.1476351043  -2.904465812 2.361670e-02 0.078743209 -4.0232431
ABI1 - P306   -0.024617302 -0.2641281272  -0.124476047 9.045418e-01 0.938126228 -7.1302850
ABI2 - P393  0.606370808  0.4189494488   4.967715915 1.772071e-03 0.012105942 -1.2516359
limma • 487 views
ADD COMMENT
1
Entering edit mode
@gordon-smyth
Last seen 5 minutes ago
WEHI, Melbourne, Australia

The problem is that your final table doesn't contain log-fold-changes at all, you're just making tables of intercepts. You're getting tables of intercepts because you set coef=1 instead of coef=2. The intercepts change between your first and second passes because you defined a different group to be the reference group. The logFCs will change sign for the same reason.

Apart from the error in specifying the correct coefficient, you seem to be going to unnecessary work. A simpler and correct workflow would be:

design <- model.matrix(~group)
fit <- lmFit(y, design)
fit <- eBayes(fit)
tt <- topTable(fit, n=Inf)
ADD COMMENT

Login before adding your answer.

Traffic: 782 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6