Question: RNAseq: how to make design matrix for limma
0
gravatar for vw
7 months ago by
vw20
vw20 wrote:

Hi all,

I want to compare the result between limma + voom with DESeq2. But I don't know how to make the design matrix for my date when I am using Limma.

I have two groups: WT, ABX. Each group has 6 samples: X1~X6(WT), X7~X12(ABX). It looks like:

head(exprSet, n=1)
                        X.1     X.2       X.3     X.4       X.5     X.6       X.7     X.8
ENSMUSG00000051951      88      113       83      101       57       45       45       50
                         X.9       X.10      X.11      X.12
ENSMUSG00000051951       31        51        59        43

The row name is the geneID, the column is the count for each sample.

My code is below:

## read count matrix from file
exprSet=read.table("ex_matrix_g1g2_h.txt", sep="\t", header=T, row.names = 1, stringsAsFactors = F)
group_list <- factor(c(rep("WT",6), rep("ABX",6)))

design <- model.matrix(~0+group_list) ## Is this correct?
colnames(design) <- levels(group_list)
rownames(design) <- colnames(exprSet)

v <- voom(exprSet,design,normalize="quantile", plot = T)
fit <- lmFit(v,design)
fit2 <- eBayes(fit)
tempOutput = topTable(fit2, n=Inf, adjust.method = 'BH', coef=2)  
## if I want to compare ABX - WT, how to set coef?

Unfortunately, the above code generated very different results compared with DESeq2. In limma's result, almost all genes are significantly different.

I think that maybe I used the wrong design matrix. But I am not familiar with limma. Could some help me to figure out what the matter of my code?

Thank you so much!

rna-seq limma deseq2 • 353 views
ADD COMMENTlink modified 7 months ago by kristoffer.vittingseerup3.0k • written 7 months ago by vw20
2
gravatar for kristoffer.vittingseerup
7 months ago by
European Union
kristoffer.vittingseerup3.0k wrote:

That is because you do not compare the two conditions but test whether ABX is expressed in the first place (whether ABX gene expression different from zero) . It's quite easy to fix - try this instead:

exprSet=read.table("ex_matrix_g1g2_h.txt", sep="\t", header=T, row.names = 1, stringsAsFactors = F)

group_list <- factor(
    x = c(rep("WT",6), rep("ABX",6)),
    levels=c("WT","ABX")                     # Set WT to be the first level
)

design <- model.matrix(~group_list)          # Remove the zero

v <- voom(exprSet,design,normalize="quantile", plot = T)
fit <- lmFit(v,design)
fit2 <- eBayes(fit)

tempOutput = topTable(
    fit2, n=Inf,
    adjust.method = 'BH',
    coef='group_listABX'                     # Easier to read in 6 months than coef=2.
)

By setting WT to be the first level and removing the zero from the design WT will be set to be the intersect - and then when you test the "group_listABX" you are comparing ABX to the intersect (aka you compare WT -> ABX). This also means that positive log2FC values are higher expressed in ABX.

I assumed that you want to do the WT -> ABX comparison. If you did not set WT to be the first levels ABX would be that instead (A is before W in the alphabet) meaning you would compare ABX -> WT, which would mean positive log2FC values would be higher expressed in WT.

Hope this helps.

Cheers Kristoffer

ADD COMMENTlink written 7 months ago by kristoffer.vittingseerup3.0k

Thank you so much! You solved my question perfectly. Now I know the meaning of the limma's design matrix.

ADD REPLYlink written 7 months ago by vw20
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1718 users visited in the last hour