diff --git a/CHANGELOG.md b/CHANGELOG.md index f123fd6..11ec15b 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index d54b172..a469559 100644 --- a/CONTRIBUTING.md +++ b/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. diff --git a/inc/rntuple/RAMFile.h b/inc/rntuple/RAMFile.h index f0ebd28..3f3d1bd 100644 --- a/inc/rntuple/RAMFile.h +++ b/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); } diff --git a/inc/rntuple/RAMNTupleRecord.h b/inc/rntuple/RAMNTupleRecord.h index d9988fa..b85bf0c 100644 --- a/inc/rntuple/RAMNTupleRecord.h +++ b/inc/rntuple/RAMNTupleRecord.h @@ -13,6 +13,7 @@ #include #include +#include #include #include #include @@ -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 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" diff --git a/src/ramcore/BamtoNTuple.cxx b/src/ramcore/BamtoNTuple.cxx index 4415d74..02f6aa3 100644 --- a/src/ramcore/BamtoNTuple.cxx +++ b/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(b->core.mpos + 1)); - rec->SetTLEN(static_cast(b->core.isize)); + rec->SetMate(static_cast(b->core.mpos + 1), static_cast(b->core.isize), + RAMMateOnSameRef(rname, rnext)); rec->SetSEQ(GetSeq(b)); rec->SetQUAL(GetQual(b)); tags.Split(GetTags(b), *rec); diff --git a/src/ramcore/SamToNTuple.cxx b/src/ramcore/SamToNTuple.cxx index ae90385..7baeb44 100644 --- a/src/ramcore/SamToNTuple.cxx +++ b/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); diff --git a/src/rntuple/RAMNTupleRecord.cxx b/src/rntuple/RAMNTupleRecord.cxx index 8d7144e..b1f40af 100644 --- a/src/rntuple/RAMNTupleRecord.cxx +++ b/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(static_cast(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. diff --git a/test/bam_conversion.cxx b/test/bam_conversion.cxx index 2fb0830..30afe27 100644 --- a/test/bam_conversion.cxx +++ b/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("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) diff --git a/test/parallel_conversion_test.cxx b/test/parallel_conversion_test.cxx index 8f786df..e8b6ce5 100644 --- a/test/parallel_conversion_test.cxx +++ b/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); } diff --git a/test/ramcoretests.cxx b/test/ramcoretests.cxx index 97935e2..f4d6252 100644 --- a/test/ramcoretests.cxx +++ b/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 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) { diff --git a/test/ramdump_test.cxx b/test/ramdump_test.cxx index c62a54d..f9c1a78 100644 --- a/test/ramdump_test.cxx +++ b/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()); diff --git a/tools/ramdump.cxx b/tools/ramdump.cxx index a0b7907..f2ca46b 100644 --- a/tools/ramdump.cxx +++ b/tools/ramdump.cxx @@ -77,8 +77,8 @@ class RecordReader { ROOT::RNTupleView m_mapq; ROOT::RNTupleView> m_cigar; ROOT::RNTupleView m_refnext; - ROOT::RNTupleView m_pnext; - ROOT::RNTupleView m_tlen; + ROOT::RNTupleView m_pnextDelta; + ROOT::RNTupleView m_tlenDelta; ROOT::RNTupleView m_seq; ROOT::RNTupleView m_qual; ROOT::RNTupleView 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();