The results (logFC) are not consistence (edgeR packages)

Question

[solarchan7]

I'm trying to find the differences between the test and control groups (genes). But the results are not consistence. For example, when there's a decrease in the test group, the fold change is positive! This happens in about 5/200 of them. Here's my codes, I don't know what's wrong... (because of the signed contact, I can't post any data here, sorry!) :

The data (df) is like:

gene Test1 Test2 Test3 Control1 Control2 Control3 Control4

A ...

B ...

C ...
...
n = 200

groups = factor(c(rep("test",3), rep("control",4)))
dge = DGEList(as.matrix(df), group = groups) 
dge = calcNormFactors(dge, method = "TMM") 
design_matrix = model.matrix(~0+groups, data = dge$samples) 
colnames(design_matrix) = c("test", "control") 
dge = estimateDisp(dge, design_matrix, robust=TRUE) 
FIT = glmQLFit(dge, design_matrix, robust=TRUE)
model = glmQLFTest(FIT, contrast = makeContrasts(contrasts = c("test-control"), levels = design_matrix)) 
output = data.frame(Gene = rownames(dge$counts), FoldChange = 2^-(model$table$logFC))

thanks in advance!

[ATpoint]

Unclear what the problem is. Do you run the analysis on only 200 genes? What do you mean by when there's a decrease in the test group?

[solarchan7]

Thanks for the reply. For example, the data is like: (https://ibb.co/8ckvMdx)

where the test group values are less than the control groups, but the fold change is 2! but the following line, the test group values are also less than the control groups but the fold change is ~0.5 which make way more sense.

[i.sudbery]

It looks like you are assuming the rows are in the same order in dge and in model$table. Are you sure this is the case?

[solarchan7]

ummmm..... not really. How can I check if they match up?

[Gordon Smyth]

They are in the same order.

Answer 1

[Gordon Smyth]

edgeR logFCs are always consistent with the counts and the data, so the title of your question can't be correct.

You can't judge expression changes just by looking at the matrix of counts, because the different library sizes for each sample need to be taken into account as well. To examine the data for each gene adjusting for library sizes, the simplest way is to compute the logCPM values:

logCPM <- cpm(dge, log=TRUE)

The logFC reported by edgeR are closely related to (but not the same as) the difference in mean logCPM between groups.

[solarchan7]

I'm a bit confuse :( . So if I want to get the fold change, I do:

logCPM = cpm(dge, log=TRUE) FC = 2^logCPM

? And this will fix the problem of different sign ?

[Gordon Smyth]

No, that's not what I mean. If you unlog logCPM you will just get CPM, not logFC.

I don't understand why you're not simply using topTags(model) which shows you the differentially expressed genes and their logFCs. I don't follow what else you are trying to do or what problem you see.

Your code FoldChange = 2^-(model$table$logFC)) is converting logFCs to inverse-fold-changes. I'm not sure why you would want to do that, but fold-changes and inverse-fold-changes are of course always positive.

[solarchan7]

I am trying to compare the test group (test1, test2, and test3) and the control group (control1, control2, control3, and control4). So yes, it is possible to get what I want with topTags(model); however, the problem is: (for example) using topTags(model), I got gene A and gene B having logFC -1.5 and 0.9 respectively, however, they SHOULD have the same sign according to their original data where both genes have a significant larger control groups than the test group.

[Gordon Smyth]

I'm not convinced that there is any problem. How have you judged that the genes SHOULD have the same logFC? None of the code you show in your question would allow you to make that judgement. What "original data" have you looked at and how?

[statmushroom]

• I DON'T know how to solve the problem, but I think I can clarify some points here. I think the OP assumes logFC = log(A/B) = log(test group/control group). So if the test group value is larger than the control group value, then it should result in log(some value > 1). Here, log(1) = 0. So if the test group value is larger than the control group value, then it SHOULD have POSITIVE value. However, OP seems to have both positive and negative values for genes that have test value > control value.

[Kevin Blighe]

solarchan7, if you believe there is something 'awry', then please paste a complete reproducible example, including code to generate sample input data and code to process such data such that the problem, as you interpret it, can be observed by all in this thread. Thanks!