Duplicate `MC` tags in samblaster output

[DonFreed]

Starting in version 0.7.16, BWA produces SAM files containing MC tags. Running samblaster with --addMateTags on BWA's output then produces a output file with duplicate MC tags.

We would like to use samblaster to add MQ tags to BWA's output, but the duplicate MC tags are causing issues for downstream tools.

[GregoryFaust]

I have known for a while about the possibility of duplicate MC tags using samblaster with recent versions of bwa. However, I have never known of any downstream tools having a problem with this input. For curiosity, what tools in your pipeline are having issues? In the mean time, I will check to see how much it will slow down samblaster to check for and correct this problem.

[DonFreed]

Sorry I missed the notification here, and I hope I didn't overstate the problem. The duplicate MC tags are causing two minor problems:

1. Occasionally, bwa and samblaster will produce different MC tags. This seems to mostly occur on secondary alignments.
2. This increases the output file size.

[GregoryFaust]

Thanks for your response. I have noticed these differing MC tags before, but have not found a definitive answer for the issue. samblaster always appends the CIGAR string from the primary alignment of the mate. A recent test run that I performed on the 50x platinum sequence data from Illumina for sample NA12878 showed no disagreement between the BWA MC tag values and those from samblaster. What version of BWA were you using? Maybe the issue has been fixed in recent BWA versions? Or did you notice where the CIGAR string used by BWA came from? I will look harder into this to see if I can again find some disagreeing data, especially around the issue of the use of the -M flag on either BWA and/or samblaster.

In the mean time, I have prepared a version of samblaster for release in the near future that will no longer produce duplicate mate tags pending a determination of what to do about the above, and also issue #20 for which there continues to be growing interest.

[DonFreed]

Excellent, I'm looking forward to trying the new version.

I'm able to produce different MC tags with BWA 0.7.17 and the the Illumina Platinum Genomes run ERR174327. Here's an example command using a short script to detect the mismatched tags, https://gist.github.com/DonFreed/9d59e758ff5780261fdeb8a03e714717:

$ bwa-0.7.17 mem -R "@RG\tID:TEST\tSM:TEST" -t 32 -K 10000000 hs37d5.fa \
    ERR174327_1.fastq.gz ERR174327_2.fastq.gz | samblaster-0.1.23/samblaster --addMateTags \
    --acceptDupMarks -i stdin -o stdout | python3 different_tags.py --check
Different MC tags: ERR174327.607        2209     5       120089979       0       41H33M27H       hs37d5  26062971        0       ATATATAATATATAATATATAATATATAACATA       JIIJJJJJJIJJJJJJJJJJJJJJJJJJJJJJI    NM:i:0  MD:Z:33 MC:Z:5M3D90M6H  AS:i:33 XS:i:32 RG:Z:TEST       SA:Z:hs37d5,28674526,-,46S42M13S,25,0;  MC:Z:5M3D90M6S  MQ:i:6

In this test case, all of the mismatched tags come from the supplementary alignments of split reads. BWA produces MC tags with hard-clips, while Samblaster has soft-clipped tags.

Here are the BWA alignments for the read-pair from the above case:

ERR174327.607   113     hs37d5  26062971        6       5M3D90M6S       =       28674526        2611500 TATAATAAATATTTATATATACTATAATATAAAACATATATAATATATATCATATAGTATATATTATATACCGTATATATTATATATCATATATTGTATAT   FFHHJIJIIJJIJIIGIIIIIIIIIHIIJJJIIFJJIIJIJJJIIGJJJJJIIJJIIJJJJJJJJJJIHGIFIJIIHIIJIJJJJJJJHHHGHFFFFFCCC   NM:i:11 MD:Z:5^TTT6G9A32A1A5T2G6A10G11  MC:Z:46S42M13S  AS:i:50 XS:i:45 RG:Z:TEST       XA:Z:hs37d5,-26065287,5M3D90M6S,12;hs37d5,+28674355,31S19M1I8M2I40M,5;hs37d5,+25148381,14M2I63M22S,8;   MC:Z:46S42M13S  MQ:i:25
ERR174327.607   177     hs37d5  28674526        25      46S42M13S       =       26062971        -2611500        TATATTTCATGTTATATCAACATATATTATGTTATATATTATATATTATATATTATATATCACATATTACATTATATATTATATATCAAATATATGTTACA   HHHJJJJJIJJJJJJJJJJIJJJJJJIIJJJJJJJJJJJJJJJJJJJJJJIJJJJJJIIJIIIJJJJJJJJJIIIJJJJJJJJFIIJJHHGHHFFFFFCCC   NM:i:0  MD:Z:42 MC:Z:5M3D90M6S  AS:i:42 XS:i:37 RG:Z:TEST       SA:Z:5,120089979,+,41S33M27S,0,0;       XA:Z:16,+16873121,15S15M2I40M29S,4;4,-169077189,34S35M3I14M15S,4;hs37d5,+670448,13S14M5I36M33S,6;       MC:Z:5M3D90M6S  MQ:i:6
ERR174327.607   2209    5       120089979       0       41H33M27H       hs37d5  26062971        0   ATATATAATATATAATATATAATATATAACATA       JIIJJJJJJIJJJJJJJJJJJJJJJJJJJJJJI       NM:i:0  MD:Z:33 MC:Z:5M3D90M6H  AS:i:33 XS:i:32 RG:Z:TEST       SA:Z:hs37d5,28674526,-,46S42M13S,25,0;  MC:Z:5M3D90M6S  MQ:i:6

Samblaster is correctly reporting the mate alignment as soft-clipped.

[GregoryFaust]

I looked into this more while working on release 0.1.25. As you suggested, I only found mismatched tags on supplementary alignments, but seemingly only when the primary alignment of the mate read had a large clip. In such cases, again as you suggest, bwa mem showed a hard-clipped mate CIGAR instead of the correct soft-clipped CIGAR. However, all downstream processing tools of which I am aware, with the exception of bam2fastq, will essentially treat the soft and the hard clipped CIGAR the same. They just want to know how the alignment matches to the reference.

In any event, I have fixed the issue of duplicate MC (or MQ) tags in release 0.1.25. samblaster does not correct any existing MC or MQ tags if it has a value that is different than what is already in the input file. In my view, samblaster is not in the business of making changes to the output from the aligner in such cases.

[DonFreed]

Thanks for the fix!

In my view, samblaster is not in the business of making changes to the output from the aligner in such cases.

Sounds quite reasonable to me.