Merge branch 'jdaw/fix-split-read-counting' into 'master'

Account for split reads during progress tracking

See merge request machine-learning/dorado!667

dorado/read_pipeline/HtsWriter.cpp

@@ -90,20 +90,33 @@ void HtsWriter::worker_thread() {
 
         auto aln = std::move(std::get<BamPtr>(message));
         write(aln.get());
-        std::string read_id = bam_get_qname(aln.get());
+
+        // For the purpose of estimating write count, we ignore duplex reads
         int64_t dx_tag = 0;
         auto tag_str = bam_aux_get(aln.get(), "dx");
         if (tag_str) {
             dx_tag = bam_aux2i(tag_str);
         }
 
-        // For the purpose of estimating write count, we ignore duplex reads
         bool ignore_read_id = dx_tag == 1;
 
         if (ignore_read_id) {
             // Read is a duplex read.
             m_duplex_reads_written++;
         } else {
+            std::string read_id;
+
+            // If read is a split read, use the parent read id
+            // to track write count since we don't know a priori
+            // how many split reads will be generated.
+            auto pid_tag = bam_aux_get(aln.get(), "pi");
+            if (pid_tag) {
+                read_id = std::string(bam_aux2Z(pid_tag));
+                m_split_reads_written++;
+            } else {
+                read_id = bam_get_qname(aln.get());
+            }
+
             m_processed_read_ids.insert(std::move(read_id));
         }
     }
@@ -150,6 +163,7 @@ stats::NamedStats HtsWriter::sample_stats() const {
     auto stats = stats::from_obj(m_work_queue);
     stats["unique_simplex_reads_written"] = m_processed_read_ids.size();
     stats["duplex_reads_written"] = m_duplex_reads_written.load();
+    stats["split_reads_written"] = m_split_reads_written.load();
     return stats;
 }
 

dorado/read_pipeline/HtsWriter.h

@@ -52,6 +52,7 @@ class HtsWriter : public MessageSink {
     size_t m_num_reads_expected;
     std::unordered_set<std::string> m_processed_read_ids;
     std::atomic<int> m_duplex_reads_written{0};
+    std::atomic<int> m_split_reads_written{0};
 };
 
 }  // namespace dorado

tests/BamWriterTest.cpp

@@ -2,6 +2,7 @@
 #include "read_pipeline/HtsReader.h"
 #include "read_pipeline/HtsWriter.h"
 #include "utils/bam_utils.h"
+#include "utils/stats.h"
 
 #include <catch2/catch.hpp>
 #include <htslib/sam.h>
@@ -35,9 +36,13 @@ class HtsWriterTestsFixture {
         writer_ref.set_and_write_header(reader.header);
 
         reader.read(*pipeline, 1000);
+        pipeline->terminate(DefaultFlushOptions());
+        stats = writer_ref.sample_stats();
         pipeline.reset();
     }
 
+    stats::NamedStats stats;
+
 private:
     fs::path m_in_sam;
     fs::path m_out_bam;
@@ -47,7 +52,7 @@ TEST_CASE_METHOD(HtsWriterTestsFixture, "HtsWriterTest: Write BAM", TEST_GROUP)
     int num_threads = GENERATE(1, 10);
     HtsWriter::OutputMode emit_fastq = GENERATE(
             HtsWriter::OutputMode::SAM, HtsWriter::OutputMode::BAM, HtsWriter::OutputMode::FASTQ);
