Skip to content

Commit 2ae4aaa

Browse files
tfennetfenne
authored
Implementation of a SAMRecord comparator that matches samtools' queryname sort order. (#1600)
* Implementation of a SAMRecord comparator that matches samtools' queryname sort order.
1 parent 1fda1a9 commit 2ae4aaa

2 files changed

Lines changed: 300 additions & 0 deletions

File tree

‎src/main/java/htsjdk/samtools/SAMRecordQueryNameNaturalComparator.java‎

Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
1+
/*
2+
* The MIT License
3+
*
4+
* Copyright (c) 2009 The Broad Institute
5+
*
6+
* Permission is hereby granted, free of charge, to any person obtaining a copy
7+
* of this software and associated documentation files (the "Software"), to deal
8+
* in the Software without restriction, including without limitation the rights
9+
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10+
* copies of the Software, and to permit persons to whom the Software is
11+
* furnished to do so, subject to the following conditions:
12+
*
13+
* The above copyright notice and this permission notice shall be included in
14+
* all copies or substantial portions of the Software.
15+
*
16+
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17+
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18+
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19+
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20+
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21+
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
22+
* THE SOFTWARE.
23+
*/
24+
package htsjdk.samtools;
25+
26+
import java.io.Serializable;
27+
28+
/**
29+
* Comparator for natural "queryname" ordering of SAMRecords, matching the order that samtools/htslib
30+
* produce for a queryname-natural sort (i.e. SAM {@code @HD SO:queryname SS:queryname:natural}).
31+
*
32+
* <p>Read names are compared so that maximal runs of digits are ordered by numeric value rather than
33+
* lexicographically, so e.g. {@code read2} sorts before {@code read10}. Everything else (the pairing,
34+
* strand, secondary/supplementary and {@code HI} tie-breaking) is inherited unchanged from
35+
* {@link SAMRecordQueryNameComparator}; only the read-name comparison differs.</p>
36+
*
37+
* <p>The read-name comparison is a faithful port of htslib's {@code strnum_cmp}. It walks the two
38+
* names digit-by-digit rather than parsing numeric runs into integers, so it produces the same result
39+
* as samtools even for numeric runs too large to fit in a {@code long}.</p>
40+
*/
41+
public class SAMRecordQueryNameNaturalComparator extends SAMRecordQueryNameComparator implements Serializable {
42+
private static final long serialVersionUID = 1L;
43+
44+
/**
45+
* Compares two read names using natural ordering, matching htslib's {@code strnum_cmp}.
46+
*
47+
* <p>Maximal runs of ASCII digits are compared as numbers; all other characters are compared by
48+
* their code point. Leading zeros are ignored, so numeric runs with the same value compare equal
49+
* (e.g. {@code 8}, {@code 08} and {@code 008}), matching samtools. Because the comparison never
50+
* parses a run into a number, it is correct for numeric runs of arbitrary length.</p>
51+
*
52+
* @return a negative number if {@code a < b}, zero if they are equal, a positive number if {@code a > b}
53+
*/
54+
public static int compareNatural(final String a, final String b) {
55+
final int lenA = a.length();
56+
final int lenB = b.length();
57+
int ia = 0;
58+
int ib = 0;
59+
60+
while (ia < lenA && ib < lenB) {
61+
final char chA = a.charAt(ia);
62+
final char chB = b.charAt(ib);
63+
64+
if (isDigit(chA) && isDigit(chB)) {
65+
// Skip leading zeros in each run so equal-value runs agree on their significant digits.
66+
while (ia < lenA && a.charAt(ia) == '0') ia++;
67+
while (ib < lenB && b.charAt(ib) == '0') ib++;
68+
69+
// Skip the digits that the two runs have in common.
70+
while (ia < lenA
71+
&& ib < lenB
72+
&& isDigit(a.charAt(ia))
73+
&& isDigit(b.charAt(ib))
74+
&& a.charAt(ia) == b.charAt(ib)) {
75+
ia++;
76+
ib++;
77+
}
78+
79+
final boolean aHasDigit = ia < lenA && isDigit(a.charAt(ia));
80+
final boolean bHasDigit = ib < lenB && isDigit(b.charAt(ib));
81+
82+
if (aHasDigit && bHasDigit) {
83+
// Both runs still have digits at the first differing position. The run with more
84+
// remaining digits is the larger number; if they are the same length, the differing
85+
// digit decides.
86+
int i = 0;
87+
while (ia + i < lenA && ib + i < lenB && isDigit(a.charAt(ia + i)) && isDigit(b.charAt(ib + i))) {
88+
i++;
89+
}
90+
if (ia + i < lenA && isDigit(a.charAt(ia + i))) return 1;
91+
if (ib + i < lenB && isDigit(b.charAt(ib + i))) return -1;
92+
return a.charAt(ia) - b.charAt(ib);
93+
} else if (aHasDigit) {
94+
return 1; // a's numeric run is longer, so it is the larger number
95+
} else if (bHasDigit) {
96+
return -1; // b's numeric run is longer, so it is the larger number
97+
}
98+
// Otherwise the two runs have the same numeric value. Leading zeros do not affect the
99+
// ordering (matching samtools), so we simply keep comparing the characters that follow.
100+
} else {
101+
if (chA != chB) return chA - chB;
102+
ia++;
103+
ib++;
104+
}
105+
}
106+
107+
// One name is a prefix of the other (or they are equal); the longer name sorts last.
108+
if (ia < lenA) return 1;
109+
if (ib < lenB) return -1;
110+
return 0;
111+
}
112+
113+
/** Returns true if the character is an ASCII digit. Read names are ASCII per the SAM spec. */
114+
private static boolean isDigit(final char c) {
115+
return c >= '0' && c <= '9';
116+
}
117+
118+
/**
119+
* Compares the read names of the two records using natural ordering. This is the only part of the
120+
* queryname comparison that differs from {@link SAMRecordQueryNameComparator}; the remaining
121+
* tie-breaking is inherited.
122+
*
123+
* @return negative if {@code samRecord1 < samRecord2}, 0 if their read names are equal, else positive
124+
*/
125+
@Override
126+
public int fileOrderCompare(final SAMRecord samRecord1, final SAMRecord samRecord2) {
127+
return compareNatural(samRecord1.getReadName(), samRecord2.getReadName());
128+
}
129+
}

