Graphing coverage() results
2
0
Entering edit mode
gen ▴ 30
@gen-7383
Last seen 9.6 years ago
United States

Hi I have a GAlignments Object from scanning 1000-1100 of chromosome x of a bam file

GAlignments object with 46 alignments and 5 metadata columns:
       seqnames strand       cigar    qwidth     start       end     width     njunc   |    rname
          <Rle>  <Rle> <character> <integer> <integer> <integer> <integer> <integer>   | <factor>
   [1]     seq1      +         35M        35       970      1004        35         0   |     seq1
   [2]     seq1      +         35M        35       971      1005        35         0   |     seq1

 

Running coverage() on the object I get an RleList that looks something like this

 

RleList of length 2
$seq1
integer-Rle of length 1575 with 52 runs
  Lengths: 969   1   1   1   1   1   1   5   1   1   1 ...   1   1   1   1   2   1   1   3   1 531
  Values :   0   1   2   3   4   5   6   7  10  12  14 ...  15  12  11  10   8   7   5   4   2   0

$seq2
integer-Rle of length 1584 with 1 run
  Lengths: 1584
  Values :    0

 

I want to graph the coverage data (of the chunk I scanned) with ggbio or ggplot2 but I'm not sure exactly what the coverage result represents. Is there a way to get a graph with chromosome location I've scanned on the x axis and read coverage on the y with this? Or are the results usually graphed a different way?

 

Thanks

ggbio genomicalignments bam • 1.6k views
ADD COMMENT
0
Entering edit mode
@james-w-macdonald-5106
Last seen 19 minutes ago
United States
Have you seen the help for stat_coverage() in ggbio? I believe that should do exactly what you want.
ADD COMMENT
0
Entering edit mode
gen ▴ 30
@gen-7383
Last seen 9.6 years ago
United States

Hello, thanks for your help. I haven't been able to get stat_coverage to fully graph the data

I can convert my GAlignments Object to a GRangesList but not a GRanges

bamranges<-as(bamaligns,"GRanges")
Error in `elementMetadata<-`(x, ..., value = value) : 
  names of metadata columns cannot be one of  "seqnames", "ranges", "strand", "seqlevels", "seqlengths", "isCircular", "start", "end", "width", "element"

 

 

bamranges

GRangesList object of length 46:
[[1]] 
GRanges object with 1 range and 0 metadata columns:
      seqnames      ranges strand
         <Rle>   <IRanges>  <Rle>
  [1]     seq1 [970, 1004]      +

[[2]] 
GRanges object with 1 range and 0 metadata columns:
      seqnames      ranges strand
  [1]     seq1 [971, 1005]      +

 

 

 

when I attempt to use the GRangesList in stat_coverage I also get an error.

> ggplot()+stat_coverage(bamranges)
Error in `elementMetadata<-`(`*tmp*`, value = <S4 object of class "DataFrame">) : 
  names of metadata columns cannot be one of  "seqnames", "ranges", "strand", "seqlevels", "seqlengths", "isCircular", "start", "end", "width", "element"

I can graph one element of the GRangesList but that is not particularly useful

> ggplot()+stat_coverage(bamranges[[1]])

 

 

 

ADD COMMENT
0
Entering edit mode

Hi,

It looks like the coercion methods (to GRanges or GRangeList) choke on some of the metadata columns of your GAlignments object. Try to remove them first with:

mcols(bamaligns) <- NULL

or coerce to GRanges or GRangesList with:

granges(bamaligns)
grglist(bamaligns)

Hope this will help,

H.

ADD REPLY
0
Entering edit mode
How about using the autopilot() function, with the which argument to define the region or regions you care about?
ADD REPLY
0
Entering edit mode
The autoplot function can use the bam file directly as input.
ADD REPLY

Login before adding your answer.

Traffic: 950 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