Adding XS flag to Novoalign reads

[zauhar]

Hi folks, here is the scenario -

I have unpaired, unstranded reads that have already been aligned using Novoalign by another lab, and I have the resulting BAM files. Our pipeline uses GSANP/cufflinks, and I want to compare DE results for the two pipelines. (I also have the reads from the other lab, and aligned using GSNAP.)

Cufflinks wants to see the XS flag for spliced reads, which are not generated by Novoalign, and throws an error for such reads, apparently ignoring them.

Nonetheless, the FPKMs found for most genes actually compare very closely for the two pipelines, which is encouraging - but there are also some big discrepancies, so I want to add the XS flag to the Novoalign reads and make sure they are being counted. (I should add that for this test I am comparing the same sample using output from the two different aligners, which used the same reference genome; so in principle all FPKMs should be close.)

I tried bamutils (-xs option) for this - but it adds the XS tag to ALL Novoalign reads, and also tinkers with the some of the other tags - now the FPKMs are wildly different from GSNAP! A quick test where I kept just the unspliced reads from both pipelines shows that either cufflinks does not know how to handle the XS flag when applied to all reads, or the other changes made by bamutils have screwed things up.

Anyone have experience with this, or suggestions for another tool to add the XS flag? I can probably figure it out myself, but I already invested a lot more time in this than expected (of course that is not unexpected, is it?)

Thanks!

Randy

[awenocur]

An idea that occurs to me is to use a regex to detect splice junctions, and append the tag to all lines that match it.

[zauhar]

Yes, that is true - I can just look for the spliced aligns, and use the SAM flag to generate the XS tag, trusting that Novoalign made use of the splice signal info.

I am curious though why the addition by bamutiils so completely bollocksed things up. The other changes to the tags look innocuous (data format changes, 'A' instead of 'Z' for a single character chromosome name).

Thanks Adam!