DESeq Error -- multivariate conditions with non-multivariate data

[muppetjones]

I just started doing RNA-seq analysis, and I'm trying to replicate some DESeq results. Long story short, my results were awful, and I can't install the 'pasilla' package (don't ask...) to test DESeq using the vignette.

So...I am now trying to validate my install using the method described in S. Anders et al. Nature Protocols. 2013. In an effort to be brief, I'll only post the results that differed from the paper. If more information is needed, please ask.

I have downloaded the data and compiled everything with 2 notable exceptions:
1) I am using a new build of the D. melanogaster genome (BDGP5.74 instead of BDGP5.70). I don't expect any major issues here.
2) I am using an older version of R, which could be causing problems. Unfortunately, I have not been able to successfully install a newer version due to serious proxy issues.

> sessionInfo()
R version 2.15.0 (2012-03-30)
Platform: x86_64-redhat-linux-gnu (64-bit)
...
other attached packages:
[1] edgeR_3.0.8         limma_3.14.4        DESeq_1.10.1       
[4] lattice_0.20-6      locfit_1.5-9.1      Biobase_2.18.0     
[7] BiocGenerics_0.4.0  BiocInstaller_1.8.3

I have managed to successfully follow the edgeR analysis without any issues. All of the data matched exactly, with the exception of the following (step 14-B-vi, Use the topTags function to present a tabular summary of the differential expression statistics):

> tt = topTags(de, n = nrow(d)) 
> head(tt$table)
                logFC   logCPM       PValue          FDR
FBgn0039155 -4.646404 5.868219 1.333990e-87 9.599391e-84
FBgn0025111  2.929648 6.855564 1.156645e-54 4.161607e-51
FBgn0003360 -3.227304 8.420833 1.052926e-52 2.525618e-49
FBgn0039827 -3.988696 4.386542 6.595629e-52 1.186554e-48
FBgn0000071  2.718231 4.724487 9.833718e-39 1.415269e-35
FBgn0034736 -3.515391 4.117020 2.467228e-38 2.959028e-35

which had some additional rows than the data given in the paper and different values for the rows that were the same. The next step (vii, Inspect the depth-adjusted reads per million for some of the top differentially expressed genes) had the same extra rows, but the values matched the paper. I imagine these changes are likely due to the newer version of the genome & feature file.

The DESeq analysis goes smoothly (estimating size factors and dispersion) until I run nbinomTest. I am getting the following error:

> res = nbinomTest(cds,"CTL","KD")
Error in nbinomTest(cds, "CTL", "KD") : 
  For CountDataSets with multivariate conditions, only the GLM-based test can be used.

The cds data frame is:

> cds
CountDataSet (storageMode: environment)
assayData: 15682 features, 7 samples 
  element names: counts 
protocolData: none
phenoData
  sampleNames: CT.SI.1 CT.PA.2 ... CT.SI.7 (7 total)
  varLabels: sizeFactor condition LibraryLayout
  varMetadata: labelDescription
featureData
  featureNames: FBgn0000003 FBgn0000008 ... FBgn0264727 (15682 total)
  fvarLabels: disp_pooled
  fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:  
> cds$condition
[1] CTL CTL CTL KD  KD  KD  CTL
Levels: CTL KD

The samples data frame:

> samples
    LibraryName LibraryLayout
1   Untreated-3        PAIRED
2   Untreated-4        PAIRED
3 CG8144_RNAi-3        PAIRED
4 CG8144_RNAi-4        PAIRED
5   Untreated-1        SINGLE
6 CG8144_RNAi-1        SINGLE
7   Untreated-6        SINGLE
                                                                                           fastq1
1                                                             SRR031715_1.fastq,SRR031714_1.fastq
2                                                             SRR031717_1.fastq,SRR031716_1.fastq
3                                                             SRR031725_1.fastq,SRR031724_1.fastq
4                                                             SRR031727_1.fastq,SRR031726_1.fastq
5 SRR031713.fastq,SRR031712.fastq,SRR031711.fastq,SRR031710.fastq,SRR031709.fastq,SRR031708.fastq
6 SRR031723.fastq,SRR031722.fastq,SRR031721.fastq,SRR031720.fastq,SRR031719.fastq,SRR031718.fastq
7                                                                 SRR031729.fastq,SRR031728.fastq
                               fastq2 condition shortname              countf
1 SRR031715_2.fastq,SRR031714_2.fastq       CTL   CT.PA.1   Untreated-3.count
2 SRR031717_2.fastq,SRR031716_2.fastq       CTL   CT.PA.2   Untreated-4.count
3 SRR031725_2.fastq,SRR031724_2.fastq        KD   KD.PA.3 CG8144_RNAi-3.count
4 SRR031727_2.fastq,SRR031726_2.fastq        KD   KD.PA.4 CG8144_RNAi-4.count
5                                           CTL   CT.SI.5   Untreated-1.count
6                                            KD   KD.SI.6 CG8144_RNAi-1.count
7                                           CTL   CT.SI.7   Untreated-6.count

which was used to create the deseq data frame:

samplesDESeq = with(samples, 
	data.frame(shortname = I(shortname), 
	countf = I(countf), 
	condition = condition, 
	LibraryLayout = LibraryLayout))
cds = newCountDataSetFromHTSeqCount(samplesDESeq,directory=".")

I have followed the code as provided by the paper, and while I know my way around code, I am new to this pipeline. Am I missing something? Any help would be greatly appreciated! Thanks!

[dpryan]

You have two types of conditions ("condition" and "LibraryLayout") in samplesDESeq, making it multivariate. Use the GLM methods as the error message suggests.