Question

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

0

Entering edit mode

rasi1983 ▴ 10

@rasi1983-7454

Last seen 6.9 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.6k views

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

score 2 · Answer 1 · 2015-03-12

2

Entering edit mode

James W. MacDonald 67k

@james-w-macdonald-5106

Last seen 4 days 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 9.8 years ago James W. MacDonald 67k

0

Entering edit mode

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

ADD REPLY • link 9.8 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 9.8 years ago Hervé Pagès 16k

0

Entering edit mode

Thanks for the clarification.

ADD REPLY • link 9.8 years ago rasi1983 ▴ 10