Question: How to specify 'contrast' in glmLRT command when treatments in design matrix are named with numbers, not letters
0
4.6 years ago by
MBWatson0
United States
MBWatson0 wrote:

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?

rnaseq edger glmlrt() • 1.4k views
modified 4.6 years ago by James W. MacDonald51k • written 4.6 years ago by MBWatson0
Answer: How to specify 'contrast' in glmLRT command when treatments in design matrix are
3
4.6 years ago by
United States
James W. MacDonald51k wrote:

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))
Answer: How to specify 'contrast' in glmLRT command when treatments in design matrix are
1
4.6 years ago by
Denali
Steve Lianoglou12k wrote:

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)
...

Not sure if it's a Safari thing, but I had a hell of a time with the simple formatting of the code in this comment -- hopefully we'll get a markdown editor back in here someday!

What was the issue? Was it that the code all looked like it was formatted as inline code? If so, are you using RStudio?

I was entering the answer all at once in the textbox itself (not copy/pasting from elsewhere). After that I'd highlight the relevant bits and change the formatting (ie. change "this" to inline code, change "that" to Code <blockstyle>, but it just wasn't behaving well. The insertion point was even jumping up to the top of the text box if I tried to delete whitespace in the middle of that endeavor.

Also: using RStudio? Puh-leeze ... emacs/ess or bust ... and I'd thank you if you refrain from repeating such an attempt to assassinate my character again in the future!

(j/k, for one reason or another I've been flirting with RStudio on occasion recently, and it's actually quite nice for some things)

1

Hmmm. A Mac guy who also uses emacs/ess? I didn't know you could do that - sort of like adding matter and antimatter together without the resulting annihilation. Do you also have a fusion reactor in your back yard?

Anyway, it seems to me that Safari is neck and neck with IE for 'most busted browser'. Does Chrome or Firefox work better?

The issue I see with RStudio (which I don't really use either - no need to change all the hard won keybinding muscle memory I have for using emacs/ess) is that a copy/paste results in everything being highlighted with what looks like inline code:

> really <- function(x) "why is this all greyed out?"
> really()
[1] "why is this all greyed out?"

Which annoys me to no end...

Answer: How to specify 'contrast' in glmLRT command when treatments in design matrix are
0
4.6 years ago by
MBWatson0
United States
MBWatson0 wrote:

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