View Single Post
Old 04-29-2012, 10:06 AM   #1
r_sitaram
Member
 
Location: Helsinki

Join Date: Apr 2012
Posts: 10
Default DEXSeq output not convincing

Hello,
I'm running DEXSeq for finding out the differential exon usage between 3 Wild type and 3 transgenic rat data. I obtained the exon count using the scripts dexseq_count.py and the gff file using dexseq_prepare_annotation.py. After following all the steps given in your vignette,When I run the DEUresultTable command, I get a lot of NA's in the pvalue column and just 1's in the padjust column. Also after executing the estimateDispersions function, the column dispersion_CR_est is showing NULL. I have a feeling that I'm making a mistake in the exoncount data object but I'm not sure. Could someone throw some light on this?

I could probably highlight the steps(in bold) I have done.

annotationfile = file.path("/triton/ics/scratch/rnaseq/exoncount/Rattus_norvegicus.RGSC3.4.66_DEXSEQ.gff")

annotationfile
[1] "/triton/ics/scratch/rnaseq/exoncount/Rattus_norvegicus.RGSC3.4.66_DEXSEQ.gff"

samples=data.frame(condition=c(rep("wild",3),rep("trans",3)),replicate=c(1,1,1,1,1,1),stringsAsFactors = TRUE,check.names = FALSE)
row.names(samples) <- c("01_5b_WT.txt","02_6b_WT.txt","03_15b_WT.txt","04_4b_TG.txt","05_7b_TG.txt","06_16b_TG.txt")


samples
condition replicate
01_5b_WT.txt wild 1
02_6b_WT.txt wild 1
03_15b_WT.txt wild 1
04_4b_TG.txt trans 1
05_7b_TG.txt trans 1
06_16b_TG.txt trans 1

countfiles = file.path("/triton/ics/scratch/rnaseq/exoncount",paste(row.names(samples))

countfiles
[1] "/triton/ics/scratch/rnaseq/exoncount/01_5b_WT.txt"
[2] "/triton/ics/scratch/rnaseq/exoncount/02_6b_WT.txt"
[3] "/triton/ics/scratch/rnaseq/exoncount/03_15b_WT.txt"
[4] "/triton/ics/scratch/rnaseq/exoncount/04_4b_TG.txt"
[5] "/triton/ics/scratch/rnaseq/exoncount/05_7b_TG.txt"
[6] "/triton/ics/scratch/rnaseq/exoncount/06_16b_TG.txt"

ecs = read.HTSeqCounts(countfiles,design=samples,flattenedfile=annotationfile)

ExonCountSet (storageMode: environment)
assayData: 236327 features, 6 samples
element names: counts
protocolData: none
phenoData
sampleNames: 01_5b_WT.txt 02_6b_WT.txt ... 06_16b_TG.txt (6 total)
varLabels: sizeFactor condition replicate countfiles
varMetadata: labelDescription
featureData
featureNames: ENSRNOG00000000001:001 ENSRNOG00000000001:002 ...
ENSRNOG00000045521:001 (236327 total)
fvarLabels: geneID exonID ... transcripts (13 total)
fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:

sampleNames(ecs) = row.names(samples)

ecs <- estimateSizeFactors(ecs)


sizeFactors(ecs)
01_5b_WT.txt 02_6b_WT.txt 03_15b_WT.txt 04_4b_TG.txt 05_7b_TG.txt
1.6794610 0.5407073 0.8215439 1.7972423 0.8516481
06_16b_TG.txt
0.9878676

ecs <- estimateDispersions(ecs)


featureData(ecs)$dispersion_CR_est
NULL

ecs <- fitDispersionFunction(ecs)

ecs <- testForDEU(ecs)

res <- DEUresultTable(ecs)
r_sitaram is offline   Reply With Quote