Question: filter BAM file for reads that map 100% identity
1
gravatar for oriolebaltimore
2.6 years ago by
United States
oriolebaltimore150 wrote:

Hello, I have a BAM file and I want to create another BAM file by filtering only reads that are 100% identical mapped. For example if my read length is 100 , I want to select CIGAR 100M and make a bam file. Could you please suggest how this can be done. Thanks Adrian

samtools bam cigar • 1.7k views
ADD COMMENTlink modified 2.6 years ago by h.mon30k • written 2.6 years ago by oriolebaltimore150
2
gravatar for h.mon
2.6 years ago by
h.mon30k
Brazil
h.mon30k wrote:

To check for a given length of total match - e.g. 100 - you can use perl and anchor the regular expression: /^100M$/. This will exclude any other total match length (e.g. 95M or 125M) and any soft clipped read. You will have to re-header the file after filtering.

samtools view potexvirus2.bam \
    | perl -lane 'print if $F[5] =~ /^100M$/;'
ADD COMMENTlink modified 2.6 years ago • written 2.6 years ago by h.mon30k
1
gravatar for Pierre Lindenbaum
2.6 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum129k wrote:

using samjdk http://lindenb.github.io/jvarkit/SamJdk.html . Check the edit-distance (NM) exists and is equals to zero.

java -jar  dist/samjdk.jar -e 'return !record.getReadUnmappedFlag() && record.getCigarString().equals("100M") &&  (record.getIntegerAttribute("NM")==null || record.getIntegerAttribute("NM").intValue()==0);' in.bam
ADD COMMENTlink modified 2.6 years ago • written 2.6 years ago by Pierre Lindenbaum129k

Depending on your aligner, NM may or may not include clipped bases so you'll probably want to check the CIGAR as well.

ADD REPLYlink written 2.6 years ago by d-cameron2.2k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1010 users visited in the last hour