Cannot reproduce the DEXSeq example

[Jerry_Zhao]

Hi, All,

I recently try to perform alternative splicing analysis on my RNA-seq data sets.
I tried DEXSeq package. However, I got almost all NAs in the padj column for my data sets.

Therefore, I planed to reproduce the example in DEXSeq library, but I have problems to reproduce the example (pasilla).
I used the exact code in the DEXSeq manual as follows, but I got an error message at the near last step.
When I tried to get the number of significant genes using "table ( dxr1$padj 0.1 )", I got an error message "Error: unexpected numeric constant in "table ( dxr1$padj 0.1" ".

I think this may due to the NAs in the padj column. Is there any way to deal with this?
When comparing my results with the results in the DEXSeq manual, the only difference I can tell is that: the numbers in the columns of exonBaseMean and exonBaseVar in the manual is integer ( such as 64 and 1251)
whereas in my analysis has four significant numbers (64.0000 and 1250.667).

Codes:
inDir = system.file("extdata",package="pasilla")
countFiles= list.files(inDir, pattern="fb.txt$",full.names=TRUE)
countFiles
flattenedFile = list.files (inDir, pattern="gff$",full.names=TRUE)
flattenedFile

sampleTable = data.frame(
row.names = c("treated1","treated2","treated3","untreated1","untreated2","untreated3","untreated4"),
condition = c("knockdown","knockdown","knockdown","control","control","control","control"),
libType = c("single-end","paired-end","paired-end","single-end","single-end","paired-end","paired-end"))

sampleTable
suppressPackageStartupMessages(library("DEXSeq"))
dxd = DEXSeqDataSetFromHTSeq(countFiles,sampleData=sampleTable,
design=~sample+exon+condition:exon,flattenedfile=flattenedFile)

genesForSubset = read.table(file.path(inDir,"geneIDsinsubset.txt"),stringsAsFactors=FALSE)[[1]]

dxd =dxd[geneIDs( dxd ) %in% genesForSubset,]
colData(dxd)
head(counts(dxd), 5)
split(seq_len(ncol(dxd)), colData(dxd)$exon)
head(featureCounts(dxd), 5)
head(rowData(dxd), 3)
sampleAnnotation(dxd)
dxd = estimateSizeFactors(dxd)
dxd = estimateDispersions(dxd)
plotDispEsts(dxd)
dxd = testForDEU(dxd)
dxd = estimateExonFoldChanges(dxd, fitExpToVar="condition")
dxr1 = DEXSeqResults(dxd)
dxr1
elementMetadata(dxr1)$description

table ( dxr1$padj 0.1 )
table(tapply(dxr1$padj 0.1, drx1$groupID, any))
plotMA(dxr1, cex=0.8)

[dpryan]

Looks like either you or Alejandro forget a "<". Try:

table(dxr1$padj < 0.1 )

[Jerry_Zhao]

Thanks!

I forgot the "<". This is such an obvious mistake.


Best,
Jerry

[dpryan]

I've long ago lost track of how many times I've made simple mistakes like that.