# Parse BAM by insertion size and get genomic coordinates

I would like to find out where (in genomic coordinates) large insertions are found within a given BAM file, file.bam.

(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 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.

You basically have the right framework already in your pysam code, which I took the liberty of tidying up a bit. You just need to add a bit:

for read in bamfile:
interesting = False
for cigarType, cigarLength in cigarLine:
if cigarType == 1:
if cigarLength > 100:  ## Arbitrary value
interesting = True
break
if interesting:

\$ java -jar dist/bioalcidaejdk.jar -e 'stream().filter(R->!R.getReadUnmappedFlag() && R.getCigar().getCigarElements().stream().anyMatch(C->C.getLength()>100 && (C.getOperator().equals(CigarOperator.N) || C.getOperator().equals(CigarOperator.D)))).forEach(R->println(R.getContig()+"\t"+R.getStart()+"\t"+R.getEnd()));' input.bam