[BioC] iSeq package\'s peakreg fails

[Lauri Lyly [guest]]

If you just look at peakreg's code and what's in here, it'll be very clear what's happening.

Using ppcut, if there's only one value above the cutoff, sum(id) will be 1, by this code:

    id = NULL
    if (ppcut) {
        id = (pp > cutoff)
    }

This leads to using matrix instead of cbind here:

    if (sum(id) != 1) {
        x = cbind(chrn[id], chrpos[id, 2:3], count[id, 1:2],
            rowID[id])
    }
    else {
        x = matrix(c(chrn[id], chrpos[id, 2:3], count[id, 1:2],
            rowID[id]), nrow = 1)
    }

Now, x has only one row at least in my case.

Which in turn makes the following piece of code throw an exception:

(More precisely the substatement x[, 2] +  x[, 3] )

    xlen = nrow(x)
    if (xlen == 1) {
        br = cbind(x[, 1:3], x[, 6], x[, 6], floor((x[, 2] +
            x[, 3])/2), round(pp[id], 2), count[id, ], sum(count[id,
            ]), round(min(count[id, ])/sum(count[id, ]), 2))
        colnames(br) = c("chr", "gstart", "gend", "rstart", "rend",
            "peakpos", "meanpp", "ct1", "ct2", "ct12", "sym")
        br = as.data.frame(br)
        br[1, 1] = chrpos[br[1, 4], 1]
        return(br)
    }


To demonstrate the code in debugger (invoked by debug(peakreag), then running an actual peakreg call):

Browse[2]> x=cbind(chrn[id], chrpos[id, 2:3], count[id, 1:2], rowID[id])
Browse[2]> x
    chrn[id]  gstart    gend ipct1 ipct2 rowID[id]
389        1 2765786 2765788     3    13       389
Browse[2]> x[,2]+x[,3]
[1] 5531574


Browse[2]> x=(matrix(c(chrn[id], chrpos[id, 2:3], count[id, 1:2], rowID[id]), nrow = 1))
Browse[2]> x[,2]+x[,3]
Error in x[, 2] + x[, 3] : non-numeric argument to binary operator 

This error is obviously very nasty because only by changing the cutoff or by feeding certain data, you can cause any iSeq analysis to fail.

 -- output of sessionInfo(): 


Browse[2]> sessionInfo()
R version 2.14.0 (2011-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=C                 LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] iSeq_1.6.0

loaded via a namespace (and not attached):
[1] tools_2.14.0



--
Sent via the guest posting facility at bioconductor.org.