Question: Filter Paired-End Sam File For Xt:A:U
0
gravatar for Steffi
7.4 years ago by
Steffi560
Germany
Steffi560 wrote:

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

filter paired-end sam bwa • 3.5k views
ADD COMMENTlink modified 5.7 years ago by Biostar ♦♦ 20 • written 7.4 years ago by Steffi560
4
gravatar for Leszek
7.4 years ago by
Leszek4.0k
IIMCB, Poland
Leszek4.0k wrote:

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 = ''
ADD COMMENTlink modified 7.4 years ago • written 7.4 years ago by Leszek4.0k

Thanks a lot! I appreciate your help!

ADD REPLYlink written 7.4 years ago by Steffi560

the way to thank on Biostar is to upvote and accept the answer ;-)

ADD REPLYlink written 7.4 years ago by Istvan Albert ♦♦ 81k
2
gravatar for Pierre Lindenbaum
7.4 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum122k wrote:

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)"
ADD COMMENTlink written 7.4 years ago by Pierre Lindenbaum122k

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.

ADD REPLYlink written 7.4 years ago by Steffi560

I see, I misunderstood your question.

ADD REPLYlink written 7.4 years ago by Pierre Lindenbaum122k

just do additional pair-end filtering after grep: samtools view -h -f 0x002 file.bam | egrep "(^@|XT:A:U)" | samtools view -Sh -f 0x002 -

ADD REPLYlink written 7.4 years ago by Leszek4.0k

This does not work, because a read can still have the flag of being properly mapped but having lost its mate at the grep-step before...

ADD REPLYlink written 7.4 years ago by Steffi560

Right, I didn't notice that. Anyway, below you have a few python lines I use for this:)

ADD REPLYlink written 7.4 years ago by Leszek4.0k
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: 1118 users visited in the last hour