filter reads based on mapped length
Hi everyone, I’m trying to filter reads based on their mapped length, for example from 1000 bp onwards. Could you suggest something? Would these commands work?
samtools view -h /home/bar/input.bam | awk 'BEGIN {OFS="\t"} !/^@/ {len=0; n=split($6, CIGAR, /[MIDNSHP]/); for(i=1; i<n; i++) {if(index($6, CIGAR[i] "M")) {len+=CIGAR[i]}}} {if(len>=1000 || /^@/) print $0}' | samtools view -b > /home/bar/output.bam
Thank you very much for your support.
• 2,065 views
•
link
updated
by
John Marshall
•
written
by
marco.barr •
•
written
by
marco.barr •
1 answer
Use the filtering query language. samtools view -e 'rlen >= 1000' input.bam should do it.
Although that's the number of reference bases aligned against. It's not quite the same as the number of query bases, but maybe close enough. The difference is basically counting "D" cigars vs "I" cigars.
• 0 views
•
link
updated
by
John Marshall
•
written
by
jkbonfield
•
written
by
jkbonfield
Log in to answer this question.
More posts like this
-
Max length reads from fastq
written by marco.barr •
-
Split large deletion D in CIGAR string and filtering reads
written by marco.barr •
-
Filter out Bam not overlapping with Bed File (keeping a read and its mate)?
written by Eliveri •
-
samtools view extract from .bed list to .bam file
written by matthews.alix •
-
not primary alignment flag filtering does not work
written by buffealo •
-
Extract reads overlapping a specific region in bam file
written by User000
-
samtools - Finding Average Fragment size for pair?
written by oars
-
Extract reads where both pairs maps properly and none of the pairs are split from a bam file
written by rimjhim.roy.ch
-
Problem Converting BAM to SAM on to Amazon-S3 Using Samtools
written by gundalav
-
How To Combine Two Samtools Commands.
written by dustar1986
Can you, in plain english, say what the filtering should do, and on which fields of the SAM/BAM you want to apply it.
I think they're trying to apply it to the CIGAR string (counting the M's, representing the alignment match).
In theory, that should work but they should verify the command on a few BAM records (I don't have time to debug the awk code).
I guess they're dealing with long reads because they're talking about >1000bp alignments. Also, not sure why one would want to filter by mapped length -- that's typically the job of the alignment tool.