RNA-seq Anova on three groups

Question

[BM]

Hi,

I would like some advice on RNA-seq analysis. I have three groups (control, mutant 1 and mutant 2). We would like perform a one-way anova/GLM and to analyse all three groups pair wise. So analyse all three groups to test for significance, then establish if any gene is changed bewteen any of the groups: control vs mutant 1; control vs mutant 2; mutant 1 vs mutant 2. However, we are running into problems.

Is it possible to do this in DESeq2 or EdgeR, or should iuse Genespring or Partek?

Thanks

[BM]

Thanks. But what is the formula for doing this in DESeq2. So will Deseq2 analyse all three groups at once or only two (which I thought was the case). And also if all three are analysed in ANOVA/ANODEV, are pairwise compaisions performed an p-values listed for all pairwise comparisons or an overall comparision?

[Devon Ryan]

Just read through the DESeq2 vignette. You would normally use a design like~group and just specify each pairwise comparison in a different contrast. You'll then get 3 separate sets of adjusted p-values, just as you would with an ANOVA. You could also do an overall comparison if you wanted to, for which you would use a likelihood ratio test with a reduced model of ~1 (again, just read through the vignette for the exact commands).

[BM]

The data has been aligned using TopHat and counts performed using Dexseq. I am trying to use Deseq2 for differential expression, but my experience of R is limited.

1. I am having difficulty constructing the metadata on each of the samples (columns of the count table). I can load the sample info, but I can't assign it to coldata. So I have used this instead:

> counts = read.delim("dexseqcount.txt", header=T, row.names=1)
> sample <- read.delim("~/sample.txt")
> column.data <- data.frame(genotype=as.factor(c("WT", "WT", "WT", "WT", "Mut1", "Mut1", "Mut1", "Mut2", "Mut2", "Mut2", "Mut2")))
> count.data.set <- DESeqDataSetFromMatrix(countData=counts, colData=column.data,design= ~ genotype)
> dds<-DESeq(count.data.set)
> res <- results(dds)

2. Also if I want to a pairwise comparison between all three groups (WT, Mut1, Mut2), as the above is a comparison between WT and Mut2 would this be ok:

> res2 <-results(dds, contrast=c("genotype","WT","Mut1"))
> res3 <-results(dds, contrast=c("genotype","Mut1","Mut2"))

3. I presume it is correct doing a pairwise comparison between all the three groups and looking to see if a gene is differentially expressed between WT v Mut2 and also in Mut1 vs Mut2 but not WT and Mut1. This will be done "by eye" by suing the three lists generated and seeing if there is any overlap

4. Final question, is it correct to do the above by doing the default test in Deseq2 i.e. Wald or should I use the alternative likelihood ratio test?

Thanks

[Michael Love]

please do not cross post

https://support.bioconductor.org/p/70285/

Answer 1

[Devon Ryan]

You can certainly do this in either DESeq2 or edgeR. You'll just specify the 3 contrasts, one for each comparison.

Answer 2

[Steven Lakin]

Post-hoc testing is required once you have constructed your GLM. The GLM will only estimate the dispersions amongst your variables; in order to get p-values, you have to perform testing based on the model. Probably the easiest thing to do is to take a look at the edgeR vignette under the section 2.9 More complex experiments (glm functionality); the package has the functions to do the pairwise post-hoc testing, and they have an example for you to follow at the end of the documentation.

http://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf