{"action":"created","issue":{"url":"https://api.github.com/repos/ComputationalProteomicsUnit/Pbase/issues/19","labels_url":"https://api.github.com/repos/ComputationalProteomicsUnit/Pbase/issues/19/labels{/name}","comments_url":"https://api.github.com/repos/ComputationalProteomicsUnit/Pbase/issues/19/comments","events_url":"https://api.github.com/repos/ComputationalProteomicsUnit/Pbase/issues/19/events","html_url":"https://github.com/ComputationalProteomicsUnit/Pbase/issues/19","id":58423381,"number":19,"title":"Protein coverage","user":{"login":"lgatto","id":384198,"avatar_url":"https://avatars.githubusercontent.com/u/384198?v=3","gravatar_id":"","url":"https://api.github.com/users/lgatto","html_url":"https://github.com/lgatto","followers_url":"https://api.github.com/users/lgatto/followers","following_url":"https://api.github.com/users/lgatto/following{/other_user}","gists_url":"https://api.github.com/users/lgatto/gists{/gist_id}","starred_url":"https://api.github.com/users/lgatto/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/lgatto/subscriptions","organizations_url":"https://api.github.com/users/lgatto/orgs","repos_url":"https://api.github.com/users/lgatto/repos","events_url":"https://api.github.com/users/lgatto/events{/privacy}","received_events_url":"https://api.github.com/users/lgatto/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":1,"created_at":"2015-02-20T23:08:39Z","updated_at":"2015-02-27T16:53:10Z","closed_at":null,"body":"I have some issues with `proteinCoverage`. I find the implementation a bit confusion and can't figure out why I end up with `NA`.\r\n\r\n```r\r\n> suppressPackageStartupMessages(library(\"Pbase\"))\r\n> data(p)\r\n> p <- proteinCoverage(p)\r\n> acols(p)$Coverage\r\n    A4UGR9     A6H8Y1     O43707     O75369     P00558     P02545     P04075 \r\n0.10758743 0.10518293 0.08232711         NA 0.16067146 0.18373494 0.59615385 \r\n  P04075-2     P60709 \r\n0.49043062         NA \r\n```\r\nbut `O75369` looks fine to me\r\n```r\r\nplot(p[4])\r\n```\r\n![rplot001](https://cloud.githubusercontent.com/assets/384198/6310992/85d2280e-b954-11e4-8888-551cde81ac8e.png)\r\n\r\nA naive implementation gives the following, which matches `proteinCoverage` ommitting `NA`s\r\n\r\n```r\r\n> prtl <- width(aa(p))\r\n> pepl <- sapply(pranges(p), function(x) sum(width(reduce(x))))\r\n> cvg <- pepl/prtl\r\n> rbind(proteinCoverage = acols(p)$Coverage, cvg)\r\n                   A4UGR9    A6H8Y1     O43707     O75369    P00558    P02545\r\nproteinCoverage 0.1075874 0.1051829 0.08232711         NA 0.1606715 0.1837349\r\ncvg             0.1075874 0.1051829 0.08232711 0.06379708 0.1606715 0.1837349\r\n                   P04075  P04075-2     P60709\r\nproteinCoverage 0.5961538 0.4904306         NA\r\ncvg             0.5961538 0.4904306 0.03466667\r\n```\r\n@sgibb, could you clarify what happens and why?"},"comment":{"url":"https://api.github.com/repos/ComputationalProteomicsUnit/Pbase/issues/comments/76427768","html_url":"https://github.com/ComputationalProteomicsUnit/Pbase/issues/19#issuecomment-76427768","issue_url":"https://api.github.com/repos/ComputationalProteomicsUnit/Pbase/issues/19","id":76427768,"user":{"login":"sgibb","id":1828443,"avatar_url":"https://avatars.githubusercontent.com/u/1828443?v=3","gravatar_id":"","url":"https://api.github.com/users/sgibb","html_url":"https://github.com/sgibb","followers_url":"https://api.github.com/users/sgibb/followers","following_url":"https://api.github.com/users/sgibb/following{/other_user}","gists_url":"https://api.github.com/users/sgibb/gists{/gist_id}","starred_url":"https://api.github.com/users/sgibb/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/sgibb/subscriptions","organizations_url":"https://api.github.com/users/sgibb/orgs","repos_url":"https://api.github.com/users/sgibb/repos","events_url":"https://api.github.com/users/sgibb/events{/privacy}","received_events_url":"https://api.github.com/users/sgibb/received_events","type":"User","site_admin":false},"created_at":"2015-02-27T16:53:10Z","updated_at":"2015-02-27T16:53:10Z","body":"You are right, it is confusing and contains at least one bug.\r\n\r\nThe first lines check the arguments and ensure that the order of the `pattern` (peptides) is the same as in the `subject` (proteins; mabye we can skip the ordering step? or implement a general function?): [ll. 19-41](https://github.com/ComputationalProteomicsUnit/Pbase/blob/master/R/proteinCoverage-functions.R#L19-L41)\r\n\r\nThe crucial step is in [ll. 43-46] (https://github.com/ComputationalProteomicsUnit/Pbase/blob/master/R/proteinCoverage-functions.R#L43-L46). Here I check whether the peptides is within a protein (I can't remember why I have implemented this check. It should ensure that the peptides sequence is not longer than the protein sequence. Could this really happen? Maybe this is completley superfluous.)\r\nWhile this seems useless it reveals a bug in `.splitIRanges` in combination with the argument `unshift = TRUE`. `unshift = TRUE` assumes an ordered `IRangesList` (which is not the case in our example). I also can't remember why I implemented the `unshift = TRUE` option. By using `unshift = TRUE` we get some `IRanges` that start at one, e.g.:\r\n```r\r\nir <- IRanges(start = c(1, 6), end = c(5, 10))\r\nir\r\n# IRanges of length 2\r\n#     start end width\r\n# [1]     1   5     5\r\n# [2]     6  10     5\r\nPbase:::.splitIRanges(ir, f = c(\"a\", \"b\"), unshift = FALSE)\r\n# IRangesList of length 2\r\n# $a\r\n# IRanges of length 1\r\n#     start end width\r\n# [1]     1   5     5\r\n#\r\n# $b\r\n# IRanges of length 1\r\n#     start end width\r\n# [1]     6  10     5\r\n\r\nPbase:::.splitIRanges(ir, f = c(\"a\", \"b\"), unshift = TRUE)\r\n# IRangesList of length 2\r\n# $a\r\n# IRanges of length 1\r\n#     start end width\r\n# [1]     1   5     5\r\n#\r\n# $b\r\n# IRanges of length 1\r\n#     start end width\r\n# [1]     1   5     5\r\n```\r\n\r\nThe last lines [ll. 48-51](https://github.com/ComputationalProteomicsUnit/Pbase/blob/master/R/proteinCoverage-functions.R#L48-L51) are exactly the same as the lines you suggested.\r\n\r\nAFAIK we don't use the `unshift` option anywhere. If we avoid the check whether a peptide is in its corresponding protein we can remove the `unshift` option completely and make the code more readable. Otherwise I need to fix `.splitIRanges`. What do you think? I tend to remove the `unshift` argument."}}