Question: Combine multiple samtools view
0
gravatar for rah
19 days ago by
rah20
rah20 wrote:

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

bam samtools alignment • 84 views
ADD COMMENTlink modified 19 days ago by GenoMax95k • written 19 days ago by rah20

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 REPLYlink written 19 days ago by Pierre Lindenbaum133k

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 REPLYlink written 19 days ago by rah20
1
gravatar for Pierre Lindenbaum
19 days ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum133k wrote:

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 COMMENTlink written 19 days ago by Pierre Lindenbaum133k
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: 1682 users visited in the last hour
_