problem with getting a subset of alignments mapping to a GRange using readGAlignmentsFromBam
1
0
Entering edit mode
rasi1983 ▴ 10
@rasi1983-7454
Last seen 5.0 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.1k views
ADD COMMENT
2
Entering edit mode
@james-w-macdonald-5106
Last seen 2 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
0
Entering edit mode

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

ADD REPLY
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
0
Entering edit mode

Thanks for the clarification. 

ADD REPLY

Login before adding your answer.

Traffic: 483 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6