This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Convert Samtools flags to R code

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.

r samtools
• 1,791 views • link
updated by Ram Ram • written by DanK • DanK

1 answer

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

• 0 views • link
updated by Ram Ram • written by Istvan Albert Istvan Albert

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?

• 0 views • link
written by DanK • DanK

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)
• 0 views • link
written by Istvan Albert Istvan Albert

You are right! Thank you.

• 0 views • link
written by DanK • DanK

Log in to answer this question.

More posts like this