Hi Elena, Thanks for posting the detail email. The p-values and p-adjusted values should be in your variable "test". You did: test <- testForDEU( ecs ) The function "testForDEU" will return an ExonCountSet object, your variable "ecs" with the columns of the fData for the pvalues and padjusted already filled. Then you will notice that "test" is an ExonCountSet object, so if you do: res1<- DEUresultTable(test) You will get your pvalues! Bext wishes, Alejandro Dear Simon, Thanks for your reply, even on a Saturday. =) You were right about the incorrect file path- that was problem. (To anybody else out there struggling with general issues getting their files into R, my code is below. Hope I save you some time). Another issue I was hoping you might help me with: I have no differentially-expressed exons (p-value = NA for all exons) at the end of the analysis!! I followed your other recommendations, and have spent some time trying to figure this out, but am still unsure. The fact that there are no integers in the p-value column has me worried. Maybe the problem is that the column names in counts(ecs) do not word-for-word match up with my rownames in the samples dataframe. I didn't know how to change the colnames in counts(ecs), so I left them as they were, noting that the sample identity seemed to have been preserved. Any advice on how to change the count(ecs) headers, or other advice about my non-existent p-values would be appreciated! =) Many thanks, Elena P.S. The graphics w/ gene models in DEXseq are very useful and convenient! ______________________________ library(DEXSeq) options(digits=3) setwd("C:\\Rdata\\DEXseq") library(DEXSeq) rm(list=ls()) # this GFF file is an output from Simon's script dexseq_prepare_annotations.py annotationfile = file.path("/Rdata/DEXseq/realFioles/DEXSeq_annotations.gff") annotationfile [1] "/Rdata/DEXseq/realFioles/DEXSeq_annotations.gff" > samples = data.frame( + condition = c("drug","drug","vehicle","vehicle"), + replicate = c(1:2,1:2), + row.names = c("drug1","drug2","veh1","veh2"), + stringsAsFactors = TRUE, + check.names = FALSE + ) > samples condition replicate drug1 drug 1 drug2 drug 2 veh1 vehicle 1 veh2 vehicle 2 > fullFilenames<- list.files("C:/Rdata/DEXseq/realFioles/",full.name s=TRUE,pattern="counts.txt") > fullFilenames [1] "C:/Rdata/DEXseq/realFioles/drug1_counts.txt" "C:/Rdata/DEXseq/realFioles/drug2_counts.txt" "C:/Rdata/DEXseq/realFioles/veh1_counts.txt" [4] "C:/Rdata/DEXseq/realFioles/veh2_counts.txt" > ecs<- read.HTSeqCounts(countfiles = fullFilenames,design = samples,flattenedfile = annotationfile) > head(counts(ecs)) C:/Rdata/DEXseq/realFioles/drug1_counts.txt C:/Rdata/DEXseq/realFioles/drug2_counts.txt C:/Rdata/DEXseq/realFioles/veh1_counts.txt 2L52.1:001 0 2 2 2L52.1:002 4 12 13 2L52.1:003 7 8 7 2L52.1:004 6 4 4 2L52.1:005 9 6 16 2L52.1:006 6 4 13 C:/Rdata/DEXseq/realFioles/veh2_counts.txt 2L52.1:001 4 2L52.1:002 20 2L52.1:003 8 2L52.1:004 7 2L52.1:005 14 2L52.1:006 12 > head(fData(ecs)) geneID exonID testable dispBeforeSharing dispFitted dispersion pvalue padjust chr start end strand transcripts 2L52.1:001 2L52.1 E001 FALSE NA 0.5145 0.5145 NA NA chrII 1867 1911 + 2L52.1 2L52.1:002 2L52.1 E002 TRUE 1.39e-01 0.0899 0.1391 NA NA chrII 2506 2694 + 2L52.1 2L52.1:003 2L52.1 E003 TRUE 9.56e-10 0.1436 0.1436 NA NA chrII 2738 2888 + 2L52.1 2L52.1:004 2L52.1 E004 TRUE 2.52e-10 0.2022 0.2022 NA NA chrII 2931 3036 + 2L52.1 2L52.1:005 2L52.1 E005 TRUE 2.16e-09 0.0979 0.0979 NA NA chrII 3406 3552 + 2L52.1 2L52.1:006 2L52.1 E006 TRUE 2.46e-09 0.1238 0.1238 NA NA chrII 3802 3984 + 2L52.1 # Size factors > ecs<- estimateSizeFactors(ecs) sizeFactors(ecs) # Dispersion > ecs<- estimateDispersions(ecs) # Fit Dispersion > ecs<- fitDispersionFunction(ecs) # Plot Individual exons via mean expression > meanvalues<- rowMeans(counts(ecs)) > plot(meanvalues, fData(ecs)$dispBeforeSharing, log="xy", main="mean vs CR dispersion") > x<- 0.01:max(meanvalues) > y<- ecs@dispFitCoefs[1] + ecs@dispFitCoefs[2] / x > lines(x, y, col="red") # Plot looks good # Test for Expression Difference > test<- testForDEU(ecs) # Seems to work > ecs<- estimatelog2FoldChanges(ecs) > res1<- DEUresultTable(ecs) head(res1) geneID exonID dispersion pvalue padjust meanBase log2fold(drug/vehicle) 2L52.1:001 2L52.1 E001 0.5145 NA NA 2.01 -0.844 2L52.1:002 2L52.1 E002 0.1391 NA NA 12.25 -0.304 2L52.1:003 2L52.1 E003 0.1436 NA NA 7.45 0.738 2L52.1:004 2L52.1 E004 0.2022 NA NA 5.22 0.595 2L52.1:005 2L52.1 E005 0.0979 NA NA 11.18 -0.264 2L52.1:006 2L52.1 E006 0.1238 NA NA 8.71 -0.586 > table(res1$padjust< 0.1) character(0) # Are there any p-values? > table(res1$pvalue< 1) character(0) [[alternative HTML version deleted]]