Exclude certain tags from BAM file

[yksikaksi]

Wishing you guys have a nice Mayday.

I want to exclude certain reference tags(genes) from BAM files and direct the output to a new BAM file. I tired with SAMtools commend as below:

foreach line (`cat tags2remain.txt `)
  samtools view 0H.ptMapped.bam | grep "$line" >>! 0H.tags2RemainMapped.bam
end

I obtained the output as

[bam_header_read] EOF marker is absent. The input is probably truncated.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
[main_samview] fail to read the header from "0H.tags2RemainMapped.bam"

Am I doing that correctly, kindly please give me some suggestion. OR, should I use Picard tools to do so? If yes, may I know what Picard tools I should use.

Thank you very much for the help.

[maubp]

Tags is the wrong word (that would be things like MD, or RG - information about each read). You are asking about filtering SAM/BAM files by read name.

It appears the for loop will try to concatenate BAM files together - that doesn't work.

How many read names are you dealing with? A single grep using a regular expression to describe all the read names at once could work (if possible given limits on command line lengths).

I'm not sure of any existing tool to tackle this particular need.

Personally I'd write a simple Python script to load all the names in memory (as a Python set for efficiency), then loop over the reads and filter them. Something similar would be easy in Perl or Ruby if you prefer, or indeed in Java using Picard.

[dawe]

bamtools can filter on sequence tags:

$ bamtools filter -h
Description: filters BAM file(s).

Usage: bamtools filter [-in <filename> -in <filename> ...] [-out <filename> | [-forceCompression]] [-region <REGION>] [ [-script <filename] | [filterOptions] ]

...
  -tag <TAG:VALUE>                  keep reads with this
...

HTH

d

[yksikaksi]

Thanks guys for your input.

It is a bit strange I don't receive feedback alert email from the forum.

Maybe I should rephrase my question. Put it in this way, I want to filter the genes (tags) which is let's say ribosomal protein from mapped reads .bam file.

I'm thinking @SQ could be the way to do, however, I'm not sure this is the right to do.

Thanks for your time and help.

[swbarnes2]

Why wouldn't this work:

samtools view unfiltered.bam | awk or grep or whatever | samtools view -bS > filtered.bam

You will likely lose your headers this way, unless you are careful with the awk or grep or whatever.

[yksikaksi]

Thanks swbarnes2 for your suggestion.