Fix loom matrix orientation

Author: pkharchenko

CHANGELOG.md

@@ -56,6 +56,8 @@ All notable changes to the C++ line of Baysor are documented here.
 
 ### Fixed
 
+- Fixed Loom matrix orientation so `/matrix` rows match `/row_attrs` genes and
+  columns match `/col_attrs` cells.
 - Reduced clustering memory by removing the per-thread dense gene-by-gene correlation matrix.
 - Reduced segmentation memory by switching persistent component gene counts away from dense `double` storage.
 - Reduced NCV memory by removing global all-molecule neighborhood materialization and the full all-molecule high-dimensional NCV matrix from the `run` and `preview` paths.

docs/output_files.md

@@ -149,7 +149,7 @@ Structure:
 HDF5 layout written by the current C++ implementation:
 - `/matrix`
   - dataset type: `float32`
-  - shape: `(C, G)`
+  - shape: `(G, C)`
 - `/attrs/LOOM_SPEC_VERSION`
   - variable-length UTF-8 string
   - value: `"3.0.0"`
@@ -171,7 +171,7 @@ HDF5 layout written by the current C++ implementation:
 
 Notes:
 - matrix values are counts
-- the matrix is written in a row-oriented block path for performance
+- rows are genes and columns are cells, matching Loom row/column attributes
 
 ### `segmentation_counts.tsv`
 

include/baysor/reporting/output.h

@@ -73,6 +73,7 @@ using LoomColAttrs = std::map<std::string,
 
 /// Save count matrix to Loom (HDF5) format.
 /// matrix  : n_cells × n_genes  (rows = cells/NCVs, cols = genes)
