3
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Well done! I obviously need to study the pysam API a bit more closely. Thanks for the help! $\endgroup$ – ShanZhengYang Feb 8 '18 at 20:18
1
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.