How to use Rsamtools to extracting reads map into single chromosome from BAM file
1
0
Entering edit mode
21 months ago
ant_genome • 0

Previously, in linux shell, I could extract the reads mapping to chrX like this:

samtools view -h -b test.bam chrX -o test.chrX.bam

How to I get the same result in R using Rsamtools, any one knows?

Thanks

Rsamtools extracting bamfile R • 659 views
ADD COMMENT
1
Entering edit mode
6 months ago
opplatek ▴ 110

Rsamtools manual suggests how to handle big bam files (section 2.1 Large bam files). Basically, you can use which option for in ScanBamParam to limit the loaded regions. From here, you can use something like this to load one whole chromosome at the time:

loadChr <- function(seqname, bamFile)
    {
      bamInfo<-seqinfo(BamFile(bamFile)) # read reference names and their lengths lenthgs
      afrom<-1 # start of a region to extract
      ato<-seqlengths(bamInfo[seqname]) # end of a region to extract 
      param <- ScanBamParam(what=scanBamWhat(),
                              which=GRanges(seqname, IRanges(afrom, ato)),
                              flag=scanBamFlag(isUnmappedQuery=FALSE)) # set parameters for loading
      x <- scanBam(file = bamFile, param=param)[[1]] # load bam with region limits
      return(x)
    }

seqnames <- "chrX"
bam <- lapply(seqnames, loadChr, "in.bam")

I am not aware of a simple way to directly export BAM using Rsamtools, though. You can read BAM header by scanBamHeader and then somehow merge it with the rest, export as SAM and then convert it to BAM with asBam (both from Rsamtools).

ADD COMMENT
0
Entering edit mode

I agree. The lack of bam export functionality in Rsamtools is frustrating and converting from sam to bam is the only way I know

ADD REPLY

Login before adding your answer.

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