Integer overflow when tabulating read coverage
1
0
Entering edit mode
jma1991 ▴ 70
@jma1991-11856
Last seen 21 months ago
Cumbernauld

I would like to create a plot of sequencing coverage. I calculate base-pair resolution coverage of my sequencing data and then run the table function to get the number of base-pairs covered by 1,2,3,4,e.t.c reads. However, I get the following integer overflow error:

> bamCoverage <- coverage(bamFile)
> table(bamCoverage)
Error in .Call2("Rle_constructor", values, lengths, check, 0L, PACKAGE = "S4Vectors") : integer overflow while summing elements in 'lengths'

I also encounter this problem in the solution to a separate question I posted here (https://support.bioconductor.org/p/90695/#90711). Is there a particular way I should approach Rle arithmetic when performing it on a genome-wide scale which would avoid the integer overflow error?

Edit

Tabulation seems to work when I apply the table function on each list element of the coverage object, but I would like to get a whole genome table, rather than a per chromosome table count:

> tableCounts <- lapply(bamCoverage, table)
> tableCounts

$chrM

    0 
16299 

$chrX

        0         1         2         3         4         5         6         7         8         9        10        11        12        13 
158355871   8387790   2316684    819266    384669    220263    136390     92388     66952     50029     36982     29198     22878     18590 
       14        15        16        17        18        19        20        21        22        23        24        25        26        27 
    14481     11929      9701      8351      6964      5768      5056      3987      3591      3249      2671      2211      1958      1723 
       28        29        30        31        32        33        34        35        36        37        38        39        40        41 
     1418      1376      1030      1013       886       628       600       516       421       386       335       249       230       205 
       42        43        44        45        46        47        48        49        50        51        52        53        54        55 
      189       134       137       140        81        94        71        63        55        35        39        18        20        16 
       56        57        58        59        60        61        62        63        64        65        66        67        68        69 
       15        21        10        20        11        33         7         4        14         8        14        11        14        19 
       70        71        72        73        74        75        76        77        78        79        80        81        82        83 
       13        15        13        21         9        15        14         9        13        23        12        15        13         6 
       84        85        86        87        88        89        90        91        92        93        94        95        96        97 
        7         3        10        13        13        14        12        12        12        16        14        12        14        15 
       98        99       100       101       102       103       104       105       106       107       108       109       110       111 
       11        10        12         5         4         7        14         5        13        10        13         3        10         8 
      112       113       114       115       116       117       118       119       120       121       122       123       124       125 
        5        13        13         9        16        18        16        13        26        17        22        18        22        25 
      126       127       128       129       130       131       132       133       134       135       136       137       138       139 
       16         6         6         7        14        15         7        10         6         9         9        10         8         8 
      140       141       142       143       144       145       146       147       148       149       150       151       152       153 
        4         3         4         6        10         8         9         7        10        14         7         3         6         4 
      154       155       156       157       158       159       160       161       162       163       164       165       166       167 
        6         4         1         6         4         6        11         4         7         9        10        15        15         9 
      168       169       170       171       172       174       176       178       179       181       182       183       186       195 
        6         1         6         2         6         2         3         1         1        23        18         3         1         1 

$chrY

       0        1        2        3        4        5        6        7        8        9       10       11       12       13       14 
91527964   170721    27817     6670     2521     1409      875      442      294      271      305      213      193      225      161 
      15       16       17       18       19       20       21       22       23       24       25       26       27       28       29 
     116      163      157      139      163      115      161       97      104      108      100       93       96      112      179 
      30       31       32       33       34       35       36       37       38       39       40       41       42       43       44 
     120      132      120       89       57       85       87       68       50      103       55       64       59       39       42 
      45       46       47       48       49       50       51       52       53       54       55       56       57       58       59 
      61       49       44       53       44       29       30       31       48       67       29       21       24       29       32 
      60       61       62       63       64       65       66       67       68       69       70       71       72       73       74 
       9       12       10       15       16       13       15       23       38       33       23       14       30       47       29 
      75       76       77       78       79       80       81       82       83       84       85       86       87       88       89 
      43       19       13        7       13       37        6       11       11        7       18       49       10       10       18 
      90       91       92       93       94       95       96       97       98      100      101      102      103      104      105 
       6        8       15       16       10        7       13       22       13       14       10       16        6        9       15 
     106      107      108      109      110      111      112      113      114      115      116      117      119      120      121 
       3        8       13       10        7        5        3        3        6        5        6        1        3        2        1 
     122      124      126      127      128      129      130      132      133      134      135      136      137      138      140 
       1        2        2        2        2        5        1        3        5        4        1        7        3        2        1 
     141      143      145      146      148      149      150      152      154      155      157      160      161      163      164 
       2        1        2        1        2        1        2        1        1        3        5        1        4        2        3 
     165      166      167      168      171      172      174      176      179      180      182      183      184 
       6        3        2        1        2        2        1        2        3        1        1        1        1 

 

 

 

genomicalignments genomicranges • 1.4k views
ADD COMMENT
2
Entering edit mode
@michael-lawrence-3846
Last seen 3.0 years ago
United States

I made this function:

setMethod("table", "SimpleAtomicList", function(...)
{
    args <- list(...)
    if (length(args) != 1L)
        stop("\"table\" method for SimpleAtomicList objects ",
             "can only take one input object")
    x <- args[[1L]]
    levs <- sort(unique(unlist(lapply(x, function(xi) {
        if (!is.null(levels(xi))) levels(xi) else unique(xi)
    }), use.names=FALSE)))
    as.table(do.call(rbind,
                     lapply(x, function(xi) {
                         if (is(xi, "Rle"))
                             runValue(xi) <- factor(runValue(xi), levs)
                         else xi <- factor(xi, levs)
                         table(xi)
                     })))
})

It will be added to trunk soon. With that, you can compute a 2-way table without entering into long vector territory, and then take the margin for each column:

margin.table(table(bamCoverage), 2)

 

 

ADD COMMENT
0
Entering edit mode

Thank you, that works perfectly

ADD REPLY

Login before adding your answer.

Traffic: 554 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6