<p>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:</p>

<pre class="lang-py prettyprint-override"><code>for read in bamfile:
    if not read.is_unmapped:
        cigarLine = read.cigar
        interesting = False
        for cigarType, cigarLength in cigarLine:
            if cigarType == 1:
                if cigarLength &gt; 100:  ## Arbitrary value
                    interesting = True
                    break
        if interesting:
            print("{}:{}-{}".format(read.reference_name, read.reference_start, read.reference_end))
</code></pre>

<p>Unless you have long reads, this will be close enough for IGV.</p>

CC BY-SA 3.0
2018-02-08T20:18:12.540Z
Well done! I obviously need to study the pysam API a bit more closely. Thanks for the help!
ShanZhengYang
CC BY-SA 3.0
2018-02-08T20:06:06.460Z
2018-02-08T20:06:06.460Z
Devon Ryan
<p>a one-liner using <strong>bioalcidaejdk</strong>: <a href="http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html" rel="nofollow noreferrer">http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html</a></p>

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

CC BY-SA 3.0
2018-02-09T14:13:31.867Z
2018-02-09T14:13:31.867Z
Pierre
<p>I would like to find out where (in genomic coordinates) large insertions are found within a given BAM file, <code>file.bam</code>.</p>

<p>(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.)</p>

<p>Based on the BAM CIGAR, I should be able to <code>grep</code> the BAM based on some cut-off of "large inserted" bases. </p>

<p>It's unclear to me how to best do this via samtools; one could parse the cigar string via <code>samtools</code> for insertions, but I don't understand how to translate this into the genomic position in the alignment. </p>

<p>This goes with various samtools wrappers, e.g. <code>pysam</code> in python:</p>

<pre><code>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
</code></pre>

<p>I'm not sure how to do this in a non-clumsy manner. Any help appreciated.</p>

<sam><samtools>
Parse BAM by insertion size and get genomic coordinates
# Parse BAM by insertion size and get genomic coordinates

(Asked by: username_0 on 2018-02-08)

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

**Tags:** sam, samtools

---

## 2 answers

### Answer by username_1 on 2018-02-08. Score: 2 | Accepted answer

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:
    if not read.is_unmapped:
        cigarLine = read.cigar
        interesting = False
        for cigarType, cigarLength in cigarLine:
            if cigarType == 1:
                if cigarLength > 100:  ## Arbitrary value
                    interesting = True
                    break
        if interesting:
            print("{}:{}-{}".format(read.reference_name, read.reference_start, read.reference_end))

Unless you have long reads, this will be close enough for IGV.

> **Comments:**
> **username_0**: Well done! I obviously need to study the pysam API a bit more closely. Thanks for the help!  

---

### Answer by username_2 on 2018-02-09. Score: 1

a one-liner using bioalcidaejdk: http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html
$ 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


---

