Handling alignments spanning chromosomes boundaries in BWA

[nirav99]

Hi All,

When BWA aligns a read that spans over two chromosomes, it sets the flag as "unmapped" but leaves mapping quality (MAPQ) and CIGAR unchanged.

Picard's validator tool ValidateSamFile throws a validation error for such reads because it expects an unmapped read to have MAPQ of zero, and CIGAR value of "*".

I would like to understand the generally acceptable way of handling this type of reads. Should I modify these reads to have MAPQ 0 and CIGAR value of * ?

I would appreciate any suggestions and best practices to follow for such reads.

Thanks in advance.

[swbarnes2]

Most people turn the stringency down on Picard to lenient. It'll yell at you about those reads, but it will go through. You'll want to do that with every Picard command.

[nirav99]

Hi swbarnes2,

Yes, we can have picard continue through these errors by the lenient stringency. But I was wondering what would be the correct way to handle the reads themselves so that no other downstream tool (not necessarily picard) will choke.

[drio]

I guess what he wants to know is if there is anyone out there that is doing more than ignoring those cases. Maybe patching BWA or rewriting the fields in the BAM. I can see the big sequencing centers doing something like before submitting their data.

[ECO]

It would be very interesting to understand what proportion of these type of reads are generated by chimerism introduced in the library prep steps, versus biologically significant events.

My guess is that the fraction caused by the library prep is fairly high, and people just haven't understood that yet because 1) the read lengths haven't been long enough to map both halves, and 2) most people are already overwhelmed by the reads they can map so throwing out a few % that map ambiguously is pretty common practice.

Cool stuff.

[nilshomer]

Spanning two chromosomes in BWA (short-read) comes from concatenating all chromosomes into one long sequence for use with the BWT and FM-index. So the case where a read is unmapped but has a chr/pos occurs when an alignment spans a chr/chr junction, where the two chromosomes are consecutive. This is different from the read spanning two random locations in the genome on different chromosomes (inter-chromosomal translocations etc.). The BWA (long-read) algorithm will actually produce the latter on longer reads (>100bp), giving two non-overlapping mappings for a read.

[drio]

Now that I think more about it, BWA is following the SAM spec. If the unmapped flag is set, it doesn't matter what you have in the other fields. Downstream tools should take that into account.

So I'd suggest you tell picard to don't report those warnings, or just ignore them. I guess that is what other people do.

[nilshomer]

The maintainers of the spec and the Picard developers have nonetheless decided to call it invalid. Better to patch BWA.

[nirav99]

I wrote a simple tool to handle this based on the API provided by Picard. Basically, it is an extension of Picard's CleanSam.jar with the change that if a read has unmapped flag, CIGAR is set to * and mapping quality to zero.

If anyone is interested in looking at it / using it, I can have it available in a public place.

[vschulz]

Hi nirav99,

Could you post or make your code available? I think that I am having this problem on the GATK GenomeAnalysisTK.jar -T IndelRealigner command, which doesn't allow the VALIDATION_STRINGENCY=LENIENT parameter.

Thanks,

Vince

[nirav99]

Hi Vince,

Download and build the java code at https://github.com/nirav99/Illumina-.../java/FixCIGAR

This should help you.

Nirav

[daphnawv]

Hi Nirav
I would also like to fix my reads to have MAPQ 0 and CIGAR value of * using your FixCIGAR tool. however Java is so overwhelming for me could you please explain in more detail how to build your script and show an example of a command line?

i also tried filtering only the mapped reads using ViewSam.jar ALIGNMENT_STATUS=Aligned but it didn't work, don't these type of reads have a flag of 4?

Thank you!
Daphna

[nirav99]

Hi Daphna,

I am assuming that you are running on Unix/ Linux and you have the latest java compiler. If so, you can follow these steps :

1) Download and untar / (unzip) latest Picard. We need two Jar files sam-*.jar and picard-*.jar.

2) Download the files from the FixCIGAR directory

3) Open GenerateJar.sh and edit the picardPath variable to the location where you have installed the picard library. You should point to the directory path where sam-*.jar and picard-*.jar are located.

4) Run GenerateJar.sh, it should build the tool for you.

5) It will create a .jar file FixCIGAR.jar. You can run it using the command

java -jar FixCIGAR.jar INPUT=bamToFix

or

java -jar FixCIGAR.jar INPUT=bamToFix OUTPUT=FixedBAM

[daphnawv]

Thank you so much for your help!