feat: add circle annotations #24
Conversation
…nload with new wrapper)
…ownload via wrappers
…directory with one table per circle in the following commits)
…ability of GFF3 format for regulatory annotations
|
Important Review skippedDraft detected. Please check the settings in the CodeRabbit UI or the You can disable this status message by setting the WalkthroughThe changes introduce new annotation processing capabilities to the workflow, including rules and scripts for retrieving, harmonizing, and annotating genomic and regulatory features. New Conda environments are defined for these steps. Configuration files are updated for Ensembl release version and chromosome handling. Minor code improvements and utility functions are also added. Changes
Sequence Diagram(s)sequenceDiagram
participant User
participant Snakemake
participant R_Scripts
participant Conda_Env
User->>Snakemake: Trigger workflow
Snakemake->>Conda_Env: Activate biomart/rtracklayer/annotatr envs as needed
Snakemake->>R_Scripts: Run create_transcripts_to_genes_mappings.R
Snakemake->>R_Scripts: Run create_annotation_gff3.R
Snakemake->>R_Scripts: Run annotate_cleaned_circles.R
R_Scripts->>Snakemake: Output harmonized annotation files and per-circle TSVs
Poem
Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. 🪧 TipsChatThere are 3 ways to chat with CodeRabbit:
SupportNeed help? Create a ticket on our support page for assistance with any issues or questions. Note: Be mindful of the bot's finite context window. It's strongly recommended to break down tasks such as reading entire modules into smaller chunks. For a focused discussion, use review comments to chat about specific files and their changes, instead of using the PR comments. CodeRabbit Commands (Invoked using PR comments)
Other keywords and placeholders
CodeRabbit Configuration File (
|
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 9
🧹 Nitpick comments (5)
workflow/envs/rtracklayer.yaml (1)
1-7: Normalize YAML formatting and add EOF newline.
Ensure consistent spacing around the=in version specs (e.g., changebioconductor-rtracklayer =1.62to= 1.62) and add a final newline to satisfy linters.🧰 Tools
🪛 YAMLlint (1.37.1)
[error] 7-7: no new line character at the end of file
(new-line-at-end-of-file)
workflow/envs/annotatr.yaml (1)
1-7: Normalize YAML formatting and add EOF newline.
Ensure consistent spacing around the=in version specs (e.g., changebioconductor-annotatr =1.28to= 1.28) and append a trailing newline for YAML lint compliance.🧰 Tools
🪛 YAMLlint (1.37.1)
[error] 7-7: no new line character at the end of file
(new-line-at-end-of-file)
workflow/envs/biomart.yaml (1)
1-7: Normalize YAML formatting and add EOF newline.
Standardize spacing around=in version pins (e.g.,bioconductor-biomart = 2.58) and include a final newline to avoid linter errors.🧰 Tools
🪛 YAMLlint (1.37.1)
[error] 7-7: no new line character at the end of file
(new-line-at-end-of-file)
.test/config/config.yaml (1)
18-18: Improve comment formatting.
Add a space after#for consistency (change#chromosome: 7to# chromosome: 7).workflow/rules/common.smk (1)
86-90:get_bioc_species_name()produces the desired abbreviation but is fragile– Fails hard if
config["ref"]["species"]is missing the expected single underscore (e.g."drosophila_melanogaster","pan_troglodytes"are fine, but"yeast"is not).
– Upper-/lower-case is preserved; Bioconductor dataset names expect e.g.Hsapiens,Mmusculus, etc. Ifconfig["ref"]["species"]starts all-lower-case ("homo_sapiens"), the function will returnhsapiens, which will not match package names.Consider adding minimal validation & normalisation:
+def get_bioc_species_name(): + parts = config["ref"]["species"].split("_") + if len(parts) != 2: + raise ValueError("`ref.species` must be of the form 'Genus_species'") + genus, species = parts + return f"{genus[0].upper()}{species.lower()}"
📜 Review details
Configuration used: CodeRabbit UI
Review profile: CHILL
Plan: Pro
📒 Files selected for processing (12)
.test/config/config.yaml(1 hunks)config/config.yaml(1 hunks)workflow/envs/annotatr.yaml(1 hunks)workflow/envs/biomart.yaml(1 hunks)workflow/envs/rtracklayer.yaml(1 hunks)workflow/rules/circle_map.smk(1 hunks)workflow/rules/common.smk(1 hunks)workflow/rules/ref.smk(1 hunks)workflow/scripts/annotate_cleaned_circles.R(1 hunks)workflow/scripts/clean_circle_map_realign_output.py(1 hunks)workflow/scripts/create_annotation_gff3.R(1 hunks)workflow/scripts/create_transcripts_to_genes_mappings.R(1 hunks)
🧰 Additional context used
🪛 YAMLlint (1.37.1)
.test/config/config.yaml
[error] 15-15: trailing spaces
(trailing-spaces)
workflow/envs/biomart.yaml
[error] 7-7: no new line character at the end of file
(new-line-at-end-of-file)
workflow/envs/annotatr.yaml
[error] 7-7: no new line character at the end of file
(new-line-at-end-of-file)
workflow/envs/rtracklayer.yaml
[error] 7-7: no new line character at the end of file
(new-line-at-end-of-file)
🔇 Additional comments (11)
workflow/envs/rtracklayer.yaml (1)
1-4: Verify thenodefaultschannel usage.
Confirm that your Conda parser supportsnodefaultsto disable default channels, or consider using theno_default_channels: truedirective for compatibility.workflow/envs/annotatr.yaml (1)
1-4: Verify thenodefaultschannel usage.
Make surenodefaultsis recognized by your environment resolver or switch to the standardno_default_channelsflag if needed.workflow/envs/biomart.yaml (1)
1-4: Verify thenodefaultschannel usage.
Double-check thatnodefaultsis a valid mechanism to disable defaults in your Conda setup, or preferno_default_channels: truefor clarity.config/config.yaml (1)
12-12: Updated Ensembl release to 112.
The release version was bumped to112to align with new annotation rules and scripts—please confirm compatibility across downstream Snakemake rules and R scripts..test/config/config.yaml (1)
12-12: Align test config with main release version.
Test configuration now usesrelease: 112; ensure test fixtures and CI pipelines are updated accordingly.workflow/scripts/create_annotation_gff3.R (4)
1-4: LGTM!Proper logging setup for capturing both standard output and error messages.
5-9: LGTM!Good practice to enable global error tracing with
rlang::global_entrace()for better debugging capabilities.
84-115: LGTM!Well-structured data transformation with proper handling of empty values and filtering of incomplete entries.
193-196: LGTM!Proper sorting of genomic ranges before export ensures consistent output.
workflow/rules/ref.smk (2)
116-116:❓ Verification inconclusive
Verify wrapper version compatibility.
The new rules use wrapper version
v3.13.6while existing rules usev1.21.x. This significant version jump may introduce compatibility issues.Please verify that:
- The newer wrapper version is compatible with your Snakemake version
- The version jump is intentional and tested
- Consider updating all wrapper versions consistently across the workflow
Also applies to: 130-130
🏁 Script executed:
#!/bin/bash # Check all wrapper versions in the workflow echo "Wrapper versions used in the workflow:" rg -A1 'wrapper:' --no-heading | grep -E 'v[0-9]+\.[0-9]+\.[0-9]+' | sort | uniq -cLength of output: 1337
Verify wrapper version consistency and compatibility
The workflow currently mixes wrapper versions across rules:
- workflow/rules/ref.smk (lines 116–117):
"v3.13.6/bio/reference/ensembl-annotation""v3.13.6/bio/reference/ensembl-regulation"- workflow/rules/datavzrd.smk:
"v3.8.0/utils/datavzrd"- All other rules (circle_map.smk, mapping.smk, trimming.smk, utils.smk, ref.smk):
- v1.21.x wrappers (e.g.,
"v1.21.2/bio/samtools/sort","v1.21.1/bio/bwa/index")Please confirm that:
- Your Snakemake version supports wrappers from both the v3.x and v1.x branches.
- The jump to v3.13.6 for the Ensembl annotation/regulation wrappers is intentional and has been tested.
- You’ve evaluated whether you want to standardize on a single wrapper major version (or explicitly pin mixed versions) to reduce maintenance risk.
158-159:⚠️ Potential issueFix missing comma after log file path.
Missing comma after the log file path will cause a Snakemake syntax error.
log: - "logs/all_annotations.harmonized.gff3.log" + "logs/all_annotations.harmonized.gff3.log", conda:Likely an incorrect or invalid review comment.
| # turn int cols into int | ||
| circles.loc[:, int_cols] = circles.loc[:, int_cols].round(0).applymap(lambda v: int(v) if not pd.isna(v) else pd.NA) | ||
| circles[int_cols] = circles[int_cols].round(0).map(lambda v: int(v) if not pd.isna(v) else pd.NA) | ||
|
|
There was a problem hiding this comment.
.map() on a DataFrame raises AttributeError – use .applymap() or column-wise logic
DataFrame objects do not implement .map; only Series do.
This line will therefore crash before any output is written.
-circles[int_cols] = circles[int_cols].round(0).map(lambda v: int(v) if not pd.isna(v) else pd.NA)
+# convert each value to integer while keeping NA
+circles[int_cols] = (
+ circles[int_cols]
+ .round(0)
+ .applymap(lambda v: pd.NA if pd.isna(v) else int(v))
+)📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| # turn int cols into int | |
| circles.loc[:, int_cols] = circles.loc[:, int_cols].round(0).applymap(lambda v: int(v) if not pd.isna(v) else pd.NA) | |
| circles[int_cols] = circles[int_cols].round(0).map(lambda v: int(v) if not pd.isna(v) else pd.NA) | |
| # turn int cols into int | |
| # convert each value to integer while keeping NA | |
| circles[int_cols] = ( | |
| circles[int_cols] | |
| .round(0) | |
| .applymap(lambda v: pd.NA if pd.isna(v) else int(v)) | |
| ) |
🤖 Prompt for AI Agents
In workflow/scripts/clean_circle_map_realign_output.py around lines 35 to 37,
the code incorrectly uses .map() on a DataFrame, which causes an AttributeError
because .map() is only available on Series. To fix this, replace .map() with
.applymap() to apply the lambda function element-wise across the DataFrame, or
alternatively apply the lambda function column-wise to each Series in int_cols.
workflow/rules/circle_map.smk
Outdated
| rule annotate_cleaned_circles: | ||
| input: | ||
| all_annotations="resources/all_annotations.harmonized.gff3.gz", | ||
| tsv="results/circle-map/{sample}.circles.cleaned.tsv", | ||
| output: | ||
| tsvs=directory( | ||
| "results/circle-map/{sample}.circles.cleaned.annotated/" | ||
| ), | ||
| log: | ||
| "logs/circle-map/{sample}.circles.cleaned.annotated.logs", | ||
| conda: | ||
| "../envs/annotatr.yaml" | ||
| params: | ||
| build=config["ref"]["build"], | ||
| script: | ||
| "../scripts/annotate_cleaned_circles.R" |
There was a problem hiding this comment.
Output directory declared, but the R script never writes to it
annotate_cleaned_circles.R creates files in the current working directory via file.path(...), yet the rule’s declared output is
results/circle-map/{sample}.circles.cleaned.annotated/.
Snakemake will mark the rule as incomplete because the declared directory is still empty.
Pass the directory path to the script and create it there, e.g.:
output:
- tsvs=directory("results/circle-map/{sample}.circles.cleaned.annotated/")
+ outdir=directory("results/circle-map/{sample}.circles.cleaned.annotated/")
...
params:
build=config["ref"]["build"],
+ outdir=lambda wc: f"results/circle-map/{wc.sample}.circles.cleaned.annotated",
script:
"../scripts/annotate_cleaned_circles.R"and in the R script use snakemake@params[["outdir"]] + dir.create(outdir, recursive=TRUE).
📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| rule annotate_cleaned_circles: | |
| input: | |
| all_annotations="resources/all_annotations.harmonized.gff3.gz", | |
| tsv="results/circle-map/{sample}.circles.cleaned.tsv", | |
| output: | |
| tsvs=directory( | |
| "results/circle-map/{sample}.circles.cleaned.annotated/" | |
| ), | |
| log: | |
| "logs/circle-map/{sample}.circles.cleaned.annotated.logs", | |
| conda: | |
| "../envs/annotatr.yaml" | |
| params: | |
| build=config["ref"]["build"], | |
| script: | |
| "../scripts/annotate_cleaned_circles.R" | |
| rule annotate_cleaned_circles: | |
| input: | |
| all_annotations="resources/all_annotations.harmonized.gff3.gz", | |
| tsv="results/circle-map/{sample}.circles.cleaned.tsv", | |
| output: | |
| outdir=directory("results/circle-map/{sample}.circles.cleaned.annotated/"), | |
| log: | |
| "logs/circle-map/{sample}.circles.cleaned.annotated.logs", | |
| conda: | |
| "../envs/annotatr.yaml" | |
| params: | |
| build=config["ref"]["build"], | |
| outdir=lambda wc: f"results/circle-map/{wc.sample}.circles.cleaned.annotated", | |
| script: | |
| "../scripts/annotate_cleaned_circles.R" |
🤖 Prompt for AI Agents
In workflow/rules/circle_map.smk around lines 81 to 96, the rule declares an
output directory results/circle-map/{sample}.circles.cleaned.annotated/ but the
R script annotate_cleaned_circles.R does not write files there, causing
Snakemake to see the output as missing. Fix this by passing the output directory
path as a parameter to the R script via params (e.g., add outdir=params.build or
a new param outdir with the output directory path), then modify the R script to
read this parameter using snakemake@params[["outdir"]], create the directory if
needed with dir.create(outdir, recursive=TRUE), and write all output files into
this directory.
| group_walk( | ||
| ~ write_tsv( | ||
| .x, | ||
| file = file.path( | ||
| str_c( | ||
| str_replace_all(.y$circle_region, "[:-]", "_"), | ||
| ".tsv" | ||
| ) | ||
| ) | ||
| ) | ||
| ) |
There was a problem hiding this comment.
Files written outside the declared output directory
The call
file = file.path(str_c(str_replace_all(.y$circle_region, "[:-]", "_"), ".tsv"))creates files in getwd(). The Snakemake rule expects them under the output directory.
Add:
outdir <- snakemake@params[["outdir"]]
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)
...
file = file.path(outdir,
str_c(str_replace_all(.y$circle_region, "[:-]", "_"), ".tsv"))🤖 Prompt for AI Agents
In workflow/scripts/annotate_cleaned_circles.R around lines 98 to 108, the code
writes output files directly to the working directory instead of the Snakemake
output directory. To fix this, retrieve the output directory path from
snakemake@params[["outdir"]], create the directory if it doesn't exist using
dir.create with recursive=TRUE and showWarnings=FALSE, and update the file path
in write_tsv to write files inside this output directory by prepending outdir to
the file path.
| circles_gr <- GRanges( | ||
| seqnames = pull(circles, region), | ||
| ecDNA_status = "ecDNA" | ||
| ) | ||
|
|
||
| genome(circles_gr) <- genome_build | ||
|
|
There was a problem hiding this comment.
GRanges construction is missing start/end coordinates → runtime error
GRanges() requires a ranges argument (IRanges). Supplying only seqnames will throw
“argument "ranges" is missing, with no default”.
Parse the region column first:
-regions = pull(circles, region),
-ecDNA_status = "ecDNA"
+seqnames = sub(":.*", "", circles$region),
+ranges = IRanges(
+ start = as.integer(sub(".*:(\\d+)-.*", "\\1", circles$region)),
+ end = as.integer(sub(".*-(\\d+)$", "\\1", circles$region))
+),
+ecDNA_status = "ecDNA"Committable suggestion skipped: line range outside the PR's diff.
🤖 Prompt for AI Agents
In workflow/scripts/annotate_cleaned_circles.R around lines 17 to 23, the
GRanges object is constructed without specifying the required start and end
coordinates, causing a runtime error. Extract the start and end positions from
the 'region' column (which likely contains genomic coordinates) and create an
IRanges object to pass as the 'ranges' argument in GRanges. This will fix the
missing argument error and correctly define the genomic ranges.
| if (wanted_build == "GRCh37") { | ||
| grch <- "37" | ||
| version <- NULL | ||
| cli_warn(c( | ||
| "As you specified build 'GRCH37' in your configuration yaml, biomart forces", | ||
| "us to ignore the release you specified ('{release}')." | ||
| )) | ||
| } else { | ||
| grch <- NULL | ||
| version <- wanted_release | ||
| } |
There was a problem hiding this comment.
cli_warn() references undefined variable {release}
{release} is not in scope (should be {wanted_release}), causing glue evaluation to fail.
-cli_warn(c(
- "As you specified build 'GRCH37' in your configuration yaml, biomart forces",
- "us to ignore the release you specified ('{release}')."
-))
+cli_warn(c(
+ "As you specified build 'GRCh37' in your configuration yaml, biomart forces",
+ "us to ignore the release you specified ('{wanted_release}')."
+))📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| if (wanted_build == "GRCh37") { | |
| grch <- "37" | |
| version <- NULL | |
| cli_warn(c( | |
| "As you specified build 'GRCH37' in your configuration yaml, biomart forces", | |
| "us to ignore the release you specified ('{release}')." | |
| )) | |
| } else { | |
| grch <- NULL | |
| version <- wanted_release | |
| } | |
| if (wanted_build == "GRCh37") { | |
| grch <- "37" | |
| version <- NULL | |
| cli_warn(c( | |
| "As you specified build 'GRCh37' in your configuration yaml, biomart forces", | |
| "us to ignore the release you specified ('{wanted_release}')." | |
| )) | |
| } else { | |
| grch <- NULL | |
| version <- wanted_release | |
| } |
🤖 Prompt for AI Agents
In workflow/scripts/create_transcripts_to_genes_mappings.R around lines 16 to
26, the cli_warn() call references an undefined variable {release} which causes
glue evaluation to fail. Replace {release} with the correctly scoped variable
{wanted_release} in the warning message to fix the issue.
| if (retrieved_build != build) { | ||
| cli_abort(c( | ||
| "Ensembl release and genome build number in your configuration are not compatible.", | ||
| "x" = "Genome build '{build}' not available via biomart for Ensembl release '{release}''.", | ||
| "i" = "Ensembl release '{release}' only provides build '{retrieved_build}'.", | ||
| " " = "Please fix your configuration yaml file's reference entry, you have two options:", | ||
| "*" = "Change the build entry to '{retrieved_build}'.", | ||
| "*" = "Change the release entry to one that provides build '{build}'. You have to determine this from biomart by yourself." | ||
| )) |
There was a problem hiding this comment.
🛠️ Refactor suggestion
cli_abort() strings also reference missing variable release
Inside get_mart() only version is available, not release; this will error at runtime.
- "x" = "Genome build '{build}' not available via biomart for Ensembl release '{release}''.",
- "i" = "Ensembl release '{release}' only provides build '{retrieved_build}'.",
+ "x" = "Genome build '{build}' not available via biomart for Ensembl release '{version}'.",
+ "i" = "Ensembl release '{version}' only provides build '{retrieved_build}'.",🤖 Prompt for AI Agents
In workflow/scripts/create_transcripts_to_genes_mappings.R around lines 42 to
50, the cli_abort() call references a variable named release which is not
defined or available in this scope; only version is accessible inside
get_mart(). To fix this, replace all instances of release in the cli_abort()
message with version or ensure the correct variable holding the Ensembl release
number is used consistently to avoid runtime errors.
| dplyr::rename( | ||
| ensemble_id = id | ||
| ) |> |
There was a problem hiding this comment.
Fix typo in column name.
The column name should be ensembl_id instead of ensemble_id. Ensembl is the name of the genomic database.
dplyr::rename(
- ensemble_id = id
+ ensembl_id = id
) |>🤖 Prompt for AI Agents
In workflow/scripts/create_annotation_gff3.R around lines 188 to 190, there is a
typo in the column rename operation where the column is renamed to "ensemble_id"
instead of the correct "ensembl_id". Change the rename argument from ensemble_id
= id to ensembl_id = id to correctly reflect the genomic database name.
| separate_wider_delim( | ||
| ID, | ||
| delim = ":", | ||
| names = c("id_type", "id") | ||
| ) |> | ||
| dplyr::select(-id_type) |> | ||
| unnest( | ||
| Parent, | ||
| keep_empty = TRUE | ||
| ) |> | ||
| separate_wider_delim( | ||
| Parent, | ||
| delim = ":", | ||
| names = c("parent_type", "parent_id") | ||
| ) |> | ||
| dplyr::rename(name = Name) |
There was a problem hiding this comment.
🛠️ Refactor suggestion
Add validation for ID and Parent format assumptions.
The code assumes IDs and Parents contain a ":" delimiter. Consider adding validation or error handling for cases where this format is not present.
Add validation before the separation:
+ # Validate ID format
+ if (!all(str_detect(.$ID, ":"))) {
+ warning("Some IDs do not contain ':' delimiter, these will be handled incorrectly")
+ }
separate_wider_delim(
ID,
delim = ":",
names = c("id_type", "id")
) |>Committable suggestion skipped: line range outside the PR's diff.
🤖 Prompt for AI Agents
In workflow/scripts/create_annotation_gff3.R around lines 47 to 62, the code
assumes that the ID and Parent fields contain a ":" delimiter for separation,
but this may not always be true. Add validation checks before calling
separate_wider_delim to confirm that these fields contain the ":" delimiter. If
the delimiter is missing, handle the case gracefully by either logging a
warning, skipping the separation, or providing a default value to avoid errors
during processing.
…rkflows/dna-seq-short-read-circle-map into feat/add-circle-annotations
…rkflows/dna-seq-short-read-circle-map into feat/add-circle-annotations
Summary by CodeRabbit
New Features
Chores
Refactor