Convert Samtools flags to R code

Question

[DanK]

Hi,

I have the commands "samtools view -f 64" and I want to convert into R code, with the Rsamtools package. Does anyone know how this is done?

I'm confused reading the Rsamtools overview.

Thank you in advance.

Answer 1

[Istvan Albert]

It will look like this, for example to select alignments on forward strand do:

library(Rsamtools)  

flag = scanBamFlag(isMinusStrand=FALSE)
param = ScanBamParam(what=scanBamWhat(), flag=flag)
aln = scanBam('alignment.bam', param=param)

see more information in the scanBamParam documentation ?scanBamParam

[DanK]

Thank you very much for the reply.

Your example is very instructive. I think that the -f 64 flag is to extract the first read in pair-end sequencing, so the R flag will be the isFirstMateRead=TRUE.

Right?

[Istvan Albert]

I believe so, I will say that you should double-check the counts as well with countBam() instead of scanBam(). This latter will return a count that you can verify from command line as well:

samtools view -cf 64 alignment.bam

ought to give you the same count as:

flag = scanBamFlag(isFirstMateRead=TRUE)
param = ScanBamParam(what=scanBamWhat(), flag=flag)
count = countBam('alignment.bam', param=param)
print(count$records)

[DanK]

You are right! Thank you.