This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Remove both pair end reads with low mapping quality using Samtools

To my understanding, when use samtools -q 20 on paired end sequencing data, if one read of a pair has quality 30, the other has 0, only the 0 quality one will be removed. So I'm wondering if samtools has a option to remove both, as I'm trying to test my program and I prefer my reads stay paired.

Thank you very much!

mapping-quality samtools
• 4,786 views • link
updated by Ram Ram • written by DVA DVA

Does that actually happen? It must be fairly rare and, what samtools utility are you running? If it's mpileup, then you don't need to worry about reads as it's doing the per-base pileup. In short, you're either creating unnecessary work for yourself or it'd be good to give more information about what you are trying to do.

• 0 views • link
updated by Ram Ram • written by brentp brentp

brentp : I've just tested a few bams from 1000G, it doesn't happen often, but it happens:

$ curl -s "ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/data/150140/alignment/150140.chrom20.ILLUMINA.bwa.CHM1.20131218.bam"  |\
samtools view  - |\
head -n 20000 | \
grep -E '(SRR642643.96013593|SRR642643.4610992|SRR642643.101471738)' | sort -k1,1 | cut -f 1-9

SRR642643.101471738    163    20    61043    37    100M    =    61133    124
SRR642643.101471738    83    20    61133    36    65S35M    =    61043    -124
SRR642643.4610992    163    20    61252    29    38M62S    =    61340    187
SRR642643.4610992    83    20    61340    37    100M    =    61252    -187
SRR642643.96013593    163    20    61234    29    40M60S    =    61322    187
SRR642643.96013593    83    20    61322    37    100M    =    61234    -187
• 0 views • link
updated by Ram Ram • written by Pierre Lindenbaum Pierre Lindenbaum

you mean because one end has a qual of 37 and the other of 29?

• 0 views • link
written by brentp brentp

It happens to me when I try to test some of my own software. I guess it is very rare in real sequencing analysis. Thank you.

• 0 views • link
written by DVA DVA

1 answer

Samtools doesn't have such an option (it would be too rarely needed to incorporate it). You can query sort the BAM file and write a little Python script with pysam to process the pairs as a unit. That'd be a simple enough way to produce the same result.

• 0 views • link
written by Devon Ryan Devon Ryan

Thanks! I believe there is no simple ways to get around with it:)

• 0 views • link
written by DVA DVA

Log in to answer this question.

More posts like this