We have an RNA-Seq data set from mouse with three conditions in triplicates. For better understanding of the reaction, each of the animal was weighted and its urea levels were measured beforehand.

We would like to apply a linear mixed-model to the data-set in order to analyze the influence of weight and urea conc. on the expression intensities of genes in my RNA-Seq experiment.

The weight.score and urea.score were calculated based on multiple measurements on different days

We're trying to identify genes, for which a correlation between the urea conc., body weight and their expression behavior exits.

Our first try was to apply a Pearson correlation to the total.score, we calculated for each of the genes (see script below). But in the review we were asked to do to apply a mixed linear-model to test the influence of the body weight and urea concentration on those expressions values and to examine the degree of contribution of each gene to the outcome score. and here I would like to ask for your help.

We would like to include both the urea.score and the weight.score in the mixed-model (Or do I need to include the actual values in the model?). I'm not sure how to do so.

Our question is - Is there an association between the body weight and urea concentration and the intensities of the gene expression?

In the model I regards the gene expression as my dependent variable, while urea.score and weight.score are the random variables. Then fitting the mixed model with these effects would be something like that lmer(Gene.X ~ (1|urea.score) + (1|weight.score), data = mergedCorr.norm).

But when running the linear model for urea and body-weight separately I can see a significance, while in the mixed-model this is no longer apparent. Am I setting the parameters correctly? Can I sue these two parameters as random variables, even though I have only three different levels?

I would appreciate any help or suggestions as to how I should apply these here.

Script used:

# prepare data

Scores.RiboM <- structurelistAnimal.ID = structure(c(5L, 6L, 7L, 4L, 9L, 10L, 
2L, 1L, 8L, 3L), .Label = c("C20", "C22", "C24", "CR4.1", "CR4.2", 
"CR4.3", "CR4.4", "HP11", "HP12", "HP14"), class = "factor"), 
group = c(3L, 3L, 3L, 3L, 2L, 2L, 1L, 1L, 2L, 1L), urea = structure(c(8L, 
10L, 3L, 2L, 1L, 4L, 6L, 5L, 9L, 7L), .Label = c("105.57", 
 "110.71", "125.89", "195.87", "321.43", "329.92", "340.81", 
 "35.59", "363.37", "94.76"), class = "factor"), weight = structure(c(9L, 
 6L, 4L, 2L, 8L, 7L, 3L, 5L, 10L, 1L), .Label = c("20", "20.1", 
 "20.7", "21.9", "22.1", "22.6", "23", "23.3", "23.5", "24"
 ), class = "factor"), urea.score = c(1L, 1L, 1L, 1L, 1L, 
2L, 3L, 4L, 4L, 4L), weight.score = c(1L, 1L, 1L, 2L, 3L, 
3L, 3L, 3L, 3L, 4L), total.score = c(2L, 2L, 2L, 3L, 4L, 
5L, 6L, 7L, 7L, 8L)), row.names = c(8L, 9L, 10L, 7L, 5L, 
6L, 2L, 1L, 4L, 3L), class = "data.frame")

# scale() centers the data (the column mean is subtracted from the values in the column) and then scales it (the centered column values are divided by the column’s standard deviation).
gene.scaled <- structure(c(1.14326057352769, -0.331767048798829, 2.0634403441122, 
0.0974979558410874, -0.477899816335822, 0.634079211640983, -0.487033114306884, 
-1.06014756199103, -0.982514529237001, -0.598916014452395), .Dim = c(10L, 
1L), .Dimnames = list(c("C20", "C22", "C24", "HP11", "HP12", 
"HP14", "CR4.1", "CR4.2", "CR4.3", "CR4.4"), "Gene.X"))

# Correlation
mergedCorr.norm <- merge(Scores.RiboM, gene.scaled, by.x=1, by.y=0)
mergedCorr.norm <- mergedCorr.norm[order(mergedCorr.norm$total.score, decreasing = FALSE),]
CorValue.pear <- cor.test(x=mergedCorr.norm$total.score, y=mergedCorr.norm[,8], method = "pearson")
CorValue.spear <- cor.test(x=mergedCorr.norm$total.score, y=mergedCorr.norm[,8], method = "spearman")

# Plotting
ggplot(mergedCorr.norm, aes(x=total.score, y=mergedCorr.norm[,8], label=mergedCorr.norm$Animal.ID)) +
        geom_point(shape=1) +
        geom_smooth(method=lm) +
        scale_y_continuous(name="norm. read counts")+
        geom_text(size=3, hjust=-0.1, vjust=-0.4,colour="darkgreen")+ 
        annotate("text", x = 8, y = max(mergedCorr.norm[,8]),label=paste("Pear. r= ", format(CorValue.pear$estimate, digits = 4), "\n Spear. r = ", format(CorValue.spear$estimate, digits = 4), sep=""), colour=ifelse(CorValue.pear$estimate>=0.9, "blue","red")) +
        theme_bw(base_size = 12, base_family = "")

# linear model 
lm.weight <- lm(Gene.X ~ weight.score, data = mergedCorr.norm)
lm.urea <- lm(Gene.X ~ urea.score, data = mergedCorr.norm)

# mixed-model
mixed.Exp <- lmer(Gene.X ~ (1|urea.score) + (1|weight.score), data = mergedCorr.norm)
stargazer(mixed.Exp, type = "text",
    digits = 3,
    star.cutoffs = c(0.05, 0.01, 0.001),
    digit.separator = "")
  • 1
    $\begingroup$ If you can display the (faked) data from one mouse, it will be helpful. $\endgroup$
    – user158565
    Aug 8, 2019 at 13:51

2 Answers 2


RNA-seq data are (normalized) count data exhibiting a mean-variance relationship. The linear mixed model works for normal data/error terms that don’t have this relationship. You could work with a Poisson or Negative Binomial mixed model instead.


Exactly, as Dimitris mentioned, RNA-Seq tends to follows a Poisson distribution or a reparameterized Negative Binomial, where the NB now models "over-dispersed Poisson" (variance > mean) data rather than the traditional "trials until failure" model.

There are a few different R packages in Bioconductor for RNA-Seq analyses that might be worth checking out. These techniques focus on finding differentially expressed genes, and some provide the flexibility to model interactions (i.e., Weight X Reaction), or paired data (i.e., before/after response for a specific mouse), or more complex relationships. They might be able to offer some functionality to address your needs.

Below are a few methods with their citation and tutorials.



Robinson, Mark D., Davis J. McCarthy, and Gordon K. Smyth. "edgeR: a Bioconductor package for differential expression analysis of digital gene expression data." Bioinformatics 26.1 (2010): 139-140.



Tarazona, S., García, F., Ferrer, A., Dopazo, J., & Conesa, A. (2011). NOIseq: a RNA-seq differential expression method robust for sequencing depth biases. EMBnet. journal, 17(B), 18-19.



Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome biology, 15(12), 550.

Hope this helps!


Not the answer you're looking for? Browse other questions tagged or ask your own question.