This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Filtering sam/bam by using CIGAR deletion sites

Hi,

I have a CLIP-Seq data and after mapping, I want to filter my sam files for reads containing deletion by using CIGAR & MAPQ >10. Is there a way in samtools to subset my sam file containing reads MAPQ>10 & with deletion?

Thank you very much!

alignment next-gen
• 6,109 views • link
updated by i.sudbery i.sudbery • written by ag1194 • ag1194

Thank you very much!!! I'll try and see samjdk first and also pysam to see the results.

• 1 views • link
written by ag1194 • ag1194

3 answers

I would do something like:

samtools view file.bam | awk '($6 ~ /D/) && ($5>10)'

($6 ~ /D/) means "deletion in CIGAR" and $5>10 means MAPQ>10

I really like awk oneliners for doing quick filters like this.

• 1 views • link
written by Geparada Geparada

With samtools view & awk tools, I cannot keep the header in my output bam. Is there a way to have filtered bam with header with awk tool? Thanks!!!

• 1 views • link
written by ag1194 • ag1194 
samtools view -h file.bam | awk '($0 ~ /^@/) || (($6 ~ /D/) && ($5>10))'
• 1 views • link
written by i.sudbery i.sudbery

Hey when I try to pipe this to a new bam like so

samtools view -h file.bam | awk '($0 ~ /^@/) || (($6 ~ /D/) && ($5>10))' > new.bam

I can't view the new.bam. Any ideas?

• 0 views • link
written by gt • gt

The result will be sam, not bam. To save it as sam, you need to pipe it to samtools view -b

• 0 views • link
written by i.sudbery i.sudbery

using samjdk: http://lindenb.github.io/jvarkit/SamJdk.html

java -jar dist/samjdk.jar -e 'return !record.getReadUnmappedFlag() && record.getMappingQuality()>10 && record.getCigar().getCigarElements().stream().map(C->C.getOperator()).anyMatch(O->O.equals(CigarOperator.D) || O.equals(CigarOperator.N)); ' input.bam
• 1 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Pysam is a good tool to manipulate SAM/BAM with various attributes like CIGAR and mapping quality.

https://pysam.readthedocs.io/en/latest/

• 1 views • link
written by Vitis Vitis 
import pysam
inbam = pysam.AlignmentFile("my_bamfile.bam")
outbam = pysam.AlignmentFile("filtered_bamfile.bam", template=inbam)
for read in inbam.fetch(until_eof=True):
    if read.mapping_quality > 30 and 'D' in read.cigarstring:
        outbam.write(read)
• 1 views • link
written by i.sudbery i.sudbery

Log in to answer this question.

More posts like this