Thanks a lot! I appreciate your help!
Dear all,
I have a sam file (BWA output, paired-end reads). I would like to retain only reads which are "properly paired". This I would do by:
samtools view -f 0x002 file.sam > file_filtered.sam
Additionally I would like to retain only those pairs of reads where both reads have the XT:A:U tag. It is important to me that after the filterting step I still have the pairs together (so read1, read2, read1, read2, ...).
Any ideas how to do so?
Thanks for any help! Stefanie
2 answers
Simple scripting will do the job:
#!/usr/bin/env python
# Report uniquely aligned pairs from BWA sam output.
import sys
i = 0
uPair = 0
lines = ''
for l in sys.stdin:
#write header info
if l.startswith('@'):
sys.stdout.write( l )
continue
#reads
lines += l
i+=1
if 'XT:A:U' in l:
uPair += 1
#every two reads
if not i % 2:
#check if both are uniquely mapped
if uPair == 2:
sys.stdout.write( lines )
uPair = 0
lines = ''
the way to thank on Biostar is to upvote and accept the answer ;-)
I would just use 'grep' to keep the headers and the line matching the tag:
/samtools view -h -f 0x002 file.bam |\
egrep "(^@|XT\:A\:U)"
But this would (possibly) destroy the pairs: It might occur that one read of a pair is mapped with XT:A:U but the other not. In that case only one read of the pair would be retained.
Log in to answer this question.