user2294316 user2294316 - 1 day ago 5
R Question

Limma package to find differentially expressed genes

I have my matrix designed in the following way which I name as mat1:

Probes sample1 sample1 sample2 sample2 sample3 sample3 sample4 sample4
rep1 rep2 rep1 rep2 rep1 rep2 rep1 rep2
------------------------------------------------------------------------
gene1 5.098 5.076 5.072 4.677 7.450 7.456 8.564 8.555
gene2 8.906 8.903 6.700 6.653 6.749 6.754 7.546 7.540
gene3 7.409 7.398 5.392 5.432 6.715 6.724 5.345 5.330
gene4 4.876 4.869 5.864 5.981 4.280 4.290 4.267 4.255
gene4 3.567 3.560 3.554 3.425 8.500 8.564 6.345 6.330
gene5 2.569 2.560 8.600 8.645 5.225 5.234 7.345 7.333


I use the limma package to find the DEG's

Group <- factor(c("p1", "p1", "p2", "p2","p3", "p4","p4")
design <- model.matrix(~0 + Group)
colnames(design) <- gsub("Group","", colnames(design))
fit <- lmFit(mat1[,1:4],design)
contrast.matrix<-makeContrasts(p1-p2,levels=design)
fit2<-contrasts.fit(fit,contrast.matrix)
fit2<-eBayes(fit2)
sel.diif<-p.adjust(fit2$F.p.value,method="fdr")<0.05
deg<-mat1[,1:4][sel.diif,]


So will "deg" just give me those genes which are significant in sample one versus two. I am interested in those genes which are differentially expressed only in first sample but not in the second sample and am not sure if this is the right approach.

Or should I try something like this:

contrast.matrix<-makeContrasts(contrasts="p1"-("p2"+"p3"+"p4")/3,levels=design)


I am not sure how I should set the contrasts matrix to obtain the DEG's from sample 1 only but not in the other three.

Answer

Your example isn't reproducible, i.e. I can't reproduce the results. However, here are a few comments:

  1. You are correct regarding deg. It will look for genes that are different between the two samples.
  2. A contrast matrix of:

    makeContrasts(contrasts="p1-(p2+p3+p4)/3", levels=design)
    

    is how I would (probably) tackle this problem. However, this may cancel out effects. For example if p2 was high and p3 was low.

  3. Alternatively, you could have something like:

    makeContrasts(contrasts=c("p1-p2", "p1-p3", "p1-p4"), levels=design)
    

    and look at overlapping genes.

Comments