Hello,

I am using edgeR to calculate differential gene expression among individuals from two populations (BR and SD) exposed to three salinity treatments (15, 35, 60 ppt), with three biological replicates each (a, b, c) for 18 samples total. When I create my design matrix, my six treatments are named with numbers, but in the edgeR documentation, the treatments are named with letters, for example see section 3.2.3 on page 24. Because of the way my treatments are named, I am unsure about how to use the 'lrt <- glmLRT(fit, contrast=c(-1,1,0))' command. If someone could suggest how to change my treatment names from numbers to letters (preferred since this will allow me to keep following the manual) OR give an example of how to use the contrast argument with numbers, I would appreciate it. Thank you! The design and my code are below.

DESIGN:

> design
          1 2 3 4 5 6
X1_BR15a  1 0 0 0 0 0
X2_BR15b  1 0 0 0 0 0
X3_BR15c  1 0 0 0 0 0
X4_BR35a  0 1 0 0 0 0
X5_BR35b  0 1 0 0 0 0
X6_BR35c  0 1 0 0 0 0
X7_BR60a  0 0 1 0 0 0
X8_BR60b  0 0 1 0 0 0
X9_BR60c  0 0 1 0 0 0
X10_SD15a 0 0 0 1 0 0
X11_SD15b 0 0 0 1 0 0
X12_SD15c 0 0 0 1 0 0
X13_SD35a 0 0 0 0 1 0
X14_SD35b 0 0 0 0 1 0
X15_SD35c 0 0 0 0 1 0
X16_SD60a 0 0 0 0 0 1
X17_SD60b 0 0 0 0 0 1
X18_SD60c 0 0 0 0 0 1
attr(,"assign")
[1] 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

MY CODE:

source("http://bioconductor.org/biocLite.R")

    biocLite("edgeR")
library("edgeR")

x <- read.delim("path_to_file", header=TRUE, row.names="Symbols") 
z <-x[ ,c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18)] 
keep <- rowSums(cpm(z)) >=18 
z<- z[keep,] 

group <- c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6) 

y <- DGEList(counts=z, group=group) 

design <- model.matrix(~0+group, data=y$samples)
colnames(design) <- levels(y$samples$group)

y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTrendedDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)

fit <-glmFit(y, design)
lrt <- glmLRT(fit, contrast=c(-1,1,0)) ## What to put for contrast?

You might consider using something different than numbers, just for sake of clarity. Something like BR_15, BR_35, etc is much easier to keep straight than trying to remember what sample/salinity combination is represented by each number.

You also don't need to estimate the dispersions separately any longer. Just use estimateDisp(), which handles all that for you. Unless of course you have some ancient R/Bioconductor install, in which case you should update and then use estimateDisp().

The contrast has to be a vector that is as long as your design matrix is wide. So if you have 6 groups, then you need a vector of six numbers. As an example,

lrt <- glmLRT(fit, contrast = c(-1,1,rep(0,4)))

will compare 2 - 1. In general you want the coefficients to add up to zero. You can also use fractions or decimals if you want to compare the mean expression of groups. In other words, if you want to test to see if the average of the BR samples are different from the SD samples, you could do

lrt <- glmLRT(fit, contrast = c(1/3,1/3,1/3,-1/3,-1/3,-1/3))

Instead of having a group vector like so:

group <- c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6)

Change the entries to meaningful names, eg. from the description of your experiment, the group entries might be something like:

group <- c(
  "BR15", "BR15", "BR15",
  "BR35", "BR35", "BR35", 
  "BR60", "BR60", "BR60", 
  "SD15", "SD15", "SD15", 
  "SD35", "SD35", "SD35",​
  "SD60", "SD60", "SD60")

Then have at it with the rest of your code:

y <- DGEList(counts=z, group=group)
design <- model.matrix(~0+group, data=y$samples)
y <- estimateGLMCommonDisp(y, design)
...

Thank you! The answer seems very obvious now so thank you for your patience with my general R inexperience.