Hi, I conducted a multi-level experiment with gene expression measures in multiple brain regions in 2 groups.

So for instance, my data for region 1 look like this (note that for this post, I renamed the rows: these are gene names in my data): (Note that I recoded missings so that I have at least 3 valid observations per group for every gene; otherwise, all observations for that gene are recoded to NA).

  Control.A Control.A.1 Control.A.2 Control.A.3    Case.A  Case.A.1  Case.A.2  Case.A.3
1        NA          NA          NA          NA       NA       NA       NA       NA
2        NA          NA          NA          NA       NA       NA       NA       NA
3        NA          NA          NA          NA       NA       NA       NA       NA
4        NA          NA          NA          NA       NA       NA       NA       NA
5  24.97941    25.50746    24.64801    25.49045 26.14554 25.77994 25.57589 25.84357
6  23.77846    23.86665    23.96921    23.60105 24.06159 23.44947 24.59412 24.40061

When I look at the topTable results for region 1, I see that there are sometimes values for logFC and AveExpr, but unaccompanied by test results.

For instance, my topTable results for the data shown above look like this:

      logFC  AveExpr        t    P.Value adj.P.Val         B
1        NA 25.06110       NA         NA        NA        NA
2        NA      NaN       NA         NA        NA        NA
3        NA 24.48188       NA         NA        NA        NA
4        NA 22.03522       NA         NA        NA        NA
5 0.6799025 25.47455 2.545508 0.01499063 0.1087153 -3.502281
6 0.3226035 24.04767       NA         NA        NA        NA

sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)

I do not understand this. First, I do not understand how I can sometimes have AVeExpr values but no logFC values. Second, if I have logFC and/or AveExpr values, why do I not obtain test results?

Any ideas on what went wrong?

kind regards
Sophie

Thank you for sending me your data offline. I am guessing that this is proteomics mass spectrometry data, which would explain the NAs.

The results you are seeing are correct, but can be improved. I see that your experiment has 6 brain regions from case and control individuals, so there are 12 distinct groups in your experiment (6 regions x 2 conditions).

1. limma will return an AveExpr value if any of the 12 groups contain any non-NA values. I think you are assuming that the non-NA values have be in the two groups that are being compared, but they do not. If you read the ?topTable help page you will see that AveExpr is the average log-expression for that gene across all 48 samples in your experiment. You will get the same AveExpr values no matter which contrast you test.

2. limma will return a logFC value if the two groups being compared both contain non-NA values.

3. If you use contrasts.fit, then limma needs all 12 groups to contain at least one non-NA value before it can compute t-statistics for any of the contrasts. You might think that limma should recognize that only some of the groups are involved in each contrast, but this is invisible to limma because of the nature of multiplication in R. Computing 0*NA in R gives NA, so limma needs all the coefficients to be non-NA in order to get an non-NA contrast value.

You can avoid the use of contrasts.fit like this. The results will be same as for your design matrix except that you will get a t-statistic and p-value whenever you have a logFC:

keep <- (rowSums(!is.na(exprs(ExprSet))) >= 4)
ExprSet <- ExprSet[keep,]
design <- model.matrix(~regionid+regionid:condition,data=pData(ExprSet))
corfit <- duplicateCorrelation(ExprSet, design, block = ExprSet$subject)
corfit$consensus.correlation
fit <- lmFit(ExprSet, design, block = ExprSet$subject, correlation = corfit$consensus)
fit <- eBayes(fit)
res.CBWM <- topTable(fit,coef="regionidA:conditionVWM",n=Inf)

It seems that your data has a huge number of NAs.

limma needs only one non-NA value to compute AveExpr but it needs at least two groups to have non-NA values before it can compute a logFC. Hence, if there are variable patterns of NA between genes, it is perfectly possible to get AveExpr without the other values.

On the other hand, I do not see how it is possible to get a logFC without a t-statistic. limma will usually give a t-statistic and p-value even for genes with no residual df (it does that by using the empirical Bayes prior). So I really would like to see the code that you used to generate the topTable.

Note that there is no need to recode any values to NA. limma will correctly use whatever non-NA values it finds and it is perfectly capable of using just 1 or 2 values in a group.

I don't know what technology you are using or what special issues there might be that lead to NAs but, in my experience, it is seldom necessary to have any NA values in a expression data set. If the data are from microarrays or RNA-seq, then NA values can almost always be avoided by appropriate processing of the data.

Correction

It is possible to get logFCs without t-statistics if you are using contrasts.fit, which I now see is happening in your code. I will elaborate on this in a separate answer.