Build a data frame in order to implement a linear mixed model

Question

I have to compare two different methods of gene quantification. I have two matrices:

• the first matrix comes from the first method and has on the column the name of the samples and on the row the name of the gene;
• the second matrix come from the second method and has on the column the name of the samples and on the row the name of the gene;

Now from these two matrices I have to create one data frame for each gene where in the column I have: first column quantification, second column method, third column sample. The row correspond to the name of the gene. In total, I have to create a data frame for each gene according to these features.

Finally I have to use this data frame to run the linear mixed model to check the variability among samples and methods as here:

# Consider different source of variability, i.e., the samples and the methods
library("lme4")

# build the data.frame
df <- data.frame(express, method, sample)

# fit the lmm
res <- lmer(express~(1|method)+(1|sample), data=df)

# variance components
var_random_effect <- as.numeric(VarCorr(res))
var_residual <- attr(VarCorr(res),"sc")^2

Do you have any idea how to create these data frames and run the linear mixed model? I need to create a data frame for each gene.

Answer 1

First, I create two data frames that look like yours:

gene <- seq(1:10)
method <- rep(1,10)
sample1 <- rnorm(10, 500,150)
sample2 <- rnorm(10, 500,150)
sample3 <- rnorm(10, 500,150)
sample4 <- rnorm(10, 500,150)
method1 <- data.frame(gene, method, sample1, sample2, sample3, sample4)
method <- rep(2,10)
sample1 <- rnorm(10, 500,150)
sample2 <- rnorm(10, 500,150)
sample3 <- rnorm(10, 500,150)
sample4 <- rnorm(10, 500,150)
method2 <- data.frame(gene, method, sample1, sample2, sample3, sample4)

Note that I added a column to each, denoting the method used in each data frame. One of the data frames look like this:

method1
   gene method  sample1  sample2  sample3  sample4
1     1      1 476.7839 227.6009 592.9634 409.8754
2     2      1 414.5094 469.0531 373.1117 521.8902
3     3      1 376.0891 321.8406 448.1987 573.8746
4     4      1 665.7464 586.7600 383.0384 453.6320
5     5      1 449.8427 609.0248 440.6923 525.9350
6     6      1 439.3713 464.3511 446.8897 609.7765
7     7      1 601.1965 593.3730 465.6389 421.5950
8     8      1 744.8830 220.0788 303.6520 461.6706
9     9      1 423.5571 101.7189 169.9164 452.0107
10   10      1 605.3401 449.9921 611.2779 613.9308

Ok, so no we need to convert them to long form instead of wide. We want all of the samples to have a row on its own:

library(reshape2)
method1_long <- melt(method1, id.vars=c("gene", "method"))
method2_long <- melt(method2, id.vars=c("gene", "method"))
df<- rbind(method1_long, method2_long)
colnames(df)[c(3,4)] <- c("sample", "express")

Now the data looks like this (only rows 35-45 shown):

df[35:45,]
   gene method  sample  express
35    5      1 sample4 525.9350
36    6      1 sample4 609.7765
37    7      1 sample4 421.5950
38    8      1 sample4 461.6706
39    9      1 sample4 452.0107
40   10      1 sample4 613.9308
41    1      2 sample1 231.2023
42    2      2 sample1 377.3416
43    3      2 sample1 551.1926
44    4      2 sample1 737.7293
45    5      2 sample1 643.7356

As you can see, we now have a factor "sample" which denotes which sample each test came from, and a variable "express" which is your dependent variable. We also have "method" and "gene". In your formula, you wanted to use method as a random effect. Unless I misunderstood you completely, you only have two methods so it makes no sense to include it as a random, rather than fixed, effect. A general recommendation is at least 5-6 different levels of a factor to include it as a random effect. Besides, comparing the methods seems to be your main goal, so it makes even less sense to use it as a random effect.

Now you want to run the mixed model per gene. There is no need to create different data frames per gene to do that:

res_gene1 <- lmer(express ~ method + (1|sample), data=df[which(df$gene==1),])
res_gene2 <- lmer(express ~ method + (1|sample), data=df[which(df$gene==2),])
res_gene3 <- lmer(express ~ method + (1|sample), data=df[which(df$gene==3),])

And so on. I hope this helps!