Using Gatk, Removing Read Contigs On Bam Level

Question

[Mdeng]

Hey everyone,

ATM I am working with GATK for the first time and a get an error, wich everybody knows:

##### ERROR MESSAGE: Input files reads and reference have incompatible contigs: Order of contigs differences, which is unsafe.
##### ERROR   reads contigs = [chr1, chr2, ..., chrUn_gl000226, chr18_gl000207_random]
##### ERROR   reference contigs = [chr1, chr2, ..., chrY, chrM]

Is there a way to remove the chr18_gl000207_random on the BAM level, I don't wanna go back to SAM level and kick them out. Because if already done some postprocessing now.

With best,

Mario

[lh3]

What if you change the reference sequence instead of the alignment?

[Mdeng]

This would work to, but I am looking for nice and fast why, Istvan Albert posted below.

Answer 1

[Istvan Albert]

You could list all the valid chromosomes as a query like so:

samtools view  -b old.bam chr1 chr2 chr3 > newbam

(old.bam file needs to be sorted and indexed)

[lh3]

I guess this does not work because the sequence dictionary is not changed. I think GATK should adopt the faidx strategy.

[Mdeng]

@lh3 Couldn't this be solved by recreating the SAM header using Picards AddOrReplaceReadGroups?

[Russh]

Do you mean ReplaceSamHeader?

Answer 2

[Russh]

I don't quite follow why you don't want to go back to SAM level. Surely making a temporary SAM, grepping that, and converting to BAM is the easiest way to do this, isn't it?

The alternative, is to use picard's ReorderSam tool, which, given that the offending contigs are absent from your reference, will drop any reads that are mapped to those contigs