‎src/test/java/htsjdk/samtools/SAMRecordQueryNameNaturalComparatorTest.java‎

Lines changed: 171 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,171 @@
1+
/*
2+
* The MIT License (MIT)
3+
*
4+
* Copyright (c) 2017 Daniel Gomez-Sanchez
5+
*
6+
* Permission is hereby granted, free of charge, to any person obtaining a copy
7+
* of this software and associated documentation files (the "Software"), to deal
8+
* in the Software without restriction, including without limitation the rights
9+
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10+
* copies of the Software, and to permit persons to whom the Software is
11+
* furnished to do so, subject to the following conditions:
12+
*
13+
* The above copyright notice and this permission notice shall be included in all
14+
* copies or substantial portions of the Software.
15+
*
16+
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17+
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18+
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19+
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20+
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21+
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
22+
* SOFTWARE.
23+
*/
24+
25+
package htsjdk.samtools;
26+
27+
import htsjdk.HtsjdkTest;
28+
import htsjdk.samtools.util.CollectionUtil;
29+
import java.util.ArrayList;
30+
import java.util.Collections;
31+
import java.util.List;
32+
import org.testng.Assert;
33+
import org.testng.annotations.Test;
34+
35+
public class SAMRecordQueryNameNaturalComparatorTest extends HtsjdkTest {
36+
private static final SAMRecordQueryNameNaturalComparator COMPARATOR = new SAMRecordQueryNameNaturalComparator();
37+
38+
/** Convenience wrapper returning the sign (-1/0/1) of a natural read-name comparison. */
39+
private static int cmp(final String a, final String b) {
40+
return Integer.signum(SAMRecordQueryNameNaturalComparator.compareNatural(a, b));
41+
}
42+
43+
@Test
44+
public void mixedNamesSortInNaturalOrder() {
45+
// A set of names with a single, unambiguous natural ordering (no leading-zero ties, which
46+
// samtools treats as equal). Numeric runs order by value: abc03 (3) < abc5 (5) < abc8 (8) < abc17 (17).
47+
final List<String> names = CollectionUtil.makeList(
48+
"abc", "abc+5", "abc- 5", "abc.d", "abc03", "abc5", "abc8", "abc17", "abc17.+", "abc17.2", "abc17.d",
49+
"abc59", "abcd");
50+
51+
final List<String> shuffled = new ArrayList<>(names);
52+
Collections.reverse(shuffled); // start from a non-sorted order so sorting does real work
53+
shuffled.sort(SAMRecordQueryNameNaturalComparator::compareNatural);
54+
55+
Assert.assertEquals(shuffled, names);
56+
}
57+
58+
@Test
59+
public void numericRunsCompareByValueNotLexicographically() {
60+
// The whole point of natural ordering: read2 < read10 even though "10" < "2" lexicographically.
61+
Assert.assertEquals(cmp("read2", "read10"), -1);
62+
Assert.assertEquals(cmp("read10", "read2"), 1);
63+
Assert.assertEquals(cmp("read9", "read10"), -1);
64+
}
65+
66+
@Test
67+
public void pureNumericNamesSortByValue() {
68+
Assert.assertEquals(cmp("2", "10"), -1);
69+
Assert.assertEquals(cmp("10", "100"), -1);
70+
Assert.assertEquals(cmp("100", "99"), 1);
71+
}
72+
73+
@Test
74+
public void leadingZerosDoNotAffectOrder() {
75+
// Numeric runs with equal value compare equal regardless of leading zeros, matching samtools.
76+
Assert.assertEquals(cmp("x008", "x08"), 0);
77+
Assert.assertEquals(cmp("x08", "x8"), 0);
78+
Assert.assertEquals(cmp("x008", "x8"), 0);
79+
// All-zero runs are likewise equal regardless of how many zeros.
80+
Assert.assertEquals(cmp("00", "0"), 0);
81+
}
82+
83+
@Test
84+
public void runsWithDifferentValuesAreOrderedRegardlessOfLeadingZeros() {
85+
// "010" (value 10) vs "09" (value 9): compared by value, not by the leading zero.
86+
Assert.assertEquals(cmp("010", "09"), 1);
87+
Assert.assertEquals(cmp("007", "8"), -1);
88+
}
89+
90+
@Test
91+
public void largeNumericRunsDoNotOverflow() {
92+
// Numeric runs beyond Long.MAX_VALUE (9223372036854775807) must still order by value.
93+
// A run parsed into a long would overflow here; the digit-by-digit comparison must not.
94+
final String twoPow64 = "read18446744073709551616"; // 2^64, 20 digits
95+
Assert.assertEquals(cmp("read9", twoPow64), -1, "single digit must sort before a 20-digit number");
96+
Assert.assertEquals(cmp(twoPow64, "read9"), 1);
97+
98+
// Two 20-digit numbers differing only in the final digit.
99+
Assert.assertEquals(cmp("read18446744073709551616", "read18446744073709551617"), -1);
100+
101+
// A 21-digit number is larger than any 20-digit number.
102+
Assert.assertEquals(cmp("read99999999999999999999", "read100000000000000000000"), -1);
103+
}
104+
105+
@Test
106+
public void equalNamesCompareEqual() {
107+
Assert.assertEquals(cmp("read1", "read1"), 0);
108+
Assert.assertEquals(cmp("abc17.2", "abc17.2"), 0);
109+
Assert.assertEquals(cmp("", ""), 0);
110+
}
111+
112+
@Test
113+
public void shorterPrefixSortsBeforeLongerName() {
114+
Assert.assertEquals(cmp("abc", "abc0"), -1);
115+
Assert.assertEquals(cmp("abc", "abcd"), -1);
116+
Assert.assertEquals(cmp("read", "read1"), -1);
117+
}
118+
119+
@Test
120+
public void emptyNameSortsBeforeAnyNonEmptyName() {
121+
Assert.assertEquals(cmp("", "a"), -1);
122+
Assert.assertEquals(cmp("a", ""), 1);
123+
}
124+
125+
@Test
126+
public void comparisonIsAntisymmetric() {
127+
final List<String> names = CollectionUtil.makeList(
128+
"abc", "abc03", "abc8", "abc08", "read2", "read10", "read18446744073709551616", "", "abc17.2");
129+
for (final String a : names) {
130+
for (final String b : names) {
131+
Assert.assertEquals(cmp(a, b), -cmp(b, a), "antisymmetry violated for '" + a + "' vs '" + b + "'");
132+
}
133+
}
134+
}
135+
136+
@Test
137+
public void recordLevelComparisonUsesNaturalReadNameOrder() {
138+
// Confirm the natural ordering is wired through compare(SAMRecord, SAMRecord) via fileOrderCompare.
139+
final SAMFileHeader header = new SAMFileHeader();
140+
final SAMRecord read2 = unmapped(header, "read2");
141+
final SAMRecord read10 = unmapped(header, "read10");
142+
143+
Assert.assertTrue(COMPARATOR.compare(read2, read10) < 0);
144+
Assert.assertTrue(COMPARATOR.compare(read10, read2) > 0);
145+
Assert.assertEquals(COMPARATOR.fileOrderCompare(read2, read10), -1);
146+
}
147+
148+
@Test
149+
public void inheritedPairTieBreakStillApplies() {
150+
// With identical read names, the inherited queryname tie-breaking still orders first-of-pair
151+
// before second-of-pair.
152+
final SAMFileHeader header = new SAMFileHeader();
153+
final SAMRecord first = unmapped(header, "q");
154+
first.setReadPairedFlag(true);
155+
first.setFirstOfPairFlag(true);
156+
final SAMRecord second = unmapped(header, "q");
157+
second.setReadPairedFlag(true);
158+
second.setSecondOfPairFlag(true);
159+
160+
Assert.assertEquals(COMPARATOR.fileOrderCompare(first, second), 0, "same read name compares equal");
161+
Assert.assertTrue(COMPARATOR.compare(first, second) < 0, "first of pair sorts before second of pair");
162+
Assert.assertTrue(COMPARATOR.compare(second, first) > 0);
163+
}
164+
165+
private static SAMRecord unmapped(final SAMFileHeader header, final String name) {
166+
final SAMRecord record = new SAMRecord(header);
167+
record.setReadName(name);
168+
record.setReadUnmappedFlag(true);
169+
return record;
170+
}
171+
}

0 commit comments

Comments
 (0)