Store PNEXT and TLEN as differences from what the record predicts

CHANGELOG.md

@@ -55,6 +55,13 @@ the README for what that means before 1.0.
   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]).
+- PNEXT and TLEN are stored as the difference from what POS, the strands and
+  the CIGAR predict, in the columns `pnext_delta` and `tlen_delta`, and still
+  come back exactly. On chr20 of NA12878 at 30x the file goes from 563.1 to
+  531.0 MB (-5.7%), and conversion and dumps take the same time. A program
+  that reads the columns directly gets the SAM values with
+  `RAMNTupleRecord::GetPNEXT()` and `GetTLEN()`; see CONTRIBUTING. The format
+  version becomes 6 ([#113]).
 
 ### Removed
 
@@ -74,6 +81,10 @@ the README for what that means before 1.0.
   `GetTags()`, `GetNOPT()` and `GetOPT()`. `tags` now holds only the tags
   without a column of their own. Converters split the tags with `RAMTagWriter`,
   and `RAMFile::AppendTags()` writes all of them back ([#111]).
+- `RAMNTupleRecord`'s `SetPNEXT()` and `SetTLEN()`, and `GetPNEXT()` and
+  `GetTLEN()` without arguments. Converters set both with `SetMate()`, and
+  readers pass `RAMFile::MateOnSameRef()` to `GetPNEXT()` and `GetTLEN()`
+  ([#113]).
 
 ## [0.1.0] - 2026-10-14
 
@@ -222,3 +233,4 @@ it stay readable by every later release.
 [#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
+[#113]: https://github.com/compiler-research/ramtools/pull/113

CONTRIBUTING.md

@@ -138,3 +138,16 @@ A new typed tag needs a member in `RAMNTupleRecord`, an entry in
 `kRAMTypedTags`, its handling in `RAMTagWriter::Split()` and
 `RAMTagReader::Append()`, and a view in `RecordReader` (`tools/ramdump.cxx`),
 which reads the record field by field and does not see new members by itself.
+
+## PNEXT and TLEN
+
+A record does not store PNEXT and TLEN as SAM has them, because other columns
+already say most of what they hold. When the mate is on the record's own
+reference, `pnext_delta` is PNEXT minus POS; otherwise it is PNEXT.
+`tlen_delta` is TLEN minus `PredictTLEN()`, the distance between the 5' ends
+of the pair as BWA and Picard compute it, with the mate's span taken to be the
+record's own. Both are differences modulo 2^32, so every value comes back
+exactly; for most records `tlen_delta` is 0. Write them with `SetMate()` and
+read them with `GetPNEXT()` and `GetTLEN()`, which take
+`RAMFile::MateOnSameRef()`. A program that reads the columns directly, such as
+an RDataFrame, has to do the same.

inc/rntuple/RAMFile.h

@@ -31,6 +31,12 @@ class RAMFile {
    [[nodiscard]] const std::string &GetRnext(int32_t refnext) const { return Name(m_meta.rnext_refs, refnext); }
    /// The refid of \p rname; -1 when the file has no such reference.
    [[nodiscard]] int32_t FindRefId(const std::string &rname) const;
+   /// Whether \p rec's mate is on \p rec's reference, which GetPNEXT() and
+   /// GetTLEN() of the record need.
+   [[nodiscard]] bool MateOnSameRef(const RAMNTupleRecord &rec) const
+   {
+      return RAMMateOnSameRef(GetRname(rec.refid), GetRnext(rec.refnext));
+   }
    /// Appends \p rec's optional tags to \p out as SAM text, each after a tab,
    /// in the order the input had them.
    void AppendTags(const RAMNTupleRecord &rec, std::string &out) const { m_tags.Append(rec, out); }

inc/rntuple/RAMNTupleRecord.h

@@ -13,6 +13,7 @@
 #include <ROOT/RNTupleView.hxx>
 
 #include <string>
+#include <string_view>
 #include <vector>
 #include <memory>
 #include <mutex>
@@ -166,7 +167,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 = 5;
+   static constexpr uint32_t kFormatVersion = 6;
 
    // Alignment data fields
    std::string qname;             // Query template NAME
@@ -176,8 +177,10 @@ class RAMNTupleRecord {
    uint8_t mapq;                  // MAPing Quality
    std::vector<uint32_t> cigar;   // CIGAR operations
    int32_t refnext;               // Reference ID of the mate/next read
-   int32_t pnext;                 // 0-based position of the mate/next read
-   int32_t tlen;                  // Observed Template LENgth
+   // PNEXT and TLEN, stored so that they compress; SetMate() fills them and
+   // GetPNEXT() and GetTLEN() give the SAM values back.
+   int32_t pnext_delta = 0;       // 0-based PNEXT - pos when the mate is on this reference, else 0-based PNEXT
+   int32_t tlen_delta = 0;        // TLEN - PredictTLEN()
    std::string seq;               // Segment sequence as uppercase IUPAC text, or "*"
    std::string qual;              // Quality scores (encoded)
 
@@ -207,8 +210,10 @@ class RAMNTupleRecord {
    void SetPOS(int32_t pos_) { pos = pos_ - 1; } // SAM is 1-based, we are 0-based
    void SetMAPQ(uint8_t mapq_) { mapq = mapq_; }
    void SetCIGAR(const std::string &cigar_str);
-   void SetPNEXT(int32_t pnext_) { pnext = pnext_ - 1; } // SAM is 1-based, we are 0-based
-   void SetTLEN(int32_t tlen_) { tlen = tlen_; }
+   /// Sets PNEXT (1-based) and TLEN, after SetFLAG(), SetPOS() and SetCIGAR():
+   /// they are stored against pos, the strands and the CIGAR. \p sameRef is
+   /// RAMMateOnSameRef() of the record's RNAME and RNEXT.
+   void SetMate(int32_t pnext, int32_t tlen, bool sameRef);
    void SetSEQ(const std::string &seq_str);
    void SetQUAL(const std::string &qual_str);
 
@@ -220,8 +225,14 @@ class RAMNTupleRecord {
    uint8_t GetMAPQ() const { return mapq; }
    std::string GetCIGAR() const;
    int32_t GetREFNEXT() const { return refnext; }
-   int32_t GetPNEXT() const { return pnext + 1; } // Convert back to 1-based for SAM
-   int32_t GetTLEN() const { return tlen; }
+   /// PNEXT (1-based) and TLEN as SetMate() was given them; \p sameRef as
+   /// there, or RAMFile::MateOnSameRef().
+   [[nodiscard]] int32_t GetPNEXT(bool sameRef) const;
+   [[nodiscard]] int32_t GetTLEN(bool sameRef) const;
+   /// The TLEN of a pair whose mate has the same reference span as this
+   /// record: the distance between their 5' ends, as BWA and Picard compute
+   /// it. 0 when either is unmapped or the mate is on another reference.
+   [[nodiscard]] int32_t PredictTLEN(int32_t pnext0, bool sameRef) const;
    [[nodiscard]] const std::string &GetSEQ() const { return seq; }
    std::string GetQUAL() const;
 
@@ -251,6 +262,13 @@ class RAMNTupleRecord {
    void SetCompressionMode(uint32_t flags) { compression_flags = flags; }
 };
 
+/// Whether a record's mate is on the record's own reference: RNEXT is "=" or
+/// RNAME, and RNAME is not "*". Writers and readers decide it the same way.
+inline bool RAMMateOnSameRef(std::string_view rname, std::string_view rnext)
+{
+   return rname != "*" && (rnext == "=" || rnext == rname);
+}
+
 // CIGAR operation codes (from BAM format)
 #include "ramcore/CigarOps.h"
 

src/ramcore/BamtoNTuple.cxx

@@ -183,8 +183,8 @@ void FillRecord(RAMNTupleRecord *rec, const bam1_t *b, const sam_hdr_t *hdr, uin
    rec->SetMAPQ(b->core.qual);
    rec->SetCIGAR(GetCigar(b));
    rec->refnext = rnexts.GetRefId(rnext);
-   rec->SetPNEXT(static_cast<int32_t>(b->core.mpos + 1));
-   rec->SetTLEN(static_cast<int32_t>(b->core.isize));
+   rec->SetMate(static_cast<int32_t>(b->core.mpos + 1), static_cast<int32_t>(b->core.isize),
+                RAMMateOnSameRef(rname, rnext));
    rec->SetSEQ(GetSeq(b));
    rec->SetQUAL(GetQual(b));
    tags.Split(GetTags(b), *rec);

src/ramcore/SamToNTuple.cxx

@@ -60,8 +60,7 @@ void FillRecordFields(const ramcore::SamRecord &sam_record, RAMNTupleRecord &rec
    rec.SetPOS(sam_record.pos);
    rec.SetMAPQ(sam_record.mapq);
    rec.SetCIGAR(sam_record.cigar);
-   rec.SetPNEXT(sam_record.pnext);
-   rec.SetTLEN(sam_record.tlen);
+   rec.SetMate(sam_record.pnext, sam_record.tlen, RAMMateOnSameRef(sam_record.rname, sam_record.rnext));
    rec.SetSEQ(sam_record.seq);
    rec.SetQUAL(sam_record.qual);
    tags.Split(sam_record.optional_fields, rec);

src/rntuple/RAMNTupleRecord.cxx

@@ -139,9 +139,45 @@ void RAMNTupleRefs::Print() const
 }
 // RAMNTupleRecord Implementation
 
-RAMNTupleRecord::RAMNTupleRecord()
-   : flag(0), refid(-1), pos(0), mapq(0), refnext(-1), pnext(0), tlen(0), compression_flags(kPhred33)
+RAMNTupleRecord::RAMNTupleRecord() : flag(0), refid(-1), pos(0), mapq(0), refnext(-1), compression_flags(kPhred33) {}
+
+// pnext_delta and tlen_delta are differences modulo 2^32, so every int32 comes
+// back exactly.
+static int32_t Wrap(int64_t value)
+{
+   return static_cast<int32_t>(static_cast<uint32_t>(value));
+}
+
+void RAMNTupleRecord::SetMate(int32_t pnext, int32_t tlen, bool sameRef)
+{
+   const int32_t pnext0 = pnext - 1; // SAM is 1-based, we are 0-based
+   pnext_delta = sameRef ? Wrap(int64_t{pnext0} - pos) : pnext0;
+   tlen_delta = Wrap(int64_t{tlen} - PredictTLEN(pnext0, sameRef));
+}
+
+int32_t RAMNTupleRecord::GetPNEXT(bool sameRef) const
 {
+   return (sameRef ? Wrap(int64_t{pnext_delta} + pos) : pnext_delta) + 1;
+}
+
+int32_t RAMNTupleRecord::GetTLEN(bool sameRef) const
+{
+   return Wrap(int64_t{tlen_delta} + PredictTLEN(GetPNEXT(sameRef) - 1, sameRef));
+}
+
+int32_t RAMNTupleRecord::PredictTLEN(int32_t pnext0, bool sameRef) const
+{
+   const int64_t span = GetRefSpan();
+   if (!sameRef || (flag & 0x4) || (flag & 0x8) || span == 0)
+      return 0;
+   const int64_t own5 = (flag & 0x10) ? pos + span - 1 : pos;
+   const int64_t mate5 = (flag & 0x20) ? pnext0 + span - 1 : pnext0;
+   const int64_t distance = mate5 - own5;
+   if (distance > 0)
+      return Wrap(distance + 1);
+   if (distance < 0)
+      return Wrap(distance - 1);
+   return 0;
 }
 
 // RAMFile::Open() reads these fields back.

test/bam_conversion.cxx

@@ -286,6 +286,28 @@ TEST_F(BamToNTupleTest, TagsOfEveryTypeComeBackAsSamtoolsPrintsThem)
                    "\tXA:A:Q\tXZ:Z:hello\tXH:H:1AE301");
 }
 
+// A mate on the same reference and one on another come back as in the BAM.
+TEST_F(BamToNTupleTest, MateFieldsComeBackAsInTheBam)
+{
+   GenerateRichBAMFile(/*bam_path=*/"test_rich.bam");
+   bamtoramntuple(/*bamfile=*/"test_rich.bam", /*treefile=*/"test_rich.ram", /*split=*/false, /*cache=*/true,
+                  /*compression_algorithm=*/505, /*quality_policy=*/0U);
+
+   const auto file = RAMFile::Open("test_rich.ram");
+   ASSERT_NE(file, nullptr);
+   auto view = file->GetReader().GetView<RAMNTupleRecord>("record");
+   std::string mates{};
+   for (auto row : file->GetReader().GetEntryRange()) {
+      const RAMNTupleRecord &rec = view(row);
+      if (rec.GetQNAME() == "read_samechrom" || rec.GetQNAME() == "read_diffchrom") {
+         const bool sameRef = file->MateOnSameRef(rec);
+         mates += rec.GetQNAME() + " " + file->GetRnext(rec.GetREFNEXT()) + " " +
+                  std::to_string(rec.GetPNEXT(sameRef)) + " " + std::to_string(rec.GetTLEN(sameRef)) + "\n";
+      }
+   }
+   EXPECT_EQ(mates, "read_diffchrom chr2 2001 0\nread_samechrom = 801 200\n");
+}
+
 // As samtools idxstats counts them: every record of the rich BAM is on chr1
 // except read_unmapped, which has no reference.
 TEST_F(BamToNTupleTest, RecordsAreCountedPerReference)

test/parallel_conversion_test.cxx

@@ -89,10 +89,11 @@ Dump ReadBack(const char *file)
    QualityBlockReader quals(*ram);
    for (auto i : reader.GetEntryRange()) {
       const RAMNTupleRecord &r = view(i);
+      const bool sameRef = ram->MateOnSameRef(r);
       std::string line = r.GetQNAME() + "\t" + std::to_string(r.GetFLAG()) + "\t" + ram->GetRname(r.GetREFID()) + "\t" +
                          std::to_string(r.GetPOS()) + "\t" + std::to_string(r.GetMAPQ()) + "\t" + r.GetCIGAR() + "\t" +
-                         ram->GetRnext(r.GetREFNEXT()) + "\t" + std::to_string(r.GetPNEXT()) + "\t" +
-                         std::to_string(r.GetTLEN()) + "\t" + r.GetSEQ() + "\t" + quals.Get(r, i);
+                         ram->GetRnext(r.GetREFNEXT()) + "\t" + std::to_string(r.GetPNEXT(sameRef)) + "\t" +
+                         std::to_string(r.GetTLEN(sameRef)) + "\t" + r.GetSEQ() + "\t" + quals.Get(r, i);
       ram->AppendTags(r, line);
       d.lines.push_back(line);
    }

test/ramcoretests.cxx

@@ -725,6 +725,52 @@ TEST_F(ramcoreTest, TagsOfADamagedRecordStillComeOut)
    EXPECT_EQ(out, "\tXA:Z:chr1,+100,150M,0;\tpa:f:0.5");
 }
 
+// PNEXT and TLEN come back exactly, however far TLEN is from what pos, the
+// strands and the CIGAR predict.
+TEST_F(ramcoreTest, MateFieldsComeBackExactly)
+{
+   struct Mate {
+      uint16_t flag;
+      int32_t pos; // 1-based, as in SAM
+      const char *cigar;
+      int32_t pnext;
+      int32_t tlen;
+      bool sameRef;
+   };
+   const std::vector<Mate> mates = {
+      {0x63, 1001, "150M", 1301, 450, true},          // forward, the mate reverse
+      {0x93, 1301, "150M", 1001, -450, true},         // that mate
+      {0x63, 1001, "150M", 1301, 448, true},          // the mate's span is 148
+      {0x63, 1001, "150M", 852, 1, true},             // the same 5' end: Picard says 1
+      {0x41, 1001, "150M", 5001, 0, false},           // the mate on another reference
+      {0x49, 1001, "150M", 1001, 0, true},            // the mate unmapped
+      {0x85, 1001, "*", 1001, 0, true},               // unmapped, placed at its mate
+      {0, 1, "4M", 0, 0, false},                      // no mate
+      {0x63, 1, "4M", 2147483647, -2147483647, true}, // the extremes
+      {0x93, 2147483647, "4M", 1, 2147483647, true},
+   };
+   for (std::size_t i = 0; i < mates.size(); i++) {
+      const Mate &m = mates[i];
+      RAMNTupleRecord rec;
+      rec.SetFLAG(m.flag);
+      rec.SetPOS(m.pos);
+      rec.SetCIGAR(m.cigar);
+      rec.SetMate(m.pnext, m.tlen, m.sameRef);
+      EXPECT_EQ(rec.GetPNEXT(m.sameRef), m.pnext) << "mate " << i;
+      EXPECT_EQ(rec.GetTLEN(m.sameRef), m.tlen) << "mate " << i;
+      if (i < 3) {
+         EXPECT_EQ(rec.pnext_delta, m.flag & 0x10 ? -300 : 300);
+         EXPECT_EQ(rec.tlen_delta, i < 2 ? 0 : -2) << "mate " << i;
+      }
+   }
+
+   EXPECT_TRUE(RAMMateOnSameRef("chr1", "="));
+   EXPECT_TRUE(RAMMateOnSameRef("chr1", "chr1"));
+   EXPECT_FALSE(RAMMateOnSameRef("chr1", "chr2"));
+   EXPECT_FALSE(RAMMateOnSameRef("chr1", "*"));
+   EXPECT_FALSE(RAMMateOnSameRef("*", "="));
+}
+
 // A converter's table: "*" is no reference, "=" gets an id like any name.
 TEST_F(ramcoreTest, RefTableMapsNamesToIds)
 {

test/ramdump_test.cxx

@@ -138,6 +138,30 @@ TEST_F(RamdumpTest, TagsComeBackAsTheyWere)
    std::remove(ram);
 }
 
+// PNEXT and TLEN are stored against pos and the CIGAR; every kind of mate
+// comes back as it was, a TLEN no aligner would write among them.
+TEST_F(RamdumpTest, MateFieldsComeBackAsTheyWere)
+{
+   const char *sam = "ramdump_test_mates.sam";
+   const char *ram = "ramdump_test_mates.root";
+   std::string records{};
+   records += "m1\t99\tchr1\t100\t60\t4M\t=\t300\t204\tACGT\tIIII\n";
+   records += "m1\t147\tchr1\t300\t60\t4M\t=\t100\t-204\tACGT\tIIII\n";
+   records += "m2\t99\tchr1\t400\t60\t2S2M\t=\t500\t103\tACGT\tIIII\n";
+   records += "m3\t99\tchr1\t600\t60\t4M\tchr1\t700\t104\tACGT\tIIII\n";
+   records += "m4\t65\tchr1\t800\t60\t4M\tchr2\t10\t0\tACGT\tIIII\n";
+   records += "m5\t73\tchr1\t900\t60\t4M\t=\t900\t0\tACGT\tIIII\n";
+   records += "m5\t133\tchr1\t900\t0\t*\t=\t900\t0\tACGT\tIIII\n";
+   records += "m6\t0\tchr2\t50\t60\t4M\t*\t0\t0\tACGT\tIIII\n";
+   records += "m7\t99\tchr2\t100\t60\t4M\t=\t200\t-5000\tACGT\tIIII\n";
+   records += "m8\t4\t*\t0\t0\t*\t*\t0\t0\tACGT\tIIII\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);
+   std::remove(sam);
+   std::remove(ram);
+}
+
 TEST_F(RamdumpTest, HeaderOnly)
 {
    EXPECT_EQ(Dump("-H"), Header());

tools/ramdump.cxx

@@ -77,8 +77,8 @@ class RecordReader {
    ROOT::RNTupleView<uint8_t> m_mapq;
    ROOT::RNTupleView<std::vector<uint32_t>> m_cigar;
    ROOT::RNTupleView<int32_t> m_refnext;
-   ROOT::RNTupleView<int32_t> m_pnext;
-   ROOT::RNTupleView<int32_t> m_tlen;
+   ROOT::RNTupleView<int32_t> m_pnextDelta;
+   ROOT::RNTupleView<int32_t> m_tlenDelta;
    ROOT::RNTupleView<std::string> m_seq;
    ROOT::RNTupleView<std::string> m_qual;
    ROOT::RNTupleView<int32_t> m_tagNm;
@@ -104,8 +104,8 @@ class RecordReader {
         m_mapq(reader.GetView("record.mapq", &m_record.mapq)),
         m_cigar(reader.GetView("record.cigar", &m_record.cigar)),
         m_refnext(reader.GetView("record.refnext", &m_record.refnext)),
-        m_pnext(reader.GetView("record.pnext", &m_record.pnext)),
-        m_tlen(reader.GetView("record.tlen", &m_record.tlen)),
+        m_pnextDelta(reader.GetView("record.pnext_delta", &m_record.pnext_delta)),
+        m_tlenDelta(reader.GetView("record.tlen_delta", &m_record.tlen_delta)),
         m_seq(reader.GetView("record.seq", &m_record.seq)),
         m_qual(reader.GetView("record.qual", &m_record.qual)),
         m_tagNm(reader.GetView("record.tag_nm", &m_record.tag_nm)),
@@ -139,8 +139,8 @@ class RecordReader {
       m_mapq(row);
       m_cigar(row);
       m_refnext(row);
-      m_pnext(row);
-      m_tlen(row);
+      m_pnextDelta(row);
+      m_tlenDelta(row);
       m_seq(row);
       m_qual(row);
       m_tagNm(row);
@@ -176,9 +176,10 @@ void WriteRecord(const RAMFile &file, const RAMNTupleRecord &rec, const std::str
    out.Tab();
    out.Str(file.GetRnext(rec.GetREFNEXT()));
    out.Tab();
-   out.Int(rec.GetPNEXT());
+   const bool sameRef = file.MateOnSameRef(rec);
+   out.Int(rec.GetPNEXT(sameRef));
    out.Tab();
-   out.Int(rec.GetTLEN());
+   out.Int(rec.GetTLEN(sameRef));
    out.Tab();
    out.Str(rec.GetSEQ());
    out.Tab();