[Bioc-devel] Why be default 'D' is not dropped from coverage? - Bioc-devel

Sat, Oct 18, 2014 10:50 AM #

Hi,

A collaborator of mine is working on a new software and we while we were
doing a sanity check to compare the base-level coverage from BAM files and
bigWig files generated from his software we realized that by default bases
corresponding to a 'D' on the CIGAR string get counted when reading
coverage from BAM files.

That is:

Method Definition:

function (x, shift = 0L, width = NULL, weight = 1L, ...)
{
    .local <- function (x, shift = 0L, width = NULL, weight = 1L,
        method = c("auto", "sort", "hash"), drop.D.ranges = FALSE)
    coverage(grglist(x, drop.D.ranges = drop.D.ranges), shift = shift,
        width = width, weight = weight, method = method)
    .local(x, shift, width, weight, ...)
}
<environment: namespace:GenomicAlignments>

Signatures:
        x
target  "GAlignments"
defined "GAlignments"

[1] ?1.2.0?


You can see that this is default elsewhere, for example:

function (cigar, pos = 1L, drop.D.ranges = FALSE, f = NULL)
{
    if (!isTRUEorFALSE(drop.D.ranges))
        stop("'drop.D.ranges' must be TRUE or FALSE")
    if (drop.D.ranges) {
        ops <- c("M", "=", "X", "I")
    }
    else {
        ops <- c("M", "=", "X", "I", "D")
    }
    cigarRangesAlongReferenceSpace(cigar, flag = NULL, pos = pos,
        f = f, ops = ops, drop.empty.ranges = FALSE, reduce.ranges = TRUE)
}
<environment: namespace:GenomicAlignments>


What is the rationale for setting `drop.D.ranges` by default to FALSE?

Thanks,
Leo

Hervé Pagès

Mon, Oct 20, 2014 2:45 PM #

Hi Leo,

On 10/18/2014 10:50 AM, Leonardo Collado Torres wrote:

Because last time I checked (but that was 3 or 4 years ago), that's what
Samtools pileup was doing.

Cheers,
H.

Herv? Pag?s

Program in Computational Biology
Division of Public Health Sciences
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N, M1-B514
P.O. Box 19024
Seattle, WA 98109-1024

E-mail: hpages at fhcrc.org
Phone:  (206) 667-5791
Fax:    (206) 667-1319

Leonardo Collado Torres

Tue, Oct 21, 2014 8:10 PM #

Thank you Herv?! I failed to realized that this was included in the docs.

Anyhow, we've been discussing internally what default value to use. We
also noted that the `genomecov` tool from BED in it's split mode
ignores the D's:
http://bedtools.readthedocs.org/en/latest/content/tools/genomecov.html

Cheers,
Leo

On Mon, Oct 20, 2014 at 5:45 PM, Herv? Pag?s <hpages at fhcrc.org> wrote:

Hi Leo,


On 10/18/2014 10:50 AM, Leonardo Collado Torres wrote:

Hi,

A collaborator of mine is working on a new software and we while we were
doing a sanity check to compare the base-level coverage from BAM files
and bigWig files generated from his software we realized that by default
bases corresponding to a 'D' on the CIGAR string get counted when
reading coverage from BAM files.

That is:

 > getMethod('coverage', 'GAlignments')

Method Definition:

function (x, shift = 0L, width = NULL, weight = 1L, ...)
{
     .local <- function (x, shift = 0L, width = NULL, weight = 1L,
         method = c("auto", "sort", "hash"), drop.D.ranges = FALSE)
     coverage(grglist(x, drop.D.ranges = drop.D.ranges), shift = shift,
         width = width, weight = weight, method = method)
     .local(x, shift, width, weight, ...)
}
<environment: namespace:GenomicAlignments>

Signatures:
         x
target  "GAlignments"
defined "GAlignments"

 > packageVersion('GenomicAlignments')

[1] ?1.2.0?


You can see that this is default elsewhere, for example:

 > extractAlignmentRangesOnReference

function (cigar, pos = 1L, drop.D.ranges = FALSE, f = NULL)
{
     if (!isTRUEorFALSE(drop.D.ranges))
         stop("'drop.D.ranges' must be TRUE or FALSE")
     if (drop.D.ranges) {
         ops <- c("M", "=", "X", "I")
     }
     else {
         ops <- c("M", "=", "X", "I", "D")
     }
     cigarRangesAlongReferenceSpace(cigar, flag = NULL, pos = pos,
         f = f, ops = ops, drop.empty.ranges = FALSE, reduce.ranges =
TRUE)
}
<environment: namespace:GenomicAlignments>


What is the rationale for setting `drop.D.ranges` by default to FALSE?


Because last time I checked (but that was 3 or 4 years ago), that's what
Samtools pileup was doing.

Cheers,
H.

Thanks,
Leo

--
Herv? Pag?s

Program in Computational Biology
Division of Public Health Sciences
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N, M1-B514
P.O. Box 19024
Seattle, WA 98109-1024

E-mail: hpages at fhcrc.org
Phone:  (206) 667-5791
Fax:    (206) 667-1319