This is a test version of Biostars. For the public version, visit https://www.biostars.org.
R: Readaligned Only Junction Reads From Bam-File

Hi,

I want to read in a BAM-file into R using the readAligned function from the ShortRead-package. However, I am just interested in the junction reads without indels and would like to filter using the CIGAR string for reads just containing M & Ns. How is that possible? I am just aware of the ScanBam argument "simpleCigar", which is not sufficient.

library("ShortRead")
param <- ScanBamParam(simpleCigar=F)

.. reads in fully aligned reads and junction reads.. I would need a filter for reading in junction reads without indels

readAligned(".", pattern=.bam,type="BAM",param=param)

Thanks for helping me out!!

r bam cigar

2 answers

You might take a look at readGappedAlignments in the GenomicRanges package. Once you read in your sequences, call ngap() on the resulting object to get the number of gaps per read.

In general if you want to filter SAM files the fastest and most efficient way is to filter from command line. If you want to avoid insertions/deletions you should filter for not having I and D in the CIGAR string like so:

# create the headers
egrep '^@' data.sam > small.sam 
# filter the CIGAR string   
cat data.sam | awk ' $6 !~ /(I|D)/ { print $0 } ' >> small.sam

Log in to answer this question.