Question: How to select aligned reads below certain mapping Quality (from BWA)?
3
gravatar for kirannbishwa01
2.6 years ago by
United States
kirannbishwa01980 wrote:

Samtools can be used to select reads above certain mapping quality.

samtools view -h -b -q 30 aligned.bam -o above.mapQ30.bam

But, how to select a read below certain mapping quality - all aligned reads below mapQ 30?

I know it can be done using awk. But, the pipeline gets lengthy and time consuming when first need to convert bam to sam - separate header - use awk for mapQ below 30 - add header - sam file - convert to bam. Really, its taking a lots of time.

Thanks,

ADD COMMENTlink modified 2.6 years ago • written 2.6 years ago by kirannbishwa01980
1

It may be possible to do this with one of the many programs from BBMap (perhaps reformat.sh). Brian Bushnell may be along later to provide that solution.

ADD REPLYlink written 2.6 years ago by genomax65k
7
gravatar for harold.smith.tarheel
2.6 years ago by
United States
harold.smith.tarheel4.3k wrote:

Use SAMtools with the -U flag to output reads that don't match your filter, e.g.

samtools view -h -b -q 30 -U below_q30.bam aligned.bam

ADD COMMENTlink written 2.6 years ago by harold.smith.tarheel4.3k

Thanks Harold. I just realized that this will only work when using samtools version 1.3. The version samtools 1.19 which is more popular and installed by default in many linux station doesn't have this flag available.

So, folks make sure you have proper version of samtools installed. i.e 1.3 http://www.htslib.org/doc/samtools.html . And, thanks much for the script.

ADD REPLYlink written 2.6 years ago by kirannbishwa01980

There is slight modification required in this script.

ADD REPLYlink written 2.6 years ago by kirannbishwa01980
2
gravatar for kirannbishwa01
2.6 years ago by
United States
kirannbishwa01980 wrote:

I haven't found any straight forward method, but if someone has similar problem in the future; below is my proposed solution:

Note: Check and make sure the mapped alignment is sorted.

samtools view -h mapped.bam -o mapped.sam

grep '^@' mapped.sam > RG.header # this is important when there is header info in the sam/bam file

grep -v '^@' mapped.sam > SAM_file.sam # retains the non header part of the mapped file

awk '($5 < 40)' SAM_file.sam > Below.mapQ40.sam

cat RG.header Below.mapQ40.sam > Below.mapQ40_mapped.sam

samtools view -h -S -b Below.mapQ40_mapped.sam -o Below.mapQ40_mapped.bam

samtools index Below.mapQ40_mapped.bam

If you don't care about header info (or if its absent):

samtools view -h mapped.bam -o mapped.sam

awk '($5 < 40)' mapped.sam > Below.mapQ40_mapped.sam

samtools view -h -S -b Below.mapQ40_mapped.sam -o Below.mapQ40_mapped.bam

samtools index Below.mapQ40_mapped.bam

ADD COMMENTlink written 2.6 years ago by kirannbishwa01980
2
gravatar for kirannbishwa01
2.6 years ago by
United States
kirannbishwa01980 wrote:

Method B - requires samtools version 1.3

first select the mapped reads only samtools view -b -h -F 4 aligned.bam > mapped_only.bam # this selects for mapped reads only.

samtools index mapped_only.bam # index the bam file

Now separate the mapped reads above mapQ 40 and below samtools -b -h -U below_mapQ40.bam -q 40 mapped_only.bam -o above_mapQ40.bam # make sure to use right samtools version (1.3)

samtools index above_mapQ40.bam # this is bam file with reads with mapQ 40 and above

samtools index below_mapQ40.bam # this is bam file with reads with mapQ below 40.

ADD COMMENTlink modified 2.6 years ago • written 2.6 years ago by kirannbishwa01980
0
gravatar for harold.smith.tarheel
2.6 years ago by
United States
harold.smith.tarheel4.3k wrote:

I don't know if it's faster, but you could subtract the ≥Q30 file from the total using BEDtools intersect with -v flag.

Edit: The OP is correct; this was a bone-headed suggestion. See other answer.

ADD COMMENTlink modified 2.6 years ago • written 2.6 years ago by harold.smith.tarheel4.3k

Bedtools rely on interval (co-ordinates; region of the genome). So, in the same interval (region of the genome) there can be reads with both high mapQ and low mapQ. So, this won't end up working. I haven't tried it but I think this is what will happen.

still thinking !!!

ADD REPLYlink written 2.6 years ago by kirannbishwa01980
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: 760 users visited in the last hour