Hi everybody,

I am trying to compare gene expression using RNA Seq data from 50 and 50 human tissue samples with two phenotypically similar diseases that may exhibit distinct pathophysiological differences on the transcriptomic level (= my hypothesis).

Several, but not all subjects in the study provided more than one sample, in one case even one from each disease type. Thus, I would want to correct for the subject as a random effect.

Additional (fixed) covariates include gender (disease B has more men that A) and a measure of disease severity (continuous integer values).

In order to be able to incorporate both my fixed and random factors (as a blocking factor), I used limma/voom

design <- model.matrix(~1+ targets$disease + targets$gender + targets$severity)
y <- voom(x,design,plot=F, normalize="quantile")
corfit <- duplicateCorrelation(y, design, block=targets$Patient)
y <- voom(x,design,plot=FALSE, block=targets$Patient, correlation=corfit$consensus, normalize="quantile")
fit=lmFit(y, design, block=targets$Patient, correlation=corfit$consensus)
fit <- eBayes(fit)

I have also tried TMM normalization, or voomWithQualityWeights.

My problem is that I get a number of DE genes for the disease B vs. A with large FCs (up to 10fold) that are mostly Y-linked, and/or ribosomal proteins. These results do not make sense biologically.

As a comparison, I have used DESeq2 to analyze the same data. Correcting for severity and patient as fixed factors (yes, I know, patient should be random...), I find very different genes to be differentially expressed, with lower FCs. Those make more sense biologically.

I know that DESeq2 is probably not appropriate for my data, and I would like to use limma here. But what am I doing wrong here? Is it a normalization problem?

Many thanks in advance for your thoughts and comments!

Have you tried examining the expression profile of some of those genes (e.g., PRKY) across all disease/sex combinations? This can provide a useful sanity check in cases like this.

1. Apply some simple normalisation, e.g., edgeR's cpm normalisation with a prior of 3.
2. For a particular gene, make a box/violin/swarm plot (depending on what visualisation method you prefer) of expression values for each disease/sex(/severity) combination.
3. Compare across combinations to see if the change in expression is genuinely present.

The reported log-fold changes from topTable are quite large, so I would expect that any changes in expression should be quite apparent in the raw data. If that's the case, then it's understandable that limma should report these genes as being DE. That fact that they don't make biological sense is hardly limma's fault, it just works with the data that's been provided. The only potential issues I can forsee are:

• You've mislabelled the male/female samples. This is easy enough to check based on XIST expression.
• The disease/sex effects are not additive. Switch to a disease/sex interaction model to deal with this.

P.S. Putting patient as a fixed effect is not wrong per se, but in an experimental set-up like this, it means that your disease comparison will be performed using only those patients who contributed samples to both disease types.

You're quite right that you can't analyse the data with patient as a fixed effect. If you do that, then the disease effect will be estimated entirely from the one patient who provided a sample for both diseases, and all other patients will be ignored.

You also can't correct for the gender effect simply by including gender as an additive effect because there is apparently an interaction between gender and disease (males having more severe disease). So I'm not surprised that the Y chromosome genes are near the top of your toptable.

There are two possible approaches. Either you can estimate the disease effect for males and females separately, for example by using the model:

gender + gender:disease

Or alternatively, you might filter the sex-linked genes (Y chromosome and Xist) from your analysis and ignore gender in your analysis. We do that frequently in our own in-house analyses when samples are not balanced by gender.

The first approach may be better if the disease very different in males and females. The second approach may be better if the transcriptional profile of the disease at a given severity is much the same in males and females even though it may have arisen somewhat differently.

There also may be complicated interactions between disease and severity, but it's hard to comment on that.

I would also include the options robust=TRUE when you run eBayes().