jmarshall
commented
1 year ago The exception encountered (Empty array passed to setUnsignedArrayAttribute for tag [tagname]) is not specific to ML but would be triggered for any tagged field. Zero-length arrays are explicitly valid in general since 2018 (see #135 and #326), and HTSJDK (and hence Picard and GATK) was adjusted accordingly in samtools/htsjdk#1194. Unfortunately setUnsignedArrayAttribute() was overlooked in that pull request, but it's clear that its check should have been removed too and this exception being thrown is an HTSJDK bug.
So much for the proximate cause. However there may be HTSJDK base modification code that would have been run for ML in particular but did not execute as this exception was thrown first, that might not handle a zero-length array; i.e., there may be additional reasons why HTSJDK might not handle a zero-length ML field — so we still need an answer to whether an empty ML array is intended to be valid.
In my opinion, such an ML value is very much valid and HTSJDK ought to accept it without issue.
SAMtags.pdf §1.7 also says
[ML] uses a byte array of type ‘C’ with the number of elements matching the summation of the number of modifications listed as being present in the MM tag
This is quite explicit that for a record that contains MM:Z:C+m; (which ¶8 of the MM description explicitly says is valid), that MM value lists a total of zero modifications, so the corresponding ML array has zero elements.
Practically speaking, it would be silly to require methylation calling software to include code to check whether their MM tag being emitted lists no modifications and suppress emitting an ML field accordingly. Far better to regularise the corner case of a correspondingly empty ML value.
So the correct interpretation is your (1), and this is simply an HTSJDK bug.
One commenter on the GATK discussion considered that it is unclear whether the SAM spec allows an empty ML value or not. We could add explicit text noting that an empty ML:B:C array would correspond to the largely empty MM:Z:C+m; value that is explicitly discussed. However personally I don't think that's necessary, as it is fairly clear (from the equality in the paragraph I quoted above) — easily derivable albeit not explicit — and so the only reason the GATK commenter felt that there was any doubt was that they were misled by the faulty HTSJDK / GATK behaviour.
jkbonfield
Currently, htsjdk will throw an exception if the ML record is present but empty. This post on the GATK forum documents an example of this exception in the wild. However, after reading section 1.7 on Modified Bases in the SAMtags spec, it is unclear to me if htsjdk is correctly interpreting this constraint. For example, is this arrangement of records out of spec?
< ... snip ...> MM:Z:C+m?; ML:B:C < ... snip ...>To quote section 1.7:
and
If I understand this correctly, the expected count of ML entries is driven by the number of list of modifications mentioned by MM. However, since
MM:Z:C+m;explicitly indicates no modifications occurred, there are no events for ML to annotate.I think there are two possible but mutually exclusive interpretations to this use case, but the question is which one?
1) Empty ML arrays are legal, and htsjdk is acting too stringently 2) Empty ML arrays are illegal, and should only be populated in records when the array is not empty