Store MD, MC, SA and XA in columns of their own

CHANGELOG.md

@@ -48,12 +48,13 @@ the README for what that means before 1.0.
   reference tables and the METADATA values, in their own objects instead of in
   static members of `RAMNTupleRecord`, so several conversions can run in one
   process at once. Before, two at once mixed their reference names ([#107]).
-- The tags NM, AS, XS, MQ, RG and PG are stored in columns of their own instead
-  of as text, RG and PG as ids into a dictionary in METADATA. The other tags
-  stay text, and `ramdump -h` still reproduces the input byte for byte. On chr20
-  of NA12878 at 30x the file is 7.3% smaller (622.5 to 576.8 MB), conversion
-  takes 5 to 13% less CPU, and dumps take the same time. The format version
-  becomes 4, and files written before must be converted again ([#111]).
+- The tags NM, AS, XS, MQ, RG, PG, MD, MC, SA and XA are stored in columns of
+  their own instead of as text, RG and PG as ids into a dictionary in METADATA.
+  The other tags stay text, and `ramdump -h` still reproduces the input byte for
+  byte. On chr20 of NA12878 at 30x the file is 9.5% smaller (622.5 to
+  563.1 MB), conversion takes 5 to 13% less CPU, and `ramdump -h` 8% less. The
+  format version becomes 5, and files written before must be converted again
+  ([#111], [#112]).
 
 ### Removed
 
@@ -220,3 +221,4 @@ it stay readable by every later release.
 [#107]: https://github.com/compiler-research/ramtools/pull/107
 [#108]: https://github.com/compiler-research/ramtools/pull/108
 [#111]: https://github.com/compiler-research/ramtools/pull/111
+[#112]: https://github.com/compiler-research/ramtools/pull/112

inc/rntuple/RAMNTupleRecord.h

@@ -166,7 +166,7 @@ class RAMNTupleRecord {
 
    /// Version of the file layout, stored in METADATA. Until 1.0 a reader opens
    /// only files of its own version; bump it with every change to the layout.
-   static constexpr uint32_t kFormatVersion = 4;
+   static constexpr uint32_t kFormatVersion = 5;
 
    // Alignment data fields
    std::string qname;             // Query template NAME
@@ -189,6 +189,10 @@ class RAMNTupleRecord {
    int32_t tag_mq = 0;
    uint32_t tag_rg = 0; // Id into METADATA's tag_rg
    uint32_t tag_pg = 0; // Id into METADATA's tag_pg
+   std::string tag_md;
+   std::string tag_mc;
+   std::string tag_sa;
+   std::string tag_xa;
    uint32_t tag_order = 0;
    std::vector<std::string> tags; // The other optional tags, as SAM text
 

inc/rntuple/RAMTags.h

@@ -17,7 +17,8 @@
 // then RG, then NM. The order also says which typed tags the record has.
 
 /// The typed tags, in the order of RAMNTupleRecord's tag_ members.
-inline constexpr std::array<const char *, 6> kRAMTypedTags = {"NM", "AS", "XS", "MQ", "RG", "PG"};
+inline constexpr std::array<const char *, 10> kRAMTypedTags = {"NM", "AS", "XS", "MQ", "RG",
+                                                               "PG", "MD", "MC", "SA", "XA"};
 /// Marks a text tag in a tag order.
 inline constexpr const char *kRAMTextTag = "**";
 
@@ -76,6 +77,8 @@ class RAMTagReader {
 private:
    /// Per order, the index into kRAMTypedTags of each tag, or -1 for text.
    std::vector<std::vector<int8_t>> m_orders;
+   /// "\tNM:i:" and the rest, by index into kRAMTypedTags.
+   std::array<std::string, kRAMTypedTags.size()> m_prefixes;
    std::vector<std::string> m_rg;
    std::vector<std::string> m_pg;
 };

src/rntuple/RAMTags.cxx

@@ -6,6 +6,7 @@
 #include <cstddef>
 #include <cstdint>
 #include <cstdlib>
+#include <initializer_list>
 #include <limits>
 #include <mutex>
 #include <string>
@@ -18,8 +19,10 @@ namespace {
 
 constexpr int kRG = 4;
 constexpr int kPG = 5;
+constexpr int kMD = 6;
 
-// The members kRAMTypedTags names, by its index; the integer ones come first.
+// The members kRAMTypedTags names, by its index: the integer tags, then RG and
+// PG, then the string tags.
 int32_t *IntTag(RAMNTupleRecord &rec, int i)
 {
    const std::array<int32_t *, 4> members = {&rec.tag_nm, &rec.tag_as, &rec.tag_xs, &rec.tag_mq};
@@ -32,6 +35,18 @@ int32_t IntTag(const RAMNTupleRecord &rec, int i)
    return members.at(static_cast<std::size_t>(i));
 }
 
+std::string *StringTag(RAMNTupleRecord &rec, int i)
+{
+   const std::array<std::string *, 4> members = {&rec.tag_md, &rec.tag_mc, &rec.tag_sa, &rec.tag_xa};
+   return members.at(static_cast<std::size_t>(i - kMD));
+}
+
+const std::string &StringTag(const RAMNTupleRecord &rec, int i)
+{
+   const std::array<const std::string *, 4> members = {&rec.tag_md, &rec.tag_mc, &rec.tag_sa, &rec.tag_xa};
+   return *members.at(static_cast<std::size_t>(i - kMD));
+}
+
 // The index into kRAMTypedTags of the tag named \p name, or -1.
 int NameIndex(std::string_view name) noexcept
 {
@@ -117,6 +132,8 @@ void RAMTagWriter::Split(const std::vector<std::string> &tags, RAMNTupleRecord &
 {
    rec.tag_nm = rec.tag_as = rec.tag_xs = rec.tag_mq = 0;
    rec.tag_rg = rec.tag_pg = 0;
+   for (std::string *value : {&rec.tag_md, &rec.tag_mc, &rec.tag_sa, &rec.tag_xa})
+      value->clear();
    rec.tags.clear();
    m_order.clear();
    std::array<bool, kRAMTypedTags.size()> present{};
@@ -125,6 +142,8 @@ void RAMTagWriter::Split(const std::vector<std::string> &tags, RAMNTupleRecord &
       bool typed = i >= 0 && !present.at(static_cast<std::size_t>(i));
       if (typed && i < kRG) {
          typed = CanonicalInt(tag, *IntTag(rec, i));
+      } else if (typed && i >= kMD) {
+         StringTag(rec, i)->assign(tag, 5);
       } else if (typed) {
          m_value.assign(tag, 5);
          if (i == kRG)
@@ -151,6 +170,8 @@ RAMTagReader::RAMTagReader(const RAMFileMetadata &meta) : m_rg(meta.tag_rg), m_p
          tokens.push_back(static_cast<int8_t>(NameIndex(std::string_view(&order[at], 2))));
       m_orders.push_back(std::move(tokens));
    }
+   for (std::size_t i = 0; i < kRAMTypedTags.size(); i++)
+      m_prefixes.at(i) = std::string("\t") + kRAMTypedTags.at(i) + (static_cast<int>(i) < kRG ? ":i:" : ":Z:");
 }
 
 void RAMTagReader::Append(const RAMNTupleRecord &rec, std::string &out) const
@@ -165,13 +186,12 @@ void RAMTagReader::Append(const RAMNTupleRecord &rec, std::string &out) const
             }
             continue;
          }
-         out += '\t';
-         out += kRAMTypedTags.at(static_cast<std::size_t>(i));
+         out += m_prefixes.at(static_cast<std::size_t>(i));
          if (i < kRG) {
-            out += ":i:";
             out += std::to_string(IntTag(rec, i));
+         } else if (i >= kMD) {
+            out += StringTag(rec, i);
          } else {
-            out += ":Z:";
             out += Value(i == kRG ? m_rg : m_pg, i == kRG ? rec.tag_rg : rec.tag_pg);
          }
       }

test/ramcoretests.cxx

@@ -644,13 +644,15 @@ TEST_F(ramcoreTest, TypedTagsComeBackInTheirOrder)
 {
    const std::vector<std::vector<std::string>> inputs = {
       {"MC:Z:150M", "PG:Z:MarkDuplicates", "MQ:i:60", "AS:i:148", "XS:i:0", "MD:Z:150", "NM:i:0", "RG:Z:lane1"},
-      {"RG:Z:lane2", "AS:i:-3", "pa:f:0.5"}, // no NM; an unknown tag
+      // no NM; an unknown tag
+      {"RG:Z:lane2", "SA:Z:chr2,500,+,100M50S,60,0;", "AS:i:-3", "pa:f:0.5", "XA:Z:chr1,+100,150M,0;"},
       {},
       {"NM:i:007", "AS:i:+5", "XS:i:4294967295", "MQ:i:-0"},                 // not written back the same as numbers
-      {"NM:Z:x", "RG:i:5"},                                                  // known tags with another type
+      {"NM:Z:x", "RG:i:5", "MD:i:5", "XA:A:Q"},                              // known tags with another type
       {"NM:i:1", "NM:i:2", "RG:Z:", "RG:Z:lane1"},                           // second copies; an empty value
       {"NM", "NM:i", "NM-i:1", "NM:i-1"},                                    // not "TG:T:value"
       {"NM:i:", "AS:i:5x", "XS:i:99999999999999999999", "MQ:i:-2147483649"}, // not an int32
+      {"MD:Z:1", "MD:Z:2", "SA:Z:"},                                         // the same for a string tag
    };
    RAMTagTables tables;
    RAMTagWriter writer(tables);
@@ -661,17 +663,30 @@ TEST_F(ramcoreTest, TypedTagsComeBackInTheirOrder)
       records.push_back(rec);
    }
 
-   EXPECT_EQ(records[0].tags, (std::vector<std::string>{"MC:Z:150M", "MD:Z:150"}));
+   EXPECT_TRUE(records[0].tags.empty());
    EXPECT_EQ(records[0].tag_mq, 60);
    EXPECT_EQ(records[0].tag_as, 148);
+   EXPECT_EQ(records[0].tag_md, "150");
+   EXPECT_EQ(records[0].tag_mc, "150M");
    EXPECT_EQ(records[1].tags, std::vector<std::string>{"pa:f:0.5"});
    EXPECT_EQ(records[1].tag_as, -3);
+   EXPECT_EQ(records[1].tag_sa, "chr2,500,+,100M50S,60,0;");
+   EXPECT_EQ(records[1].tag_xa, "chr1,+100,150M,0;");
    EXPECT_EQ(records[3].tags, inputs[3]);
    EXPECT_EQ(records[4].tags, inputs[4]);
    EXPECT_EQ(records[5].tags, (std::vector<std::string>{"NM:i:2", "RG:Z:lane1"}));
    EXPECT_EQ(records[5].tag_nm, 1);
    EXPECT_EQ(records[6].tags, inputs[6]);
    EXPECT_EQ(records[7].tags, inputs[7]);
+   EXPECT_EQ(records[8].tags, std::vector<std::string>{"MD:Z:2"});
+   EXPECT_EQ(records[8].tag_md, "1");
+
+   // A record the converter fills again keeps nothing of the one before.
+   RAMNTupleRecord reused;
+   writer.Split(inputs[0], reused);
+   writer.Split(inputs[2], reused);
+   EXPECT_EQ(reused.tag_mq, 0);
+   EXPECT_EQ(reused.tag_md, "");
 
    RAMFileMetadata meta;
    tables.Store(meta);

test/ramdump_test.cxx

@@ -126,10 +126,11 @@ TEST_F(RamdumpTest, TagsComeBackAsTheyWere)
    std::string records{};
    records += "t1\t0\tchr1\t100\t60\t4M\t*\t0\t0\tACGT\tIIII\tMC:Z:4M\tPG:Z:MarkDuplicates\tMQ:i:60\tAS:i:4\t"
               "XS:i:0\tMD:Z:4\tNM:i:0\tRG:Z:lane1\n";
-   records += "t2\t0\tchr1\t110\t60\t4M\t*\t0\t0\tACGT\tIIII\tRG:Z:lane2\tAS:i:-3\tpa:f:0.5\n";
+   records += "t2\t0\tchr1\t110\t60\t4M\t*\t0\t0\tACGT\tIIII\tRG:Z:lane2\tSA:Z:chr2,500,+,2M2S,60,0;\tAS:i:-3\t"
+              "pa:f:0.5\tXA:Z:chr1,+300,4M,0;\n";
    records += "t3\t0\tchr1\t120\t60\t4M\t*\t0\t0\tACGT\tIIII\n";
-   records += "t4\t0\tchr1\t130\t60\t4M\t*\t0\t0\tACGT\tIIII\tNM:i:007\tXS:i:4294967295\tNM:Z:x\tRG:i:5\n";
-   records += "t5\t0\tchr1\t140\t60\t4M\t*\t0\t0\tACGT\tIIII\tNM:i:1\tNM:i:2\tRG:Z:\tRG:Z:lane1\n";
+   records += "t4\t0\tchr1\t130\t60\t4M\t*\t0\t0\tACGT\tIIII\tNM:i:007\tXS:i:4294967295\tNM:Z:x\tRG:i:5\tMD:i:4\n";
+   records += "t5\t0\tchr1\t140\t60\t4M\t*\t0\t0\tACGT\tIIII\tNM:i:1\tNM:i:2\tRG:Z:\tRG:Z:lane1\tMD:Z:4\tMD:Z:3A0\n";
    std::ofstream(sam) << Header() << records;
    ASSERT_EQ(Quiet(std::string(SAMTORAMNTUPLE_BIN) + " " + sam + " " + ram + " -threads 2"), 0);
    EXPECT_EQ(Output(Command("-h", ram, "")), Header() + records);

tools/ramdump.cxx

@@ -87,6 +87,10 @@ class RecordReader {
    ROOT::RNTupleView<int32_t> m_tagMq;
    ROOT::RNTupleView<uint32_t> m_tagRg;
    ROOT::RNTupleView<uint32_t> m_tagPg;
+   ROOT::RNTupleView<std::string> m_tagMd;
+   ROOT::RNTupleView<std::string> m_tagMc;
+   ROOT::RNTupleView<std::string> m_tagSa;
+   ROOT::RNTupleView<std::string> m_tagXa;
    ROOT::RNTupleView<uint32_t> m_tagOrder;
    ROOT::RNTupleView<std::vector<std::string>> m_tags;
    ROOT::RNTupleView<uint32_t> m_compressionFlags;
@@ -110,6 +114,10 @@ class RecordReader {
         m_tagMq(reader.GetView("record.tag_mq", &m_record.tag_mq)),
         m_tagRg(reader.GetView("record.tag_rg", &m_record.tag_rg)),
         m_tagPg(reader.GetView("record.tag_pg", &m_record.tag_pg)),
+        m_tagMd(reader.GetView("record.tag_md", &m_record.tag_md)),
+        m_tagMc(reader.GetView("record.tag_mc", &m_record.tag_mc)),
+        m_tagSa(reader.GetView("record.tag_sa", &m_record.tag_sa)),
+        m_tagXa(reader.GetView("record.tag_xa", &m_record.tag_xa)),
         m_tagOrder(reader.GetView("record.tag_order", &m_record.tag_order)),
         m_tags(reader.GetView("record.tags", &m_record.tags)),
         m_compressionFlags(reader.GetView("record.compression_flags", &m_record.compression_flags))
@@ -141,6 +149,10 @@ class RecordReader {
       m_tagMq(row);
       m_tagRg(row);
       m_tagPg(row);
+      m_tagMd(row);
+      m_tagMc(row);
+      m_tagSa(row);
+      m_tagXa(row);
       m_tagOrder(row);
       m_tags(row);
       m_compressionFlags(row);