Dear Francesco, I think that you may be mis-judging the nature of a linear model with multiple factors and covariates. Limma computes log-fold-changes by fitting a linear model and extracting the estimated coefficients. The interpretation of each coefficient is that it represents the change in the reponse (log-expression in this case) for each unit change in the covariate, assuming that all other covariates in the model can be held constant. For example, the coefficient for alpha in your hypothetical model is the change in the log-expression of the gene that would occur if alpha was changed (from 0 to 1) but all other factors in the model were held constant. This is not the same as fitting a simple model with just alpha and no other factors. The coefficient that you would get from the simple model with just alpha will not be the same as the coefficient you would get for alpha from the full model. In the simple model, the other factors are not held constant. It would easily be that the full model could give a negative coefficient but the simple model gives a positive coefficient. This occurs when there are other factors in the model that are positively correlated with alpha but have the opposite effect on expression. When you make a simple summary of log-expression values for wild-type KRAS vs. mutated, you are looking at the marginal effect of KRAS when all other factors are allowed to vary as well. In other words, you are examining the simple model with just KRAS as a factor. There is no reason that this should be the same as the logFC for KRAS when all other factors have been adjusted for. The simple effect could easily be driven by factors other than KRAS. You might find it helpful to consult a textbook on multiple regression that discusses the subtleties of unbalanced anova or correlated covariates. Best wishes Gordon > Date: Mon, 14 Apr 2014 06:22:45 -0700 (PDT) > From: "Francesco Gatto [guest]" <guest at="" bioconductor.org=""> > To: bioconductor at r-project.org, gatto at chalmers.se > Subject: [BioC] Interpretation of fold.changes for covariates in a > linear model using voom > > I have implemented a linear model to fit RNAseq read counts for ~1000 > samples using ~200 covariates plus the intercept. The pipeline follows > the example in the user manual, i.e. v=voom(y,design) -> > fit=lmFit(v,design) -> fit2=eBayes(fit). Then I am using decideTests() > to assess the effect of each covariate on gene expression. However, the > output is less complete than topTable(), so I used this function to > obtain values on fold-changes and individual p-values. At this stage, > values become difficult to interpret. If the LM is: > > gene_i,j = intercept_i + alfa_i*Factor1_j + beta_i*Factor2_j + (~200 other factors) > > and, to simplify, let's say covariates are binary. Then what do the > fold-changes for topTable(fit2,coef="alfa") mean? > My interpretation was the up/down-regulation due to the presence of > covariate alfa, but if I boxplot or summarize the expression values for > the two groups (alfa = 1 vs alfa = 0) I get quite contradictory results. > For example this is the most DE gene when the sample has a mutation in > KRAS (the covariate, 0 if wild-type, 1 if mutated): > >> topTable(fit2,coef="metadataKRAS",number=1) > > hgnc_symbol logFC AveExpr t P.Value adj.P.Val B > GOLGA6A -1.314213 -5.9000494 -7.375175 3.800569e-13 7.204358e-09 18.80961 > > The summary of gene expression values for wild-type KRAS vs. mutated is: > >> by(E["GOLGA6A",],design[,"metadataKRAS"],summary) > > design[, "metadataKRAS"]: 0 > Min. 1st Qu. Median Mean 3rd Qu. Max. > -7.882 -6.865 -6.461 -6.021 -5.770 1.548 > ---------------------------------------------------------------- > design[, "metadataKRAS"]: 1 > Min. 1st Qu. Median Mean 3rd Qu. Max. > -7.657 -6.247 -5.487 -5.239 -4.582 -1.311 > > This is an up-regulation, and of magnitude way lower than the one > indicated by topTable. Is this because I am misusing the function > topTable? I cannot wrap my head around this. Thanks for your help! > > -- output of sessionInfo(): > > R version 3.0.2 (2013-09-25) > Platform: x86_64-apple-darwin10.8.0 (64-bit) > > locale: > [1] sv_SE.UTF-8/sv_SE.UTF-8/sv_SE.UTF-8/C/sv_SE.UTF-8/sv_SE.UTF-8 > > attached base packages: > [1] stats graphics grDevices utils datasets methods base > > other attached packages: > [1] edgeR_3.4.2 limma_3.18.9 > > loaded via a namespace (and not attached): > [1] tools_3.0.2 ______________________________________________________________________ The information in this email is confidential and intend...{{dropped:4}}