Question: problem with getting a subset of alignments mapping to a GRange using readGAlignmentsFromBam
0
gravatar for rasi1983
4.7 years ago by
rasi198310
United States
rasi198310 wrote:

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
 
 
ADD COMMENTlink modified 4.7 years ago by James W. MacDonald52k • written 4.7 years ago by rasi198310
Answer: problem with getting a subset of alignments mapping to a GRange using readGAlign
2
gravatar for James W. MacDonald
4.7 years ago by
United States
James W. MacDonald52k wrote:

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 COMMENTlink written 4.7 years ago by James W. MacDonald52k

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

ADD REPLYlink written 4.7 years ago by rasi198310
1

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 REPLYlink written 4.7 years ago by Hervé Pagès ♦♦ 14k

Thanks for the clarification. 

ADD REPLYlink written 4.7 years ago by rasi198310
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 291 users visited in the last hour