EdgeR: Contrasts for both within and between subject comparison

[fh_fingolfin]

I have two sets of samples, one set being sensitive to a compound treatment, the other resistant. For each group I have patient samples that were left untreated, treated 2h with a compound and 24h with a compound.

I am trying to do resistant vs sensitive comparisons on the data using edgeR. The manual has an example of a comparison both between and within subjects in chapter 3.5 which to me applies 1:1 on my design.

Following the manual, I have set up the following design

design = model.matrix(~resist+resist:sub+resist:treat)

(Intercept) resistY resistN:sub2 resistY:sub2 resistN:sub3 resistY:sub3
1 1 0 0 0 0 0
2 1 0 0 0 0 0
3 1 0 0 0 0 0
4 1 0 1 0 0 0
5 1 0 1 0 0 0
6 1 0 1 0 0 0
7 1 0 0 0 1 0
8 1 0 0 0 1 0
9 1 0 0 0 1 0
10 1 1 0 0 0 0
11 1 1 0 0 0 0
12 1 1 0 0 0 0
13 1 1 0 1 0 0
14 1 1 0 1 0 0
15 1 1 0 1 0 0
16 1 1 0 0 0 1
17 1 1 0 0 0 1
18 1 1 0 0 0 1
resistN:treat24h resistY:treat24h resistN:treat2h resistY:treat2h
1 0 0 0 0
2 0 0 1 0
3 1 0 0 0
4 0 0 0 0
5 0 0 1 0
6 1 0 0 0
7 0 0 0 0
8 0 0 1 0
9 1 0 0 0
10 0 0 0 0
11 0 0 0 1
12 0 1 0 0
13 0 0 0 0
14 0 0 0 1
15 0 1 0 0
16 0 0 0 0
17 0 0 0 1
18 0 1 0 0

Now I understand that with the contrasts c(0,0,0,0,0,0,0,0,-1,1) and c(0,0,0,0,0,0,-1,1,0,0) you can test for differentially expressed genes between resistant and sensitive cells after 2h and 24h treatment, respectively.

What is unclear to me is for example how to test for differentially expressed genes between resistant and sensitive cells in the control?

It cannot be c(0,1,0,0,0,0,0,0,0,0), since that would give me differentially expressed genes between resistant and sensitive cells in any state right?

The only workaround I have come up so far is to relevel the treatment factors and take treat2h as base level in order to be able to have now the control in the contrast in its position. But I tried it and this changes the outcome of the treat24h contrast, so apparently there's something I'm missing.

Does anyone have any experience on that? Any help would be greatly appreciated...

[Gordon Smyth]

This question has since been answered on the Bioconductor mailing list:

https://stat.ethz.ch/pipermail/bioconductor/2014-July/060577.html

[vkartha]

On a similar note - I am trying to perform differential expression contrasts for 2 different populations - one resistant to a given drug and one wild type (i.e. sensitive). Within each population, I have samples that were treated at just one given concentration of the drug, and samples that were treated with a drug vehicle (control).

So the basic study design for this experiment is :

Population A (resistant):Treatment (Y) vs Control (X)
Population B (wild type): Treatment (Y) vs Control (X)

In addition to population and treatment variables, I also have to adjust for a very strong batch effect (3 batches).

So the formula for my design matrix at the moment is:

~ Batch + Population + Treatment + Population:Treatment

and the resulting design matrix looks like so:

(Intercept) Batch2 Batch3 PopulationB TreatmentY PopulationB:TreatmentY
Sample_1 1 0 0 0 0 0
Sample_2 1 0 0 0 0 0
Sample_3 1 0 0 0 0 0
Sample_4 1 1 0 0 0 0
Sample_5 1 1 0 0 0 0
Sample_6 1 0 0 0 1 0
Sample_7 1 0 0 0 1 0
Sample_8 1 0 0 0 1 0
Sample_9 1 1 0 0 1 0
Sample_10 1 1 0 0 1 0
Sample_11 1 1 0 1 0 0
Sample_12 1 1 0 1 0 0
Sample_13 1 0 1 1 0 0
Sample_14 1 0 1 1 0 0
Sample_15 1 0 1 1 0 0
Sample_16 1 1 0 1 1 1
Sample_17 1 1 0 1 1 1
Sample_18 1 0 1 1 1 1
Sample_19 1 0 1 1 1 1
Sample_20 1 0 1 1 1 1

What I don't understand is given this complete linear model formula, how do I specify contrasts (both edgeR and DESEq) to compare the two treatment methods WITHIN each population?

i.e. I want to test the effect of treatment separately for population A and for population B, adjusting for Batch effect, using this complete model.

Specifying coef=4 would give me the overall difference in expression between the two populations and specifying coef=5 will give me the overall difference in expression between the two treatment methods, but that is less interesting a question given my study design

Any help would be greatly appreciated


Regards

[TRA]

Hi all,
I mapped my RNA-seq reads using edgeR and got around 90% of mapping, but when I entered the feature count commands I got only 70% successfully assigned reads. Does any one have an idea?
Note: I have 18 samples, 3 per group. I did single end read sequencing.

Any help would be appreciated

[dpryan]

Just because a read aligns to the genome doesn't mean that it aligned to a gene. Perhaps it came from an intron or an intergenic region (e.g., due to transcriptional noise).

[TRA]

Thank you appreciate it

[TRA]

Hi all,
I am using edgeR flow to get gene lists of multifactorial experiment. I have 6 conditions with 3 samples per each. The problem is some of the gene lists I got end at FDR much less than what I entered and so I did not get the whole gene list. For example, the upregulated genes were completely missing in some lists. I tried increasing the number of genes the program gives using the
< gene_table<-topTags(lrt, n=40000)
but did not help. Then I tried very restreignent FDR as follow
< DE<-tp$table[tp$table$FDR<0.000001,]$GeneID
But again did not help. The last gene in one of my lists was with
Symbol Chr logFC logCPM LR PValue FDR
Hpx 7 -7.186931303 13.83923493 57.6010426 3.21E-014 3.23E-014

Any help would be appreciated!