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.
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):
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:
The cds data frame is:
The samples data frame:
which was used to create the deseq data frame:
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!
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.
Code:
> 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
Code:
> 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
The DESeq analysis goes smoothly (estimating size factors and dispersion) until I run nbinomTest. I am getting the following error:
Code:
> res = nbinomTest(cds,"CTL","KD") Error in nbinomTest(cds, "CTL", "KD") : For CountDataSets with multivariate conditions, only the GLM-based test can be used.
Code:
> 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
Code:
> 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
Code:
samplesDESeq = with(samples, data.frame(shortname = I(shortname), countf = I(countf), condition = condition, LibraryLayout = LibraryLayout)) cds = newCountDataSetFromHTSeqCount(samplesDESeq,directory=".")
Comment