-    REQUIRE_NOTHROW(generate_bam(emit_fastq, num_threads));
+    CHECK_NOTHROW(generate_bam(emit_fastq, num_threads));
 }
 
 TEST_CASE("HtsWriterTest: Output mode conversion", TEST_GROUP) {
@@ -56,3 +61,10 @@ TEST_CASE("HtsWriterTest: Output mode conversion", TEST_GROUP) {
     CHECK(HtsWriter::get_output_mode("fastq") == HtsWriter::OutputMode::FASTQ);
     CHECK_THROWS_WITH(HtsWriter::get_output_mode("blah"), "Unknown output mode: blah");
 }
+
+TEST_CASE_METHOD(HtsWriterTestsFixture, "HtsWriter: Count reads written", TEST_GROUP) {
+    CHECK_NOTHROW(generate_bam(HtsWriter::OutputMode::BAM, 1));
+
+    CHECK(stats.at("unique_simplex_reads_written") == 6);
+    CHECK(stats.at("split_reads_written") == 2);
+}

tests/data/bam_reader/small.sam

@@ -2,8 +2,8 @@
 @SQ	SN:NZ_CP062160.1	LN:5498557
 @SQ	SN:NZ_CP062161.1	LN:104844
 @PG	ID:minimap2	PN:minimap2	VN:2.24-r1155-dirty	CL:minimap2 -k 19 -w 19 -H -a --MD -o test2.sam /Users/joyjit.daw/references/ecoli-GCF_014825845.1_ASM1482584v1_genomic.fn /Users/joyjit.daw/Downloads/temp.fasta
-f25635d3-3b81-4ff2-8510-f28330590801	4	*	0	0	*	*	0	0	TGCTCAGATGGGGGAATGAAAACAATACTCGACAGGAACTACAACGCCACAGCGGCGATTGTATCAGAAAGTTGACACTGGAAGATAAACTGGTGGAGAAACCTGTCACGGGGAAGAGCAGAGCTTGACGCGCTTCGCTTGTCCCGGCCTCACGACCCGCGAAAGAGTTTCGCGAACCGTTCACCGTTTTTCAATATGATTCTGTCTGCTTTGGTCATTAAAGCAGACGATAGTGTCCAACTTCCCCGTATATAACGGTACGTGCTCATCGGTAAGAGCTGGCGGTATTGTGTTTATCAGGGGAAGTTTTCAGATACTCACGAAAATGCCTGGCAGTACCGCCGTCTACAATGCCATCGATAATATGTAGGTGTTTTCACCTAATTTTGCCATGATGATAATCCTTATATTACGGTTTGTGAATAGGTCTGCTTGTCAGGGTTTGTGTAACCTTTTAACCTTAGCAATGCGTAACGCGTGTACTTTGCGTTCAATGGTTACGTATTGCTGGGTTAGAGAGGACAAAGGTTTCAACTGCCCCCAGCACCTCTGGTTAATTGGATAGCGACACTTCCTGTATCATTCGGTGGCGGAAAATCAGTGGAATCGCCGAAATGGCTGGCGCGATTAACCATTTTGCGCCTGGCACCTGCTCTATTTTACCGATTTATCAACTGGCCTCGCCTGCGGCGAATTGTTAGTGGAATGATAGCTGAAGAGCTGCAAAGTTCAGG	*	rl:i:0
-95e6aa45-efe0-4ed9-b7b9-2f35164f9836	0	NZ_CP062160.1	1012347	60	1S43M1D8M1D11M1D5M1I10M1D2M5D2M1I37M1I102M2D14M2D25M3D90M2D7M1I35M3D8M1D73M1I28M1D9M1D2M1D9M3D5M3D28M1D21M4D62M1D6M1D5M1D52M1D2M1D39M1D45M1I30M2D119M1D2M2D20M2I76M1D12M1D4M1I176M1D22M780S	*	0	0	TGTAATCACGCCGTCTGCCAGTAGGATCACTTTGTTGGCGGTTGAGCGGGATAATGCCAATCGTTCTGCTCAATGCTGCCCGGCATCTCAGCAAAGGCCGCGACACGGGTGCTGAGCTGCCACGACATCCGAGGTGGAGTAGTCTGTTTCAATAGCGATATAAGGAATGTTGTGCTGCTGGCGCACATGACGTTTAATCGCCAGCGATTCCACCGCGTAGGTATCGCGCCTGCAAAATCATCAACTACGCCATCGACCTGATATTCAACCATCTGGCTGAGCATTTGCGGGCGCTGATCGTTCGGCGAAACACAGGAGCAGCCAATCGCCAGATATTTATCCGCCAGCGCGTCATACATCACCCCGTTTCCGCCACGCATTGCTCGGTGGCTTTCGCGGTGCAGTTTCATAACCGACAACCCAGCCGCCATTCTCTTCAATCGCGCGCACCACTTTCTCTGCTGCGCCGCCAATTGGAGCAGCCGGTGATTAAAATACGTGGGCGCGGTCCTGTCATGCCTTCTTCCTGCTGAACGCGGGCGGTCATTGCATCCAGTTCGTGATCAACGCCTCTTTATCGGGCAGGTTGCGCCGTAGACTACTTTCAGAATGTCGCTGCCGCTAAGCGCTGGAGGATTTAACGCACAAATGATAAAGTTGGCCAGTGCGCGACGTTCGCGGTTTTTCAGCGCAATGGCATCGCGCGATATCTTCGCTAATCTCGTGCCCAAAACGTTCTTCTACTGTTTTTGCAAGCGCAGCATCTCGGCTTTCCATAATGCACGCGAGGCGCTCGTCCTTAACGCTGTTTGGCAACTGCATCATGAACGGGCTTAAATTCCGCCATGTATTCATACATTTTCTTTTTGCCGTCACAGGTGGTTTCACCGACCACCAGATCAGAAAAACAGAAGTAGGGGCATTTATCGGTTTTGCCGAAGCGGCTGCTTTTAATCAGCGGGCAGAGAGGTTGCGCGGCAGATCTTTCTCGGCTTCTTCAATGGTTTCATCAGAGGTGGAACATAGCGAAACCACAACCACCCGGCTGCCATCGGATACTCTTGCGGCATAAAGGTGCAGTAGGTGCCAACCAGCGGAATGCCGCGCTCCTTGAGATCCATGACGGTGAGAAAGCCTTTCTGGCGAGCTTCAGAGAACTGATCGAAAATGGCGGGTAGCTCGGTAACAAGTGACATGATTTTCCTTCCCCGTACCACGGGAGATAATGAAAAGAGCCACATTATTACCATTCTTTTGTTGTACTTCGTTCGATTGCATTACTGGGTGGCAGCGGCGTTCGAACGGCACTGATATCAACTGAAACAGACCATGGCGATGAGCTGTTGGGCGAGCAGGAGGAGCTGGTGGGCAACAAGGTCGCGAAAAAGCTGACCACTGACGCTTCTGAATGGAAGTCTTGCTGGATGGCGAAAAGGTAAAGTGAAATGGTCGCTGGTGACGAGCATAATATGCGTAGCTGATGGCGATTGCGGCGGCTCGCCATGTTGGTGTAGCGGCGGATACTCACCAGGTTCGTTTTAACGCTCGTCGCCGTCTGGAGTTGCACCATCCGCCAGCACGGATCTGCTCCAGCTACCTGGCTGTGGTTCTGAATAATGGTGAGCATGTTTTCCCGGCGTTAACCGTAATGGTCATCGGCTCGCGCCATCTGTCGCGGTGACTGACTTGCGGAATATCGGCGTTAGCCATAAAAGCAGAGAACAGCGCAGCCGCATAACTGCAATACGACAGTTAATGGTCTTGACGAGACTATCCGTGTATGTTGTAACTAAACAGAGTTTATCGTATTTGCCCACGGGCGGGACGCAGTCGCTAAGTGGCAGAACAAATGCCACCAGCGCCGGGATAGATAAATAAACCAGACGGCGCGCTTCGCTGACGCTCGGTGCTTCCTGATAGCCAAATCCCTTCAGTGTGCCAAACAGCGAAATTGTGAGTTGAGCACCGCGCTCATATCGAAGGCGATTTCCTGAAGATGGTTCCACAATCCGGCTTCATGAAATGCGCGAATGGCA	*	NM:i:124	ms:i:1853	AS:i:1832	nn:i:0	tp:A:P	cm:i:24	s1:i:396	s2:i:0	de:f:0.0802	MD:Z:17G4A4C5G7A1^A1C1A1A2^C11^A0A1A5C4A1^T2^TTACA8A0T18A11A56A17T25^GG1A12^CA25^TTC3C19T1A64^GT1G5G8T17C7^CCC8^T70C20G2C2A3^G5A3^G0C1^G9^CCA5^GAC3T12C2G8^A0T20^AACC2A13A2C29C8G3^T2C3^G5^A3A48^A2^G0C12T20C2C1^C6T26C11A5T16G6^CA16C35G26C5G0T32^C2^TA3A41C32G15G1^A12^G28A24A70A5G49^G22	rl:i:0
+f25635d3-3b81-4ff2-8510-f28330590801	4	*	0	0	*	*	0	0	TGCTCAGATGGGGGAATGAAAACAATACTCGACAGGAACTACAACGCCACAGCGGCGATTGTATCAGAAAGTTGACACTGGAAGATAAACTGGTGGAGAAACCTGTCACGGGGAAGAGCAGAGCTTGACGCGCTTCGCTTGTCCCGGCCTCACGACCCGCGAAAGAGTTTCGCGAACCGTTCACCGTTTTTCAATATGATTCTGTCTGCTTTGGTCATTAAAGCAGACGATAGTGTCCAACTTCCCCGTATATAACGGTACGTGCTCATCGGTAAGAGCTGGCGGTATTGTGTTTATCAGGGGAAGTTTTCAGATACTCACGAAAATGCCTGGCAGTACCGCCGTCTACAATGCCATCGATAATATGTAGGTGTTTTCACCTAATTTTGCCATGATGATAATCCTTATATTACGGTTTGTGAATAGGTCTGCTTGTCAGGGTTTGTGTAACCTTTTAACCTTAGCAATGCGTAACGCGTGTACTTTGCGTTCAATGGTTACGTATTGCTGGGTTAGAGAGGACAAAGGTTTCAACTGCCCCCAGCACCTCTGGTTAATTGGATAGCGACACTTCCTGTATCATTCGGTGGCGGAAAATCAGTGGAATCGCCGAAATGGCTGGCGCGATTAACCATTTTGCGCCTGGCACCTGCTCTATTTTACCGATTTATCAACTGGCCTCGCCTGCGGCGAATTGTTAGTGGAATGATAGCTGAAGAGCTGCAAAGTTCAGG	*	rl:i:0	pi:Z:aaaaaaaaaa
+95e6aa45-efe0-4ed9-b7b9-2f35164f9836	0	NZ_CP062160.1	1012347	60	1S43M1D8M1D11M1D5M1I10M1D2M5D2M1I37M1I102M2D14M2D25M3D90M2D7M1I35M3D8M1D73M1I28M1D9M1D2M1D9M3D5M3D28M1D21M4D62M1D6M1D5M1D52M1D2M1D39M1D45M1I30M2D119M1D2M2D20M2I76M1D12M1D4M1I176M1D22M780S	*	0	0	TGTAATCACGCCGTCTGCCAGTAGGATCACTTTGTTGGCGGTTGAGCGGGATAATGCCAATCGTTCTGCTCAATGCTGCCCGGCATCTCAGCAAAGGCCGCGACACGGGTGCTGAGCTGCCACGACATCCGAGGTGGAGTAGTCTGTTTCAATAGCGATATAAGGAATGTTGTGCTGCTGGCGCACATGACGTTTAATCGCCAGCGATTCCACCGCGTAGGTATCGCGCCTGCAAAATCATCAACTACGCCATCGACCTGATATTCAACCATCTGGCTGAGCATTTGCGGGCGCTGATCGTTCGGCGAAACACAGGAGCAGCCAATCGCCAGATATTTATCCGCCAGCGCGTCATACATCACCCCGTTTCCGCCACGCATTGCTCGGTGGCTTTCGCGGTGCAGTTTCATAACCGACAACCCAGCCGCCATTCTCTTCAATCGCGCGCACCACTTTCTCTGCTGCGCCGCCAATTGGAGCAGCCGGTGATTAAAATACGTGGGCGCGGTCCTGTCATGCCTTCTTCCTGCTGAACGCGGGCGGTCATTGCATCCAGTTCGTGATCAACGCCTCTTTATCGGGCAGGTTGCGCCGTAGACTACTTTCAGAATGTCGCTGCCGCTAAGCGCTGGAGGATTTAACGCACAAATGATAAAGTTGGCCAGTGCGCGACGTTCGCGGTTTTTCAGCGCAATGGCATCGCGCGATATCTTCGCTAATCTCGTGCCCAAAACGTTCTTCTACTGTTTTTGCAAGCGCAGCATCTCGGCTTTCCATAATGCACGCGAGGCGCTCGTCCTTAACGCTGTTTGGCAACTGCATCATGAACGGGCTTAAATTCCGCCATGTATTCATACATTTTCTTTTTGCCGTCACAGGTGGTTTCACCGACCACCAGATCAGAAAAACAGAAGTAGGGGCATTTATCGGTTTTGCCGAAGCGGCTGCTTTTAATCAGCGGGCAGAGAGGTTGCGCGGCAGATCTTTCTCGGCTTCTTCAATGGTTTCATCAGAGGTGGAACATAGCGAAACCACAACCACCCGGCTGCCATCGGATACTCTTGCGGCATAAAGGTGCAGTAGGTGCCAACCAGCGGAATGCCGCGCTCCTTGAGATCCATGACGGTGAGAAAGCCTTTCTGGCGAGCTTCAGAGAACTGATCGAAAATGGCGGGTAGCTCGGTAACAAGTGACATGATTTTCCTTCCCCGTACCACGGGAGATAATGAAAAGAGCCACATTATTACCATTCTTTTGTTGTACTTCGTTCGATTGCATTACTGGGTGGCAGCGGCGTTCGAACGGCACTGATATCAACTGAAACAGACCATGGCGATGAGCTGTTGGGCGAGCAGGAGGAGCTGGTGGGCAACAAGGTCGCGAAAAAGCTGACCACTGACGCTTCTGAATGGAAGTCTTGCTGGATGGCGAAAAGGTAAAGTGAAATGGTCGCTGGTGACGAGCATAATATGCGTAGCTGATGGCGATTGCGGCGGCTCGCCATGTTGGTGTAGCGGCGGATACTCACCAGGTTCGTTTTAACGCTCGTCGCCGTCTGGAGTTGCACCATCCGCCAGCACGGATCTGCTCCAGCTACCTGGCTGTGGTTCTGAATAATGGTGAGCATGTTTTCCCGGCGTTAACCGTAATGGTCATCGGCTCGCGCCATCTGTCGCGGTGACTGACTTGCGGAATATCGGCGTTAGCCATAAAAGCAGAGAACAGCGCAGCCGCATAACTGCAATACGACAGTTAATGGTCTTGACGAGACTATCCGTGTATGTTGTAACTAAACAGAGTTTATCGTATTTGCCCACGGGCGGGACGCAGTCGCTAAGTGGCAGAACAAATGCCACCAGCGCCGGGATAGATAAATAAACCAGACGGCGCGCTTCGCTGACGCTCGGTGCTTCCTGATAGCCAAATCCCTTCAGTGTGCCAAACAGCGAAATTGTGAGTTGAGCACCGCGCTCATATCGAAGGCGATTTCCTGAAGATGGTTCCACAATCCGGCTTCATGAAATGCGCGAATGGCA	*	NM:i:124	ms:i:1853	AS:i:1832	nn:i:0	tp:A:P	cm:i:24	s1:i:396	s2:i:0	de:f:0.0802	MD:Z:17G4A4C5G7A1^A1C1A1A2^C11^A0A1A5C4A1^T2^TTACA8A0T18A11A56A17T25^GG1A12^CA25^TTC3C19T1A64^GT1G5G8T17C7^CCC8^T70C20G2C2A3^G5A3^G0C1^G9^CCA5^GAC3T12C2G8^A0T20^AACC2A13A2C29C8G3^T2C3^G5^A3A48^A2^G0C12T20C2C1^C6T26C11A5T16G6^CA16C35G26C5G0T32^C2^TA3A41C32G15G1^A12^G28A24A70A5G49^G22	rl:i:0	pi:Z:aaaaaaaaaa
 13d1b2d7-b6a8-4c58-a709-49f0f341c26f	4	*	0	0	*	*	0	0	CTACCCAGATTTTTGTGCAAAACCGTGACGGTACTGATCGGCATGGGGAAATTCTTCGCTATCCTGGTCACCGGTATCGACCTTTCGGTTAGCGCGATTCTGGCGCTTTCCGGTATGGTGACCGCTAAACTGATGTTGGCAGGTATTGATCCGTTTCTCGCGGCGCTAATTGGCGGTGTACTGGTTGGCGGCGCACTGGGGCGATCAACGGCTGCCTGGTCAACTGGACGGGACTACACCCGTTCATTATCACCTTAGCACCAACGCTATTTTCCGTGGGATCACGCTGGTGATCTCCGATGCCAACTCGGTAGCCGACTTCTCATTTGACCGTGAACTTCTTTTGCCGCCAGCGTAATTGGGATACCTGTCCCCGTTATCTTCTCGCTAATTGCTTCGCGCTCATCCTTTGGTTTCTGACAACGCGTATGCGCTCGGGCGCAACATCTACGCACTGGGCGGCAGCAAAAACTCGGCGTTATTCCGGGATTGACGTGAAATAGGTGCTGAAGCGTTGAGTATCCTCTTGGCGCCTTAGCAGTACGTAGCCTTGTGTTGCCGCGTTTCAATTTACGTATTTGCTAAAAGGTTAAACAGGCGACTACAAACAGAATCGACGCTTTACCCGCAAACAATGCCCGTTGAAATTTATCCGCCGTTCATAGCCTCACCTCCGCAAATAACGGATGGCGTAGTTTTACACTGAGAAATGAAAGGATTTGAAAAAACCGCAAAGCGGGCGAAACGATATATACAGAAAGAAAGGAAAGCACTCTATCCAACAAACACACCCACAGTTGATCGGAATAAAAAGCAGA	*	rl:i:0
 489f7224-7d18-41d4-9129-d28c3e6aa780	16	NZ_CP062160.1	1629814	60	2S9M1I95M1D3M1D53M1I17M1D75M1I46M1I6M1I23M1I22M168S	*	0	0	CTTATCATGTCTTTTTGGGTGAGTATTCATCCATAATGCGTCCTCTTCTTAGCGGTTGAACTAACGGACACCTTTCGGGATGGAAAAAACTTACTGACCTGGACTTGCCTTCGTTTGTTAGCTTAACTATAAGCCACTCTTTGCAGGTTTTCATCGCATTATTAATGAAAAATTGCAATTAGAAGGAAGCGGGACAGGACTCACTTTTTCAGCCTAGCACCAACCCGCAGCAGGTAAAAGCAGTTTGCCCGAACTGTTTCATAACATTTCAGAATTATCGCCAGAAGAAATTGGTGACACTAAGTACTGGTGATTGCTCTGAATATGACGTAATAAATCGAGGAATGAATAAAGAATATGTAACGTCCGGCTTCATCCGGATTTTCTGAGGCTACCGTGTTCACCACAACCGTGGTGCTGTTACGTTTTGTTTTGAGCTGAATGTGCAGTTTCAGTCGGTTTCCTGCACCGTCTTTCAGCACGCCTGAAATCTGTACTGCCATATTCACTCCACAAATAAAAAAG	*	NM:i:22	ms:i:566	AS:i:566	nn:i:0	tp:A:P	cm:i:6	s1:i:121	s2:i:0	de:f:0.0615	MD:Z:33C6C63^C3^C54C8A6^G17G12T0A35T22G12T20C0A26C19	rl:i:0
 d7500028-dfcc-4404-b636-13edae804c55	16	NZ_CP062160.1	3164467	1	174M2D20M329S	*	0	0	CGCCGCCGGTGCGGCAATCCGGAACGATACCGATGCCGGATCGCCCCGCTGCCCCCACGCATTTACTGCCCGGACTGTCAGCGTGTACCGCCCCAGCGCCAGTTGCGTGAAGCGGTATGTGGTTTCCGTCGTCCGGGCCGTGCTGACCAGCCGCTCACTGCCGTCATCCGCTGCCGGTCAGGCGAAGCATAAAGGTGCTGTCGATTCCGTTTGTAGTCGTCTGTTTAACCTTAGCAATACGTAACTGAACGAAGTACAACACATTTTTTTTTTTTTTTTTTCAGTTACGTATTGCCAAGGTTAAGAGAGGTCAGCGTTTCAACGCTCAGCACCTATCGCTCCGGAAGTGGACGTTGACGACGAGCCAGAAGAAGAATAATTTCATCGCTTTCATGCCAGGGAAAAGGGAGCCATCTCCTCTTTGAATTGAAAAGTCCAGGCTGTAAAGTCTGGGCTTTTGTCGTATTAGGCGCGGTGTTTGGCTGTGCCTCGTAAAAAATGGCTGGCTATACACAAGGAATAG	*	NM:i:12	ms:i:321	AS:i:320	nn:i:0	tp:A:P	cm:i:4	s1:i:70	s2:i:77	de:f:0.0564	SA:Z:NZ_CP062160.1,4112332,+,2S187M334S,19,20;	MD:Z:3T29C12T19C15C4G0T0T4G76G2^CA20	rl:i:0