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
•
written
by
ant_genome •
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).
I agree. The lack of bam export functionality in Rsamtools is frustrating and converting from sam to bam is the only way I know
Log in to answer this question.
