How to use Rsamtools to extracting reads map into single chromosome from BAM file

Question

[ant_genome]

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

Answer 1

[opplatek]

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).

[changxu.fan]

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