Combine multiple samtools view
1
0
Entering edit mode
3.9 years ago
rah ▴ 30

Can anyone help me with combining two samtools view commands into a single command.

I need to remove some alignment tags "XA and SA" and afterwards i also need to keep only the aligned reads to the "canonical" chromosomes from a reference genome.

It is possible to do all of this in multiple steps as so (with index created in between):

1) remove tag - taken from other biostar posts.

samtools view file.bam | grep 'XA\|SA' | cut -f 1 > discard_ID.txt

samtools view -h file.bam | fgrep -vf discard_ID.txt | samtools view -Sb - > out.bam

2) extract aligned reads

samtools view -F 4 -q 30 out.bam chr{1..22} chrX chrY chrM -b > final.bam

But i would like to combine these into a single command. So i've tried the following

samtools view -h file.bam | fgrep -vf bad_names.txt | samtools view -F 4 -q 30 - chr{1..22} chrX chrY chrM -b > test.bam

the problem is the different intermediate bam files are of course not index, so is there a way to combine 1) and 2) into a single command. I've also tried adding the -h option to samtools view to see if the output with header is able to be passed on in the pipe.

Thanks a lot

samtools bam alignment • 1.3k views
ADD COMMENT
0
Entering edit mode

samtools view file.bam | grep 'XA\|SA' | cut -f 1 > discard_ID.txt

that's wrong. you would delete any read with a name containing XA.

ADD REPLY
0
Entering edit mode

yeah sorry for the mistake, thats what i meant by filtering out. I want to remove all the reads with XA and SA tag. Afterwards i want to extract all the reads which aligns to chr 1-22, x,y,m

ADD REPLY
1
Entering edit mode
3.9 years ago

using samjdk : http://lindenb.github.io/jvarkit/SamJdk.html

samtools view -F 4 -q 30 - chr{1..22} chrX chrY chrM -u | \
java -jar dist/samjdk.jar -e 'Object att = record.getAttribute("XA"); if(att==null)  att = record.getAttribute("SA");  return att==null ;'
ADD COMMENT

Login before adding your answer.

Traffic: 1555 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6