Question

problem with getting a subset of alignments mapping to a GRange using readGAlignmentsFromBam

0

Entering edit mode

rasi1983 ▴ 10

@rasi1983-7454

Last seen 7.7 years ago

United States

I am unable to read a subset of alignments from a Bamfile as specified by a GRange. I am new to Bioconductor - am I missing something obvious?

779 Alignments in the Bamfile

> readGAlignmentsFromBam(Bamfile)
GAlignments object with 779 alignments and 0 metadata columns:
        seqnames strand       cigar    qwidth     start       end     width     njunc
           <Rle>  <Rle> <character> <integer> <integer> <integer> <integer> <integer>
    [1]       12      -         24M        24  25359522  25359545        24         0
    [2]       12      -         30M        30  25359527  25359556        30         0
    [3]       12      -         28M        28  25359755  25359782        28         0
    [4]       12      -         27M        27  25362687  25362713        27         0
    [5]       12      -         28M        28  25362734  25362761        28         0
    ...      ...    ...         ...       ...       ...       ...       ...       ...
  [775]       12      - 22M5355N12M        34  25398308  25403696      5389         1
  [776]       12      -         31M        31  25403736  25403766        31         0
  [777]       12      -         25M        25  25403737  25403761        25         0
  [778]       12      -         38M        38  25403745  25403782        38         0
  [779]       12      -         24M        24  25403745  25403768        24         0
  -------
  seqinfo: 25 sequences from an unspecified genome

The my_grange object overlaps with only a fraction of the above alignments (92/779)

> my_grange
GRanges object with 1 range and 1 metadata column:
      seqnames               ranges strand |         tx_id
         <Rle>            <IRanges>  <Rle> | <IntegerList>
  [1]       12 [25358180, 25362845]      - |   47893,47894
  -------
  seqinfo: 1 sequence from an unspecified genome

Using the which param after subseting with my_grange still returns all the 779 alignments

> readGAlignmentsFromBam(Bamfile, which = ScanBamParam(which = my_grange))
GAlignments object with 779 alignments and 0 metadata columns:
        seqnames strand       cigar    qwidth     start       end     width     njunc
           <Rle>  <Rle> <character> <integer> <integer> <integer> <integer> <integer>
    [1]       12      -         24M        24  25359522  25359545        24         0
    [2]       12      -         30M        30  25359527  25359556        30         0
    [3]       12      -         28M        28  25359755  25359782        28         0
    [4]       12      -         27M        27  25362687  25362713        27         0
    [5]       12      -         28M        28  25362734  25362761        28         0
    ...      ...    ...         ...       ...       ...       ...       ...       ...
  [775]       12      - 22M5355N12M        34  25398308  25403696      5389         1
  [776]       12      -         31M        31  25403736  25403766        31         0
  [777]       12      -         25M        25  25403737  25403761        25         0
  [778]       12      -         38M        38  25403745  25403782        38         0
  [779]       12      -         24M        24  25403745  25403768        24         0
  -------
  seqinfo: 25 sequences from an unspecified genome

genomicalignments scanbamparam • 1.8k views

ADD COMMENT • link updated 10.6 years ago by James W. MacDonald 68k • written 10.6 years ago by rasi1983 ▴ 10

score 2 · Answer 1 · 2015-03-12

2

Entering edit mode

James W. MacDonald 68k

@james-w-macdonald-5106

Last seen 8 hours ago

United States

You want to pass the GRanges in using the 'param' argument, not the (non-existant) 'which' argument.

readGAlignmentsFromBam(Bamfile, param = ScanBamParam(which = my_grange))

ADD COMMENT • link 10.6 years ago James W. MacDonald 68k

0

Entering edit mode

Thanks, I missed an `R` gotcha: you can pass nonexistent parameters and not get any errors!

ADD REPLY • link 10.6 years ago rasi1983 ▴ 10

1

Entering edit mode

Well, that's because readGAlignmentsFromBam() has the ... argument. In BioC devel, readGAlignmentsFromBam() is deprecated in favor of readGAlignments(). And readGAlignments() doesn't have the ... argument anymore so you'll get an error if you try to pass a non existent argument to it. So it's not R but the way we design our functions and how we handle the ...

H.

ADD REPLY • link 10.6 years ago Hervé Pagès 16k

0

Entering edit mode

Thanks for the clarification.

ADD REPLY • link 10.6 years ago rasi1983 ▴ 10