Refactor/output direct to h5ad#250
Conversation
There was a problem hiding this comment.
Code Review
This pull request replaces the external Python-based pycyto convert step with native Rust writing of .h5ad files directly from cyto-ibu-count using the anndata and anndata-hdf5 crates, eliminating the need to stage through MTX files or spawn a Python subprocess. To prevent spurious file-locking issues under concurrent Rayon execution, subprocesses for cell-filter and geomux are now invoked with HDF5_USE_FILE_LOCKING=FALSE. The review feedback focuses on improving the robustness and performance of the new write_counts_h5ad function. Specifically, it is recommended to dynamically determine the number of observations (n_obs) from the populated obs_names vector to prevent potential mismatches, and to pass the features vector by value (Vec) instead of a slice to avoid cloning the entire list of feature strings.
| let n_var = features.len(); | ||
| let n_obs = counts.get_num_barcodes(); | ||
| let nnz = counts.get_nnz(); | ||
|
|
||
| let mut obs_names = Vec::with_capacity(n_obs); | ||
| let mut bc_idx_map = HashMap::with_capacity(n_obs); |
There was a problem hiding this comment.
Using counts.get_num_barcodes() as n_obs before the loop can lead to a mismatch or failure if some barcodes are filtered out or if there is any discrepancy with the actual number of unique barcodes decoded in the loop. Since obs_names is populated dynamically during the loop, it is much safer and more robust to define n_obs as obs_names.len() after the loop. We can still use counts.get_num_barcodes() to pre-allocate the capacities of obs_names and bc_idx_map.
| let n_var = features.len(); | |
| let n_obs = counts.get_num_barcodes(); | |
| let nnz = counts.get_nnz(); | |
| let mut obs_names = Vec::with_capacity(n_obs); | |
| let mut bc_idx_map = HashMap::with_capacity(n_obs); | |
| let n_var = features.len(); | |
| let n_barcodes = counts.get_num_barcodes(); | |
| let nnz = counts.get_nnz(); | |
| let mut obs_names = Vec::with_capacity(n_barcodes); | |
| let mut bc_idx_map = HashMap::with_capacity(n_barcodes); |
| let coo = CooMatrix::try_from_triplets(n_obs, n_var, row_idx, col_idx, vals) | ||
| .map_err(|e| anyhow::anyhow!("Unable to build sparse count matrix: {e}"))?; |
There was a problem hiding this comment.
Use obs_names.len() as n_obs to ensure the number of rows in the sparse matrix matches the length of obs_names exactly, preventing any potential out-of-bounds errors or mismatches.
| let coo = CooMatrix::try_from_triplets(n_obs, n_var, row_idx, col_idx, vals) | |
| .map_err(|e| anyhow::anyhow!("Unable to build sparse count matrix: {e}"))?; | |
| let n_obs = obs_names.len(); | |
| let coo = CooMatrix::try_from_triplets(n_obs, n_var, row_idx, col_idx, vals) | |
| .map_err(|e| anyhow::anyhow!("Unable to build sparse count matrix: {e}"))?; |
| pub fn write_counts_h5ad<P: AsRef<Path>>( | ||
| path: P, | ||
| counts: &BarcodeIndexCounts, | ||
| features: &[String], | ||
| header: Header, | ||
| suffix: Option<&str>, | ||
| ) -> Result<()> { |
There was a problem hiding this comment.
To avoid cloning the entire vector of feature strings (which can contain tens of thousands of elements), we can pass features by value as Vec<String> instead of a slice &[String]. Since write_counts_h5ad is the final step in the counting process and consumes the features, passing by value is highly efficient and avoids unnecessary heap allocations.
| pub fn write_counts_h5ad<P: AsRef<Path>>( | |
| path: P, | |
| counts: &BarcodeIndexCounts, | |
| features: &[String], | |
| header: Header, | |
| suffix: Option<&str>, | |
| ) -> Result<()> { | |
| pub fn write_counts_h5ad<P: AsRef<Path>>( | |
| path: P, | |
| counts: &BarcodeIndexCounts, | |
| features: Vec<String>, | |
| header: Header, | |
| suffix: Option<&str>, | |
| ) -> Result<()> { |
| let adata = AnnData::<H5>::new(path.as_ref())?; | ||
| adata.set_x(csr)?; | ||
| adata.set_obs_names(DataFrameIndex::from(obs_names))?; | ||
| adata.set_var_names(DataFrameIndex::from(features.to_vec()))?; |
There was a problem hiding this comment.
Since features is now owned by the function, we can pass it directly to DataFrameIndex::from without calling .to_vec(), completely eliminating the need to clone the feature strings.
| adata.set_var_names(DataFrameIndex::from(features.to_vec()))?; | |
| adata.set_var_names(DataFrameIndex::from(features))?; |
| features | ||
| .expect("Must provide a feature file to write h5ad") | ||
| .as_slice(), |
No description provided.