[Bioc-devel] Some suggestion DEXSeq

[Alejandro Reyes]

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
>