Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
43 changes: 30 additions & 13 deletions src/rntuple/RAMNTupleRecord.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -460,20 +460,32 @@ std::string DecodeSequence(const std::string &encoded_seq, size_t length)

std::string EncodeQuality(const std::string &qual, uint32_t compression_flags)
{
std::string encoded;

if (compression_flags & RAMNTupleRecord::kDrop) {
encoded = "*";
} else if (compression_flags & RAMNTupleRecord::kIlluminaBinning) {
encoded.resize(qual.size());
return "*";
}

if (compression_flags & RAMNTupleRecord::kIlluminaBinning) {
// "*" means "quality not available". It is a sentinel, not a Phred
// string, so it must not be fed through the binning table
if (qual == "*")
return {};

std::string encoded(qual.size(), '\0');
for (size_t i = 0; i < qual.size(); i++) {
encoded[i] = kIlluminaBinning[static_cast<uint8_t>(qual[i])];
// SAM stores quality as Phred+33 ASCII, but kIlluminaBinning is
// indexed by the Phred VALUE. Without the -33 every lookup lands 33
// slots too far right
// Clamp to 0..93: 93 is SAM's maximum Phred, and it also keeps the
// index inside the initialised part of the table (entries 110..255
// are zero-filled, so an out-of-range value would silently decode as
// Q0 -- the opposite error, but still an error).
const int phred = std::clamp(static_cast<int>(static_cast<unsigned char>(qual[i])) - 33, 0, 93);
encoded[i] = static_cast<char>(kIlluminaBinning[phred]);
}
} else {
encoded = qual;
return encoded;
}

return encoded;
return qual;
}

std::string DecodeQuality(const std::string &encoded_qual, uint32_t compression_flags)
Expand All @@ -482,15 +494,20 @@ std::string DecodeQuality(const std::string &encoded_qual, uint32_t compression_
return "*";
}

std::string qual = encoded_qual;

if (compression_flags & RAMNTupleRecord::kIlluminaBinning) {
// Empty is the "quality not available" sentinel written by
// EncodeQuality; restore the "*" it stood for.
if (encoded_qual.empty())
return "*";

std::string qual = encoded_qual;
for (auto &q : qual) {
q += 33;
q = static_cast<char>(static_cast<unsigned char>(q) + 33);
}
return qual;
}

return qual;
return encoded_qual;
}

std::vector<uint32_t> ParseCIGAR(const std::string &cigar_str)
Expand Down
87 changes: 76 additions & 11 deletions test/ramcoretests.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include <Rtypes.h>
#include <TFile.h>
#include <TTree.h>
#include <array>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
Expand Down Expand Up @@ -323,19 +324,83 @@ TEST_F(ramcoreTest, RecordGetters)
dropRecord.SetQUAL("IIIII");
EXPECT_EQ(dropRecord.GetQUAL(), "*");

// kIlluminaBinning, ASCII bins 0, 1, 6, 15, 22, 27, 33, 37, 40
// kIlluminaBinning maps a Phred VALUE to one of 0,1,6,15,22,27,33,37,40.
// SAM writes quality as Phred+33 ASCII, so the encoder must subtract 33
// before the lookup. These expectations are stated in Phred space and
// converted, so they cannot silently drift back to indexing by ASCII.
RAMNTupleRecord binRecord;
binRecord.SetBit(RAMNTupleRecord::kIlluminaBinning);
binRecord.SetQUAL("\""); // ASCII 34 → bin 33 → 'B'
EXPECT_EQ(binRecord.GetQUAL(), "B");
binRecord.SetQUAL("$"); // ASCII 36 → bin 37 → 'F'
EXPECT_EQ(binRecord.GetQUAL(), "F");
binRecord.SetQUAL("'"); // ASCII 39 → bin 37 → 'F'
EXPECT_EQ(binRecord.GetQUAL(), "F");
binRecord.SetQUAL("("); // ASCII 40 → bin 40 → 'I'
EXPECT_EQ(binRecord.GetQUAL(), "I");
binRecord.SetQUAL("2"); // ASCII 50 → bin 40 → 'I'
EXPECT_EQ(binRecord.GetQUAL(), "I");

// Two characters, not one: Phred 9 encodes to ASCII 42, which is '*'. A
// one-base read whose quality is "*" is ambiguous in SAM itself (sentinel
// vs. Q9), so single-character quality strings are a bad test vector.
auto phred = [](int q) { return std::string(2, static_cast<char>(q + 33)); };
auto roundTrip = [&](int q) {
binRecord.SetQUAL(phred(q));
const std::string out = binRecord.GetQUAL();
EXPECT_EQ(out.size(), 2U);
return static_cast<int>(static_cast<unsigned char>(out[0])) - 33;
};

struct Bin {
int in;
int want;
};
const std::array<Bin, 16> kBins = {{
{0, 0},
{1, 1},
{2, 6},
{9, 6},
{10, 15},
{19, 15},
{20, 22},
{24, 22},
{25, 27},
{29, 27},
{30, 33},
{34, 33},
{35, 37},
{39, 37},
{40, 40},
{93, 40},
}};
for (const auto &c : kBins)
EXPECT_EQ(roundTrip(c.in), c.want) << "Q" << c.in << " should bin to Q" << c.want;

// NOTE: binning is NOT monotonically downward -- Illumina maps each bin to
// a representative value near its middle, so Q2 legitimately becomes Q6.
// "never raises a quality" is therefore the wrong invariant. What the
// ASCII-indexing bug actually violated is captured below.

// Q0 means "no usable base". It must survive as Q0; the old code rewrote it
// as Q33 (0.05% error), fabricating confidence a variant caller would trust.
EXPECT_EQ(roundTrip(0), 0) << "a zero-quality base must not be upgraded";

int previous = -1;
for (int q = 0; q <= 93; ++q) {
const int got = roundTrip(q);
EXPECT_GE(got, previous) << "binning must be monotonic; broke at Q" << q;
previous = got;
EXPECT_TRUE(got == 0 || got == 1 || got == 6 || got == 15 || got == 22 || got == 27 || got == 33 || got == 37 ||
got == 40)
<< "Q" << q << " produced Q" << got << ", not a legal Illumina bin";
}

// "*" means "no quality available"; it is a sentinel, not a Phred string,
// and must survive rather than be run through the table. Feeding it through
// mapped '*' (ASCII 42) to bin 40, producing a one-character QUAL against a
// full-length SEQ -- a malformed record, not merely a wrong one.
binRecord.SetQUAL("*");
EXPECT_EQ(binRecord.GetQUAL(), "*");

// The lossless path must pass the sentinel through untouched too.
RAMNTupleRecord losslessRecord;
losslessRecord.SetQUAL("*");
EXPECT_EQ(losslessRecord.GetQUAL(), "*");

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: integer literal has suffix 'u', which is not uppercase [readability-uppercase-literal-suffix]

Suggested change
EXPECT_EQ(binRecord.GetQUAL().size(), 10U);

// Length must be preserved for real quality strings.
binRecord.SetQUAL("IIIIIIIIII");
EXPECT_EQ(binRecord.GetQUAL().size(), 10U);
}

} // namespace
Expand Down
Loading