revisit downsample as a better SM workflow #173
Conversation
📝 Walkthrough📝 WalkthroughWalkthroughThe pull request includes the deletion of three Snakemake workflow files and a command-line interface module, which collectively facilitated genomic data processing and HPC job submission. Additionally, a new utility for extracting barcodes from SAM/BAM files is introduced, alongside significant modifications to the downsampling function and the addition of a new Snakemake workflow for downsampling. Updates to GitHub Actions workflow configuration and dependency management in YAML files are also included. Changes
Possibly related PRs
Suggested labels
Thank you for using CodeRabbit. We offer it for free to the OSS community and would appreciate your support in helping us grow. If you find it useful, would you consider giving us a shout-out on your favorite social media? 🪧 TipsChatThere are 3 ways to chat with CodeRabbit:
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: 1
🧹 Outside diff range and nitpick comments (3)
harpy/bin/extract_bxtags.py (1)
37-40: Simplify nestedifstatementsYou can simplify the nested
ifstatements by combining the conditions into a single line.Here's the suggested change:
try: barcode = record.get_tag(bx_tag) - if invalid_pattern.search(barcode): - # invalid barcode retention - if rng.random() > inv_prop: - continue + if invalid_pattern.search(barcode) and rng.random() > inv_prop: + continue🧰 Tools
🪛 Ruff (0.8.0)
37-39: Use a single
ifstatement instead of nestedifstatements(SIM102)
harpy/snakefiles/sv_leviathan_pop.smk (1)
Line range hint
69-70: Adjust thread allocation to prevent over-allocationUsing
max(10, workflow.cores)may allocate more threads than available cores, leading to resource contention. Consider usingmin(10, workflow.cores)to ensure the number of threads does not exceed the available cores.Here's the suggested change:
threads: - max(10, workflow.cores) + min(10, workflow.cores).github/workflows/tests.yml (1)
836-836: Consider using a more flexible file pattern.The current pattern
test/fastq/sample1.*gzonly matches gzipped files. While this works for the test case, consider using the original patterntest/fastq/sample1.*to maintain flexibility for both gzipped and non-gzipped files, especially since the PR aims to enhance format support.- run: harpy downsample -d 1 --quiet test/fastq/sample1.*gz + run: harpy downsample -d 1 --quiet test/fastq/sample1.*
📜 Review details
Configuration used: CodeRabbit UI
Review profile: CHILL
📒 Files selected for processing (11)
.deprecated/align-minimap.smk(0 hunks).deprecated/hpc.py(0 hunks).deprecated/phase.bak.smk(0 hunks).github/workflows/tests.yml(1 hunks)harpy/bin/extract_bxtags.py(1 hunks)harpy/downsample.py(2 hunks)harpy/snakefiles/downsample.smk(1 hunks)harpy/snakefiles/sv_leviathan.smk(1 hunks)harpy/snakefiles/sv_leviathan_pop.smk(1 hunks)resources/harpy.yaml(1 hunks)resources/meta.yaml(1 hunks)
💤 Files with no reviewable changes (3)
- .deprecated/phase.bak.smk
- .deprecated/align-minimap.smk
- .deprecated/hpc.py
🧰 Additional context used
🪛 Ruff (0.8.0)
harpy/bin/extract_bxtags.py
37-39: Use a single if statement instead of nested if statements
(SIM102)
🔇 Additional comments (6)
resources/harpy.yaml (1)
9-9: Consider adding version constraints for lrez
Other dependencies like bcftools, conda, and samtools have specific version constraints. Consider adding version constraints for lrez to ensure reproducibility and prevent potential compatibility issues.
Let's verify the package availability and latest version:
resources/meta.yaml (2)
Line range hint 38-46: Dependencies are consistently declared
The run requirements mirror the dependencies in harpy.yaml, maintaining consistency across configuration files.
39-39: Verify lrez usage in the downsampling workflow
The addition of lrez aligns with the PR's objective to improve the downsampling workflow. However, let's verify its integration:
✅ Verification successful
Let me check the content of the downsampling workflow file to verify the lrez usage.
lrez is correctly added as a requirement and properly integrated
The verification confirms that lrez is actively used in the downsampling workflow (harpy/snakefiles/downsample.smk) for both BAM and FASTQ processing:
- Used for indexing files via
LRez indexcommand - Used for querying and extracting data via
LRez querycommand - Also utilized in other workflows (sv_leviathan.smk and sv_leviathan_pop.smk)
🏁 Scripts executed
The following scripts were executed for the analysis:
Script:
#!/bin/bash
# Description: Check how lrez is used in the downsampling workflow
# Look for lrez imports or usage in Python files
rg -l 'import.*lrez|from.*lrez'
# Look for lrez in Snakemake rules
rg -l 'lrez' --type snakemake
# Look for workflow files containing downsample
fd -e smk -e snakefile -e snake -X rg -l 'downsample'
Length of output: 185
Script:
#!/bin/bash
# Check the content of downsample.smk for lrez usage
cat harpy/snakefiles/downsample.smk
# Also search for any other files mentioning lrez to ensure complete coverage
rg -i "lrez"
Length of output: 6720
harpy/downsample.py (1)
29-36: Function parameters updated appropriately
The updates to the function parameters and command-line options accurately reflect the new workflow and enhance the script's flexibility.
harpy/snakefiles/sv_leviathan_pop.smk (1)
111-112: Verify availability of required tools when container is set to None
By setting container: None, this rule will not use any container or conda environment. Ensure that all necessary software dependencies for this rule are available in the execution environment.
harpy/snakefiles/sv_leviathan.smk (1)
64-65: LGTM: Container directive change aligns with best practices.
The change from conda to container directive is appropriate, utilizing the default container specified at the top of the file (containerized: "docker://pdimens/harpy:latest").
✅ Verification successful
LGTM: Container configuration is correctly implemented
The verification confirms that:
- The file correctly sets
containerized: "docker://pdimens/harpy:latest"at the top level - Setting
container: Nonein theindex_barcoderule is intentional, as this pattern is consistently used across all rules in the workflow - All Snakemake files in the harpy project follow the same containerization pattern with the top-level
containerizeddirective
🏁 Scripts executed
The following scripts were executed for the analysis:
Script:
#!/bin/bash
# Verify container usage consistency across rules
rg -A 1 'rule.*:$' | rg -B 1 'container:'
Length of output: 41
Script:
#!/bin/bash
# Let's first check if the file exists and view its content
cat harpy/snakefiles/sv_leviathan.smk
Length of output: 7474
Script:
#!/bin/bash
# Search for container directives in all Snakemake files
fd -e smk -e snakefile -I . -x rg -B2 -A2 'container:'
Length of output: 12207
Script:
#!/bin/bash
# Search for containerized directive at the top level
rg 'containerized:' -B2 -A2
Length of output: 3796
There was a problem hiding this comment.
Actionable comments posted: 0
🧹 Outside diff range and nitpick comments (3)
harpy/reports/hapcut.Rmd (3)
54-60: Add input validation and optimize the NX function.The NX function implementation is correct but could be enhanced for robustness and efficiency.
Consider these improvements:
NX <- function(lengths, X=50){ + if(length(lengths) == 0) return(0) + if(X < 0 || X > 100) stop("X must be between 0 and 100") + if(!is.numeric(lengths)) stop("lengths must be numeric") lengths <- as.numeric(sort(lengths, decreasing=TRUE)) index <- which(cumsum(lengths) / sum(lengths) >= X/100)[1] + if(is.na(index)) return(0) # Handle case when threshold is never reached return(lengths[index]) }
Line range hint
167-180: Consider adding N75 and N90 metrics to per-contig view.For consistency with the overall statistics, consider including N75 and N90 metrics in the per-contig table.
Here's how to modify the code:
percontig <- df %>% group_by(contig) %>% summarise( n_haplo = round(length(n_snp)), mean_snps = round(mean(n_snp), digits = 0), median_snps = median(n_snp), mean_blocksize = round(mean(block_length), digits = 0), N50 = NX(block_length, 50), + N75 = NX(block_length, 75), + N90 = NX(block_length, 90), max_blocksize = max(block_length) ) DT::datatable( percontig, rownames = F, extensions = 'Buttons', options = list( dom = 'Brtip', buttons = list(list(extend = "csv",filename = "phasing_per_contig")), scrollX = TRUE ), - colnames = c("Contig", "Total Haplotypes", "Mean SNPs", "Median SNPs", "Mean Haplotype Length", "Haplotype N50", "Largest Haplotype"), + colnames = c("Contig", "Total Haplotypes", "Mean SNPs", "Median SNPs", "Mean Haplotype Length", "Haplotype N50", "Haplotype N75", "Haplotype N90", "Largest Haplotype"), fillContainer = T )
Line range hint
263-276: Consider adding N75 and N90 metrics to per-sample view.For consistency with the overall statistics and to maintain parity with the suggested per-contig improvements, consider including N75 and N90 metrics in the per-sample table.
Here's how to modify the code:
persample <- df %>% group_by(sample) %>% summarise( n_haplo = sum(length(n_snp)), mean_snps = round(mean(n_snp), digits = 0), median_snps = median(n_snp), mean_blocksize = round(mean(block_length), digits = 0), N50 = NX(block_length, 50), + N75 = NX(block_length, 75), + N90 = NX(block_length, 90), max_blocksize = max(block_length) ) DT::datatable( persample, rownames = F, extensions = 'Buttons', options = list( dom = 'Brtip', buttons = list(list(extend = "csv",filename = "phasing_per_sample")), scrollX = TRUE ), - colnames = c("Sample", "Haplotypes", "Mean SNPs", "Median SNPs", "Mean Haplotype Length", "Haplotype N50", "Largest Haplotype"), + colnames = c("Sample", "Haplotypes", "Mean SNPs", "Median SNPs", "Mean Haplotype Length", "Haplotype N50", "Haplotype N75", "Haplotype N90", "Largest Haplotype"), fillContainer = T )
📜 Review details
Configuration used: CodeRabbit UI
Review profile: CHILL
📒 Files selected for processing (1)
harpy/reports/hapcut.Rmd(7 hunks)
🔇 Additional comments (2)
harpy/reports/hapcut.Rmd (2)
81-84: LGTM! Improved metrics for haplotype length distribution.
The addition of N50, N75, and N90 metrics provides a more comprehensive view of haplotype length distribution, which is standard practice in genomics analysis.
102-114: LGTM! Well-structured UI elements for new metrics.
The value boxes for N50, N75, and N90 are consistently styled and properly formatted, maintaining the dashboard's visual hierarchy.
This reluctantly changes
downsamplefrom a self-contained script to a proper snakemake workflow that has better support for fastq vs bam and a much more sensible method of invalid barcode handlingSummary by CodeRabbit
New Features
Bug Fixes
Chores
lrezto configuration files.