← Back to branch summary

~ubuntu-branches/ubuntu/wily/r-bioc-genomicranges/wily-proposed

« back to all changes in this revision

Viewing changes to R/coverage-methods.R

  • Committer: Package Import Robot
  • Author(s): Andreas Tille
  • Date: 2013-10-18 10:40:04 UTC
  • Revision ID: package-import@ubuntu.com-20131018104004-ktm4ub0pcoybnir6
Tags: upstream-1.12.4
Import upstream version 1.12.4

Show diffs side-by-side

added added

removed removed

Lines of Context:
R/coverage-methods.R
 
1
### =========================================================================
 
2
### coverage methods
 
3
### -------------------------------------------------------------------------
 
4
###
 
5
 
 
6
### Should work if 'arg' is an ordinary vector (i.e. atomic vector or list),
 
7
### or a List object.
 
8
.coverage.recycleAndSetNames <- function(arg, argname, seqlevels)
 
9
{
 
10
    k <- length(seqlevels)
 
11
    if (length(arg) < k)
 
12
        arg <- rep(arg, length.out = k)
 
13
    if (is.null(names(arg)))
 
14
        names(arg) <- seqlevels
 
15
    if (!all(seqlevels %in% names(arg)))
 
16
        stop("some seqnames missing from names(", argname, ")")
 
17
    arg
 
18
}
 
19
 
 
20
.coverage.normargShiftOrWeight <- function(arg, argname, x)
 
21
{
 
22
    if (is.list(arg) || is(arg, "List"))
 
23
        return(.coverage.recycleAndSetNames(arg, argname, seqlevels(x)))
 
24
    if (is.numeric(arg)) {
 
25
        if (length(arg) == 1L) {
 
26
            arg <- recycleNumericArg(arg, argname, length(seqlevels(x)))
 
27
            arg <- as.list(arg)
 
28
            names(arg) <- seqlevels(x)
 
29
        } else {
 
30
            arg <- recycleNumericArg(arg, argname, length(x))
 
31
            arg <- split(arg, as.factor(seqnames(x)))
 
32
        }
 
33
        return(arg)
 
34
    }
 
35
    stop("'", argname, "' must be a numeric vector, or a list, ",
 
36
         "or a List object")
 
37
}
 
38
 
 
39
.coverage.normargWidth <- function(width, x)
 
40
{
 
41
    if (is.null(width))
 
42
        width <- seqlengths(x)
 
43
    else if (!is.numeric(width))
 
44
        stop("'width' must be NULL or a numeric vector")
 
45
    .coverage.recycleAndSetNames(width, "width", seqlevels(x))
 
46
}
 
47
 
 
48
### 'circle.length' must be NA (if the underlying sequence is linear) or the
 
49
### length of the underlying circular sequence (integer vector of length 1
 
50
### with the name of the sequence).
 
51
### 'rg' must be an IRanges object.
 
52
.coverage.circle <- function(circle.length, rg, shift, width, weight)
 
53
{
 
54
    if (is.na(circle.length))
 
55
        return(coverage(rg, shift=shift, width=width, weight=weight))
 
56
    cvg <- fold(coverage(rg, weight=weight),
 
57
                circle.length, from=1L-shift)
 
58
    if (is.null(width))
 
59
        return(cvg)
 
60
    if (width > length(cvg))
 
61
        stop("invalid width (", width, ") ",
 
62
             "for circular sequence ", names(circle.length))
 
63
    cvg[seq_len(width)]
 
64
}
 
65
 
 
66
.coverage.seq <- function(seqlength, isCircular, rg, shift, width, weight)
 
67
{
 
68
    ## A shortcut for a common situation. Is it worth it?
 
69
    if (length(rg) == 0L && is.na(width))
 
70
        return(Rle(weight, 0L))
 
71
    if (isCircular %in% TRUE)
 
72
        circle.length <- seqlength
 
73
    else
 
74
        circle.length <- NA
 
75
    if (is.na(width))
 
76
        width <- NULL
 
77
    .coverage.circle(circle.length, rg, shift, width, weight)
 
78
}
 
79
 
 
80
setMethod("coverage", "GenomicRanges",
 
81
    function(x, shift=0L, width=NULL, weight=1L, ...)
 
82
    {
 
83
        if (any(start(x) < 1L))
 
84
            stop("'x' contains ranges starting before position 1. ",
 
85
                 "coverage() currently doesn't support this.")
 
86
        shift <- .coverage.normargShiftOrWeight(shift, "shift", x)
 
87
        width <- .coverage.normargWidth(width, x)
 
88
        if (isSingleString(weight)) {
 
89
            x_mcols <- mcols(x)
 
90
            if (!is(x_mcols, "DataTable")
 
91
             || sum(colnames(x_mcols) == weight) != 1L)
 
92
                stop("'mcols(x)' has 0 or more than 1 \"",
 
93
                     weight, "\" columns")
 
94
            weight <- x_mcols[[weight]]
 
95
        }
 
96
        weight <- .coverage.normargShiftOrWeight(weight, "weight", x)
 
97
        x_seqlevels <- seqlevels(x)
 
98
        x_seqlengths <- seqlengths(x)
 
99
        x_isCircular <- isCircular(x)
 
100
        x_ranges_by_seq <- split(unname(ranges(x)), seqnames(x))
 
101
        ii <- seq_len(length(x_seqlevels))
 
102
        names(ii) <- x_seqlevels
 
103
        ans_listData <- lapply(ii, function(i)
 
104
                               .coverage.seq(x_seqlengths[[i]],
 
105
                                             x_isCircular[[i]],
 
106
                                             x_ranges_by_seq[[i]],
 
107
                                             shift[[i]],
 
108
                                             width[[i]],
 
109
                                             weight[[i]]))
 
110
        IRanges:::newList("SimpleRleList", ans_listData)
 
111
    }
 
112
)
 
113
 
 
114
setMethod("coverage", "GRangesList",
 
115
    function(x, shift=0L, width=NULL, weight=1L, ...)
 
116
    {
 
117
        coverage(x@unlistData, shift=shift, width=width, weight=weight, ...)
 
118
    }
 
119
)
 
120
 
 
121
setMethod("coverage", "GappedAlignments",
 
122
    function(x, shift=0L, width=NULL, weight=1L, drop.D.ranges = FALSE, ...)
 
123
        coverage(grglist(x, drop.D.ranges = drop.D.ranges),
 
124
                 shift=shift, width=width, weight=weight, ...)
 
125
)
 
126
 
 
127
setMethod("coverage", "GappedAlignmentPairs",
 
128
    function(x, shift=0L, width=NULL, weight=1L, drop.D.ranges = FALSE, ...)
 
129
          coverage(grglist(x, drop.D.ranges = drop.D.ranges),
 
130
                   shift=shift, width=width, weight=weight, ...)
 
131
)
 
132