+/// on disk : /matrix is n_genes × n_cells to match Loom row/column attributes
 /// col_attrs: optional extra per-cell attributes (e.g. ncv_color, confidence)
 void save_matrix_to_loom(const Eigen::SparseMatrix<float>& matrix,
                          const std::vector<std::string>& gene_names,

src/reporting/output.cpp

@@ -241,7 +241,7 @@ static nlohmann::json polygons_to_geojson_json(
 // Loom spec: https://linnarssonlab.org/loompy/format/index.html
 //
 // Layout produced (mirrors Julia's save_matrix_to_loom):
-//   /matrix            float32  n_cells × n_genes  (chunked, shuffle+deflate)
+//   /matrix            float32  n_genes × n_cells  (chunked, shuffle+deflate)
 //   /row_attrs/Name    variable-length UTF-8 strings  [n_genes]
 //   /col_attrs/Name    variable-length UTF-8 strings  [n_cells]
 //   /col_attrs/CellID  float64  [n_cells]   (1-based indices)
@@ -312,18 +312,26 @@ void save_matrix_to_loom(
         throw std::runtime_error("save_matrix_to_loom: gene_names length mismatch");
     if (static_cast<int>(cell_names.size()) != n_cells)
         throw std::runtime_error("save_matrix_to_loom: cell_names length mismatch");
+    for (const auto& [key, val] : col_attrs) {
+        size_t attr_len = std::holds_alternative<std::vector<std::string>>(val)
+            ? std::get<std::vector<std::string>>(val).size()
+            : std::get<std::vector<double>>(val).size();
+        if (attr_len != static_cast<size_t>(n_cells)) {
+            throw std::runtime_error("save_matrix_to_loom: col_attrs '" + key + "' length mismatch");
+        }
+    }
 
     // Open file
     hid_t fid = H5Fcreate(path.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
     check(fid, "H5Fcreate");
 
     // ---- /matrix (float32, chunked, shuffle + deflate-3) ----
     {
-        hsize_t dims[2]  = { static_cast<hsize_t>(n_cells),
-                             static_cast<hsize_t>(n_genes) };
-        constexpr int row_block = 128;
-        hsize_t chunk[2] = { static_cast<hsize_t>(std::min(n_cells, row_block)),
-                             static_cast<hsize_t>(std::min(n_genes, 256)) };
+        hsize_t dims[2]  = { static_cast<hsize_t>(n_genes),
+                             static_cast<hsize_t>(n_cells) };
+        constexpr int cell_block = 128;
+        hsize_t chunk[2] = { static_cast<hsize_t>(std::min(n_genes, 256)),
+                             static_cast<hsize_t>(std::min(n_cells, cell_block)) };
 
         hid_t space  = H5Screate_simple(2, dims, nullptr);
         hid_t plist  = H5Pcreate(H5P_DATASET_CREATE);
@@ -335,27 +343,27 @@ void save_matrix_to_loom(
                                H5P_DEFAULT, plist, H5P_DEFAULT);
         check(ds, "create /matrix");
 
-        // Write in row blocks to reduce HDF5 call overhead without
-        // materializing the full dense matrix.
-        std::vector<float> block_buf(static_cast<size_t>(row_block) * static_cast<size_t>(n_genes), 0.0f);
+        // Write Loom rows as genes. The input matrix stays cells x genes to
+        // match the other count-matrix writers.
+        std::vector<float> block_buf(static_cast<size_t>(n_genes) * static_cast<size_t>(cell_block), 0.0f);
         hsize_t block_offset[2] = {0, 0};
 
-        for (int block_start = 0; block_start < n_cells; block_start += row_block) {
-            int n_block_rows = std::min(row_block, n_cells - block_start);
-            std::fill(block_buf.begin(), block_buf.begin() + static_cast<size_t>(n_block_rows) * static_cast<size_t>(n_genes), 0.0f);
-            for (int local_row = 0; local_row < n_block_rows; ++local_row) {
-                int row = block_start + local_row;
-                float* row_ptr = block_buf.data() + static_cast<size_t>(local_row) * static_cast<size_t>(n_genes);
-                for (Eigen::SparseMatrix<float, Eigen::RowMajor>::InnerIterator it(matrix, row); it; ++it) {
-                    row_ptr[it.col()] = it.value();
+        for (int block_start = 0; block_start < n_cells; block_start += cell_block) {
+            int n_block_cols = std::min(cell_block, n_cells - block_start);
+            std::fill(block_buf.begin(), block_buf.begin() + static_cast<size_t>(n_genes) * static_cast<size_t>(n_block_cols), 0.0f);
+            for (int local_col = 0; local_col < n_block_cols; ++local_col) {
+                int cell = block_start + local_col;
+                for (Eigen::SparseMatrix<float, Eigen::RowMajor>::InnerIterator it(matrix, cell); it; ++it) {
+                    int gene = it.col();
+                    block_buf[static_cast<size_t>(gene) * static_cast<size_t>(n_block_cols) + static_cast<size_t>(local_col)] = it.value();
                 }
             }
 
             hsize_t block_dims[2] = {
-                static_cast<hsize_t>(n_block_rows),
-                static_cast<hsize_t>(n_genes)
+                static_cast<hsize_t>(n_genes),
+                static_cast<hsize_t>(n_block_cols)
             };
-            block_offset[0] = static_cast<hsize_t>(block_start);
+            block_offset[1] = static_cast<hsize_t>(block_start);
 
             hid_t file_space = H5Dget_space(ds);
             H5Sselect_hyperslab(file_space, H5S_SELECT_SET,

tests/test_main.cpp

@@ -2038,6 +2038,92 @@ TEST(Output, SaveMatrixTo10xH5WritesExpectedStructure) {
     H5Fclose(fid);
 }
 
+TEST(Output, SaveMatrixToLoomWritesGeneByCellMatrix) {
+    Eigen::SparseMatrix<float, Eigen::RowMajor> matrix(2, 3);  // cells x genes
+    std::vector<Eigen::Triplet<float>> trips = {
+        {0, 0, 2.0f},
+        {0, 2, 5.0f},
+        {1, 1, 7.0f}
+    };
+    matrix.setFromTriplets(trips.begin(), trips.end());
+    matrix.makeCompressed();
+
+    const std::vector<std::string> gene_names = {"G1", "G2", "G3"};
+    const std::vector<std::string> cell_names = {"cell_1", "cell_2"};
+    const std::string path = write_temp_csv("", ".loom");
+
+    baysor::save_matrix_to_loom(matrix, gene_names, cell_names, path);
+
+    hid_t fid = H5Fopen(path.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+    ASSERT_GE(fid, 0);
+
+    hid_t matrix_ds = H5Dopen2(fid, "/matrix", H5P_DEFAULT);
+    ASSERT_GE(matrix_ds, 0);
+    hid_t matrix_space = H5Dget_space(matrix_ds);
+    hsize_t matrix_dims[2] = {0, 0};
+    ASSERT_EQ(H5Sget_simple_extent_ndims(matrix_space), 2);
+    ASSERT_EQ(H5Sget_simple_extent_dims(matrix_space, matrix_dims, nullptr), 2);
+    EXPECT_EQ(matrix_dims[0], 3);
+    EXPECT_EQ(matrix_dims[1], 2);
+
+    std::vector<float> values(6, -1.0f);
+    ASSERT_GE(H5Dread(matrix_ds, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, values.data()), 0);
+    const std::vector<float> expected = {
+        2.0f, 0.0f,
+        0.0f, 7.0f,
+        5.0f, 0.0f
+    };
+    EXPECT_EQ(values, expected);
+    H5Sclose(matrix_space);
+    H5Dclose(matrix_ds);
+
+    hid_t row_name_ds = H5Dopen2(fid, "/row_attrs/Name", H5P_DEFAULT);
+    ASSERT_GE(row_name_ds, 0);
+    hid_t row_name_space = H5Dget_space(row_name_ds);
+    hsize_t row_dims[1] = {0};
+    ASSERT_EQ(H5Sget_simple_extent_ndims(row_name_space), 1);
+    ASSERT_EQ(H5Sget_simple_extent_dims(row_name_space, row_dims, nullptr), 1);
+    EXPECT_EQ(row_dims[0], 3);
+    H5Sclose(row_name_space);
+    H5Dclose(row_name_ds);
+
+    hid_t col_name_ds = H5Dopen2(fid, "/col_attrs/Name", H5P_DEFAULT);
+    ASSERT_GE(col_name_ds, 0);
+    hid_t col_name_space = H5Dget_space(col_name_ds);
+    hsize_t col_dims[1] = {0};
+    ASSERT_EQ(H5Sget_simple_extent_ndims(col_name_space), 1);
+    ASSERT_EQ(H5Sget_simple_extent_dims(col_name_space, col_dims, nullptr), 1);
+    EXPECT_EQ(col_dims[0], 2);
+    H5Sclose(col_name_space);
+    H5Dclose(col_name_ds);
+
+    hid_t cell_id_ds = H5Dopen2(fid, "/col_attrs/CellID", H5P_DEFAULT);
+    ASSERT_GE(cell_id_ds, 0);
+    hid_t cell_id_space = H5Dget_space(cell_id_ds);
+    hsize_t cell_id_dims[1] = {0};
+    ASSERT_EQ(H5Sget_simple_extent_ndims(cell_id_space), 1);
+    ASSERT_EQ(H5Sget_simple_extent_dims(cell_id_space, cell_id_dims, nullptr), 1);
+    EXPECT_EQ(cell_id_dims[0], 2);
+    H5Sclose(cell_id_space);
+    H5Dclose(cell_id_ds);
+
+    H5Fclose(fid);
+}
+
+TEST(Output, SaveMatrixToLoomRejectsMismatchedColAttrLength) {
+    Eigen::SparseMatrix<float> matrix(2, 3);  // cells x genes
+    const std::vector<std::string> gene_names = {"G1", "G2", "G3"};
+    const std::vector<std::string> cell_names = {"cell_1", "cell_2"};
+    const std::string path = write_temp_csv("", ".loom");
+
+    baysor::LoomColAttrs col_attrs;
+    col_attrs["confidence"] = std::vector<double>{0.5};
+
+    EXPECT_THROW(
+        baysor::save_matrix_to_loom(matrix, gene_names, cell_names, path, col_attrs),
+        std::runtime_error);
+}
+
 TEST(Xenium, ManifestHelpers) {
     EXPECT_TRUE(baysor::is_xenium_manifest_path("experiment.xenium"));
     EXPECT_FALSE(baysor::is_xenium_manifest_path("transcripts.parquet"));