[Bioc-devel] Some suggestion DEXSeq

Hi Steve,

Thanks for your email and advice, I added the next code in newExonCountSet:


rownames(countData) <- paste(geneIDs, exonIDs, sep=":")

if( any( duplicated(rownames(countData) ) ) ) {
   stop("The geneIDs or exonIDs are not unique")
}

And in the validity function:

if( ! all(featureNames(object) == paste(geneIDs(object), exonIDs(object),
sep=":") ) )
   return( "The featureNames do not match with the geneIDs:exonIDs" )
if( ! all(rownames( counts(object) ) == featureNames(object) ) )
   return( "The rownames of the count matrix do not coincide with the
featureNames" )
if( ! all(rownames( fData( object ) ) == featureNames( object ) ) )
   return( "The rownames of the featureData do not coincide with the
featureNames" )

I think this would handle any possible problem!
Thanks again,

Alejandro

Hi,

=== 1: set rownames() on the countData object in `newExonCountSet` ===
I've built my own ExonCountSet by stitching the required pieces
together for a successful call to `newExonCountSet`.

Everything has been fine and dandy up until a call to
`estimatelog2FoldChanges`. At some point during its run, I am greeted
with:

Error in toadd[rownames(alleffects), colnames(alleffects)]<- alleffects :
   subscript out of bounds

So I've been stepping through this function to debug and I realize
that the problem is that my original `countData` matrix I passed into
`newExonCountSet` didn't have its rownames set, and it appears that
w/in the `estimatelog2FoldChanges` function, you expect it to via a
call to:

rownames(toadd)<- featureNames(ecs)

So, would it be a good idea to add some code in `newExonCountSet` to
fix this from the get go? Eg. if `is.null(rownames(countData))`, then
`rownames(countData)<- paste(geneIDs, exonIDs, sep=":")`

And add a check in the validity function for an ExonCountSet that
ensures `all(rownames(countData) == paste(geneIDs, exonIDs, sep=":"))`

=== 2: Guard against user error in exon ID names ===

As I constructed my ExounCountSet, I made a mistake in my `exonIDs`
that took me a while to figure out. I essentially did A Dumb Thing
(tm) and somehow magically had some duplicate exonIDs for the same
geneIDs.

Perhaps another check in ExonCountSet validity function to ensure that:

`all(tapply(eid, gid, function(x) length(x) == length(unique(x))))`

Thanks,
-steve