Feature head mapping - #45
Conversation
Introduce new function to map features to head instead of protein clusters and change the approach of reading from parquet_duckdb to only parquets.
Some reshuffling but mostly just adding support for manifest and modifying the way parameters are handled to account for manifest detection. Should be able to dynamically add databases through this, but accounting for different column names in those new DBs will take additional effort.
epbrenner
left a comment
There was a problem hiding this comment.
Works for me! But I also added the manifest stuff to this, so somebody else needs to approve before merging.
The merge-base changed after approval.
|
The integration of manifest JSON in the This PR can be merged. |
epbrenner
left a comment
There was a problem hiding this comment.
I still think this PR is cool. I still also think it could use one more review.
The merge-base changed after approval.
…database name column for HMMER generated matrices.
jananiravi
left a comment
There was a problem hiding this comment.
Few quick comments for @AbhirupaGhosh @eboyer221 @epbrenner
| #' ) | ||
| #' } | ||
| #' | ||
| #' @export |
There was a problem hiding this comment.
renegenerate namespace devtools::document()
| #' #' @details | ||
| #' The function performs the following steps: | ||
| #' \enumerate{ | ||
| #' \item Creates a protein-gene mapping from | ||
| #' \code{genome_gene_protein.parquet}. | ||
| #' \item Constructs a unique protein-gene dyad identifier of the form | ||
| #' \code{"protein|gene"}. | ||
| #' \item Optionally loads structural and HMMER annotations. | ||
| #' \item Generates a network edge list linking each dyad to its associated | ||
| #' features. | ||
| #' \item Writes the resulting edge list to a compressed Parquet file. | ||
| #' } |
There was a problem hiding this comment.
check rendering in roxygen2 documentation block. fix otherwise
| protein_ids AS protein, | ||
| REPLACE(Gene, '~', '.') AS gene | ||
| FROM read_parquet('%s') | ||
| WHERE protein_ids IS NOT NULL |
There was a problem hiding this comment.
Check: missing empty-string filter, inconsistent with data_processing.R's value != "" convention; can produce bogus "|Gene" dyads.
| DBI::dbReadTable(con, count_table) |> | ||
| tidyr::pivot_longer( | ||
| -genome_id, | ||
| names_to = "annotation", | ||
| values_to = "value" | ||
| ) |> | ||
| dplyr::rename(!!database := annotation) |> | ||
| dplyr::filter(!is.na(value) & value != "") |> | ||
| dplyr::mutate(value = as.integer(value)) |> | ||
| writeCompressedParquet(count_parquet) |
There was a problem hiding this comment.
any new features added here (this PR) that are not part of the PR description -- please go back and edit before merging. thanks!
| if ("struct" %in% addtnl_feature_scales) { | ||
|
|
||
| struct_path <- file.path( | ||
| parquet_dir, | ||
| "struct.parquet" | ||
| ) | ||
|
|
||
| if (!file.exists(struct_path)) { | ||
|
|
||
| message( | ||
| "Skipping struct: no parquet found. Generate struct parquet first." | ||
| ) | ||
|
|
There was a problem hiding this comment.
Duplicates the skip-pattern already generalized in create_feature_view() - minor.
| fixed = TRUE | ||
| ) | ||
| } | ||
|
|
There was a problem hiding this comment.
for each short block, worth adding a brief comment on what the succeeding block does. 'brief' would suffice (not just here -- this is true across scripts.
| @@ -0,0 +1,578 @@ | |||
| #' Build a protein-gene dyad feature network using DuckDB | |||
| #' | |||
| #' Constructs a bipartite network linking protein-gene dyads to biological | |||
There was a problem hiding this comment.
you are talking about many kinds of nodes but are calling it bipartite. 🤔
| #' protein|gene --> pfam:PFXXXXX | ||
| #' protein|gene --> cog:COGXXXX | ||
| #' protein|gene --> amr:GENE_NAME | ||
| #' protein|gene --> defensecas:FEATURE |
There was a problem hiding this comment.
make sure the generic FEATURE call doesn't lead to any downstream ambiguity. PD_FEATURE (for phage defense)?
|
@AbhirupaGhosh, please update the current PR description with accurate details on what the PR adds/fixes and how to run/check. |
Co-authored-by: Janani Ravi <janani.ravi@cuanschutz.edu>
| #' } | ||
| #' | ||
| #' @export | ||
| buildDyadFeatureMap <- function( |
There was a problem hiding this comment.
@eboyer221 check if new functions have tests -- if not, add that to your list of bioc to-dos.
| #' @export | ||
| buildDyadFeatureMap <- function( | ||
| duckdb_path, | ||
| addtnl_feature_scales = NULL, |
There was a problem hiding this comment.
is this used anywhere else? all other parameters seem to be fully spelled out
| # Generic HMMER feature view | ||
| # ========================= | ||
|
|
||
| create_feature_view <- function( |
There was a problem hiding this comment.
internal function starting with . or external with @export? add min. documentation?
@AbhirupaGhosh
| #' | ||
| #' @export | ||
| buildDyadFeatureMap <- function( | ||
| duckdb_path, |
There was a problem hiding this comment.
check for consistency
- when is this duckdb_parquet_path vs. duckdb_path?
| #' ) | ||
| #' } | ||
| #' | ||
| #' @export |
There was a problem hiding this comment.
Missing @import DBI duckdb here — related functions like buildClusterFeatureMap() in R/feature_to_cluster.R declares it for equivalent DBI/duckdb usage. Harmless, but inconsistent.
| edge_queries, | ||
| " | ||
| SELECT DISTINCT | ||
| pgd.dyad AS source, |
There was a problem hiding this comment.
what's pgd.dyad?
protein-gene-dyad.dyad?
…ation - Updated documentation references in `.defenseHMMER`, `.docker_path`, `.id_checker`, `.isValidHmmFile`, `.parseHMMEROutput`, `.parse_hmmer_profiles`, `.pick_shell`, `.prepareHmmerDatabases`, `.runHMMER`, and `.runHmmerJob` to point to the correct source files. - Removed obsolete documentation files for `.getDfIPRColNames` and `.getDfIPRColTypes`. - Added new internal functions for manifest management: `.manifest_event`, `.manifest_file_info`, `.manifest_finish`, `.manifest_git_info`, `.manifest_package_versions`, `.manifest_resume`, `.manifest_run_id`, `.manifest_stage`, and `.manifest_start`. - Introduced `.preserve_export_id_text` helper function to maintain trailing zeroes in exported tables. - Updated `runDataProcessing` and `runPanaroo2Duckdb` documentation to reflect changes in the processing pipeline, including the integration of HMMER. - Enhanced `exportProcessedData` and `exportTables` functions to include new parameters for better control over output formats and table exports. - Removed `.process_chunk` documentation as the function has been deprecated. - Added new function `buildDyadFeatureMap` for constructing a protein-gene dyad feature network.
# Conflicts: # R/helpers.R
…ature_scales, add @import, empty-string filter on protein_gene view, rework comments, add tests
|
Thanks @jananiravi, the requested changes have been addressed in the latest push. Docs / roxygen
Comments
Correctness
Tests
Merge conflict with main (R/helpers.R roxygen) resolved. PR description still needs to be updated. |
There was a problem hiding this comment.
@AbhirupaGhosh
This looks good to me! I did push a commit to address some of @jananiravi's review comments (roxygen @details rendering, addtnl_feature_scales → additional_feature_scales, @import DBI duckdb, empty-string filter on the protein_gene view, .create_feature_view rename, reworded comments, + a test file for the input-validation paths). Also merged main to clear the conflict that was showing up.
Tested: ran buildDyadFeatureMap() on the Staphylococcus_epidermidis sample Parquets (COG + AMRFinder + struct) and it reproduces the committed dyad_feature.parquet exactly - 483,686 edges, zero row-level diff. load_all()/document() clean, full test suite passes.
Still open: I think the only thing left is just updating the PR description, and a couple of the questions that are calls for you (the PD_/DefenseCas target-prefix naming; whether to rename feature_to_cluster.R's duckdb_parquet_path for consistency).
|
@eboyer221 One more thing to note is duckdb_path vs duckdb_parquet_path indicates two different database. One contains tables and other contains parquet views. |
jananiravi
left a comment
There was a problem hiding this comment.
I'll let Abhirupa go through the most recent changes, too. Looks OK otherwise.
Thanks for pointing that out, @AbhirupaGhosh, I hadn't taken note that duckdb_path (tables) and duckdb_parquet_path (views) are genuinely different DBs. I just look looked at PR #47. It settles the naming question because epbrenner's branch already adds buildDyadFeatureMap(duckdb_path = duckdb_path, ...) into runDataProcessing() (at the end, after cleanData()). So let's keep duckdb_path in #45; renaming to duckdb_parquet_path would just break #47's call site. Functionally it doesn't matter anyway because buildDyadFeatureMap() never opens the DB, it only uses dirname(duckdb_path) to find the sibling parquet files + manifest. Two things for merge coordination:
|
Description
Introduce feature-dyad mapping in the format of two columns, where dyad is protein|gene combination.
What kind of change(s) are included?
Checklist
Please ensure that all boxes are checked before indicating that this pull request is ready for review.