Question: Count reads with package "bamsignals" from a female bam file (no chrY)
1
gravatar for steige_e
4.5 years ago by
steige_e10
Germany
steige_e10 wrote:

Hi!

I got my gene annotation as a GRanges-object from Biomart. Now I try to count the reads from a bamfile using this annotation with the bamCount function from the bamsignals package. Unfortunately this otherwise excellent function gives me an

Error: chromosome chrY not present in the bam file

(no surprise, since the data comes from a woman). Since I look at more bamfiles which are from men and women, I'd like something like NA for the chrY gene counts. What would be the best way to tackle this problem? I'm sure I'm not the first one to encounter a woman.

Thanks, Edgar

genomicfeatures bamsignals • 744 views
ADD COMMENTlink modified 4.5 years ago by mammana20 • written 4.5 years ago by steige_e10
Answer: Count reads with package "bamsignals" from a female bam file (no chrY)
2
gravatar for Martin Morgan
4.5 years ago by
Martin Morgan ♦♦ 24k
United States
Martin Morgan ♦♦ 24k wrote:

Usually reads are aligned to reference genomes, and reference genomes contain chry, so probably your up-stream analysis has done something different from usual and excluded chrY for this sample.

Here's a bam file, and a way to discover the 'seqlevels' (chromsomes) it knows about

> fl = system.file("extdata", "ex1.bam", package="Rsamtools")
> seqlevels(BamFile(fl))
[1] "seq1" "seq2"

Here's a GRanges with another chromosome

> gr
GRanges object with 3 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     seq1   [1, 10]      *
  [2]     seq2   [1, 10]      *
  [3]     seq3   [1, 10]      *
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

And finally an update to the GRanges to keep only those seqlevels present in the BAM file

> keepSeqlevels(gr, seqlevels(BamFile(fl)))
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     seq1   [1, 10]      *
  [2]     seq2   [1, 10]      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

 

 

ADD COMMENTlink written 4.5 years ago by Martin Morgan ♦♦ 24k
Answer: Count reads with package "bamsignals" from a female bam file (no chrY)
1
gravatar for mammana
4.5 years ago by
mammana20
European Union
mammana20 wrote:

Thanks for the post. To insert NAs for those GRanges in chromosomes that are not in the header of the bam file is a solution I already thought about. It is probably not a very clean solution, but on the other hand changing the header of bam files is much more time consuming. So I will add this feature in the development version of bamsignals.

To comment on the already excellent answer given by Martin, and to illustrate the future behaviour of bamsignals, this is a workaround:

>gr
GRanges object with 3 ranges and 0 metadata columns:
      seqnames       ranges strand
         <Rle>    <IRanges>  <Rle>
  [1]     chr3 [5000, 5009]      *
  [2]     chr3 [6000, 6009]      *
  [3]     chr4 [1000, 1009]      *

> seqlevels(BamFile(bampath))
[1] "chr1" "chr2" "chr3"
> counts <- rep(NA, length(gr))
> valid <- as.logical(seqnames(gr) %in% seqlevels(BamFile(bampath)))
> counts[valid] <- bamCount(bampath, gr[valid])

 

Thanks for the suggestion!

Best,

Ale

ADD COMMENTlink modified 4.5 years ago • written 4.5 years ago by mammana20
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: 166 users visited in the last hour