I would like to find out where (in genomic coordinates) large insertions are found within a given BAM file,
(In terms of genomic coordinates, I just mean I would like to have a rough idea where to look at the BAM in IGV, chromosome and position.)
Based on the BAM CIGAR, I should be able to
grep the BAM based on some cut-off of "large inserted" bases.
It's unclear to me how to best do this via samtools; one could parse the cigar string via
samtools for insertions, but I don't understand how to translate this into the genomic position in the alignment.
This goes with various samtools wrappers, e.g.
pysam in python:
import pysam bamfile = pysam.Samfile(bamFile, "rb"); for read in bamfile: if not read.is_unmapped: ## only get mapped reads cigarLine = read.cigar for cigarType, cigarLength in cigarLine: if cigarType == 1: ## this means it's an insertion... ## somehow filter only large insertions, and then get coordinate information
I'm not sure how to do this in a non-clumsy manner. Any help appreciated.