New module: custom/collectfeaturecounts - #12846
Conversation
Collects one or more subread/featureCounts raw output files (typically one per sample, for the same feature type/ORF caller) into a single long-format table with per-sample TPM computed, dropping zero-count rows. Ported from two now-identical, independently-hardened local modules (nf-core/magmap's and nf-core/metatdenovo's own COLLECT_FEATURECOUNTS) as part of the cross-repo summary-table consolidation tracked in nf-core/magmap#237. Includes the tpm-rounding fix from nf-core/metatdenovo#483/nf-core#484 (rounds tpm to 6 decimal places, since independently-computed tpm values for the same counts can otherwise round the last significant digit differently depending on which R backend does the division). Test data (raw subread/featurecounts fixtures, real per-sample featureCounts.tsv output, not yet present on the modules branch) in nf-core/test-datasets#2247. Generated by Claude Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01UvGYU6kuXJeVv9ffdSzdzH
Joon-Klaps
left a comment
There was a problem hiding this comment.
The R script should not use a globs search but a direct pointer to the input file.
|
|
||
| setDTthreads($task.cpus) | ||
|
|
||
| tibble(f = Sys.glob('*.featureCounts.tsv')) %>% |
There was a problem hiding this comment.
Don't think this should be the way to go.
Should be defined as "$inputfiles"
Also thinking if in general it wouldn't be better to have a Rscript executable stored in templates/collectfeaturecounts.R
Might also be a good example to use https://github.com/nf-core/r-nf-core-utils to facilitate maintanability.
There was a problem hiding this comment.
Thanks for the push on this, @Joon-Klaps — this comment led to a bigger cleanup than the one-line fix it sounds like, and a special thanks is due since it surfaced how suboptimal the state of this code we inherited from our own local pipeline modules actually was.
What changed, in order:
- Read files from the declared
path(inputfiles)input instead ofSys.glob('*.featureCounts.tsv'). - Moved the R body into
templates/collectfeaturecounts.Rper your suggestion. Turned out a Groovy variable computed in thescript:block isn't visible inside atemplatefile (confirmed againstcustom/matrixfilter's own comment about this), so ended up staging every input file into a fixedinput/dir viastageAs: "input/*"(same patternnacho/qcuses) and globbing that from R — no Groovy→template variable passing needed at all. - While re-reading the R itself: the per-file
fread → melt → lazy_dt() → dplyr chain → as_tibble(), stitched back together withtidyr::unnest(), ran N independent small dtplyr query plans instead of one. Consolidated intorbindlist()+ a singlelazy_dt()pipeline over the combined data. - That consolidation removed the only uses of
tidyr/purrr, so dropped both fromenvironment.ymland regenerated the container + conda-lock files to match.
Each step verified with a real (non-stub) Docker nf-test run against the existing snapshot, plus a clean nf-core modules lint.
Still open: your nf-core/r-nf-core-utils suggestion — haven't looked into what it offers yet, will take a look separately rather than block this on it. Let me know if you think it should hold this PR up.
There was a problem hiding this comment.
Ok fair point, won't force you to use that R package.
Sys.glob('*.featureCounts.tsv') was carried over unchanged from the
local COLLECT_FEATURECOUNTS modules this was ported from, and only
worked by coincidence of subread/featurecounts' output naming rather
than using the files the process already declares and stages via
`path(inputfiles)`. Build the R file vector from that input instead,
per review comment on nf-core#12846.
Verified with a real Docker nf-test run; output is byte-identical to
before (existing snapshot matches unchanged).
Generated by Claude
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01QVaft273VdxeEKN4F76S5K
Per review comment on nf-core#12846, also considering the file-list handling. A Groovy variable computed in the `script:` block (like the previous commit's `r_files`) turns out not to be visible inside a `template` file -- confirmed against two existing nf-core/modules precedents (custom/matrixfilter's script: comment says so explicitly; nacho/qc sidesteps the same problem entirely). So instead of passing a file list through Groovy, stage every input file into a fixed `input/` directory (`stageAs: "input/*"`, the same pattern nacho/qc uses) and glob that directory from R -- no Groovy/template variable-passing needed at all, and it keeps working regardless of how many files are passed. `prefix` is likewise reconstructed directly from `$task.ext.prefix`/`$meta.id` inside the template, matching matrixfilter's approach, rather than relying on the Groovy-side `prefix` variable being visible there. Verified with a real Docker nf-test run (existing snapshot matches unchanged) and a clean `nf-core modules lint`. Generated by Claude Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QVaft273VdxeEKN4F76S5K
Each input file was independently pushed through fread -> melt -> lazy_dt() -> a full dplyr chain -> as_tibble(), then the N per-file tibbles were stitched back together with tidyr::unnest(). Verified (https://claude.ai/code/session_01QVaft273VdxeEKN4F76S5K) that: - lazy_dt() was already getting a genuine data.table, contrary to how it read at a glance (the outer tibble(f = ...) was an unrelated per-file iteration driver, never itself passed through dtplyr). - the per-file as_tibble() wasn't dead weight either: lazy_dt()'s dplyr chain is an unevaluated proxy (class dtplyr_step_group) until materialized, and tidyr::unnest() errors ("Input must be list of vectors") if given raw un-materialized elements. So the correctness was fine, but the structure defeats the point of using data.table/dtplyr: N independent small query plans instead of one, and setDTthreads($task.cpus) barely matters when no operation ever spans more than one input file's rows. Now: fread + melt each file (still real per-file data.table work), rbindlist() them into one data.table, then a single lazy_dt() chain over the combined data, materialized once at the very end via as_tibble() (confirmed write_tsv() needs a real data frame -- it rejects a raw dtplyr_step -- but doesn't otherwise care whether that's a tibble or a data.table). tidyr::unnest() and purrr::map() are gone, so library(tidyr) is dropped too (environment.yml still lists r-tidyr/r-purrr; left alone pending a decision on regenerating the container/conda-lock for that separately). Verified with a real Docker nf-test run (existing snapshot matches unchanged, so output is byte-identical) and a clean nf-core modules lint. Generated by Claude Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QVaft273VdxeEKN4F76S5K
Neither library is used in the R code any more since the previous commit removed tidyr::unnest()/purrr::map() in favour of rbindlist(). Rebuilt via 'nf-core modules containers create' (Wave, both amd64 and arm64) and 'nf-core modules containers conda-lock', which updates main.nf's container digest, meta.yml's container/lock-file metadata, and the .conda-lock files to match the pruned environment.yml. The regenerated container ternary in main.nf came out with zero indentation on the '?'/':' lines (still valid Groovy -- nextflow lint passes either way -- but inconsistent with every other module's formatting, e.g. custom/matrixfilter); reindented it by hand to match. Verified with a real Docker nf-test run against the newly built image (existing snapshot matches unchanged) and a clean nf-core modules lint. Generated by Claude Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QVaft273VdxeEKN4F76S5K
Co-authored-by: Joon Klaps <joon.klaps@kuleuven.be>
Co-authored-by: Joon Klaps <joon.klaps@kuleuven.be>
Joon-Klaps
left a comment
There was a problem hiding this comment.
The stub versions, should be resolved not harcoded.
The R nf-core package does allow you to do this easily
| "${task.process}": | ||
| R: 4.1.0 | ||
| dplyr: 1.0.7 | ||
| readr: 2.0.0 | ||
| stringr: 1.4.0 | ||
| dtplyr: 1.1.0 | ||
| data.table: 1.14.0 | ||
| END_VERSIONS |
There was a problem hiding this comment.
Think this should be resolved and not hard coded.
There was a problem hiding this comment.
Fair, fixed — the stub's versions.yml was hardcoding stale literal numbers left over from the original local module (R: 4.1.0 etc.), unrelated to what's actually pinned in environment.yml. Now resolves each value via Rscript -e ... at run time, same pattern custom/matrixfilter's stub already uses for its r-base line.
Verified with a real Docker nf-test run — resolved values now match environment.yml's pins exactly (R 4.5.3, dplyr 1.2.1, readr 2.2.0, stringr 1.6.0, dtplyr 1.3.3, data.table 1.17.8) — and updated the stub test's snapshot accordingly.
Per review comment on nf-core#12846: the stub block's versions.yml hardcoded stale literal version numbers left over from the original local module, unrelated to what's actually pinned in environment.yml (and already wrong -- e.g. R: 4.1.0 vs the container's actual 4.5.3). Query each value live via `Rscript -e ...` command substitution instead, same pattern custom/matrixfilter's stub already uses for its own r-base version line. The real (non-stub) script already did this correctly via packageVersion() in the template, so only the stub needed fixing. Verified with a real Docker nf-test run: resolved values now match environment.yml's actual pins exactly (R 4.5.3, dplyr 1.2.1, readr 2.2.0, stringr 1.6.0, dtplyr 1.3.3, data.table 1.17.8). Updated the stub test's snapshot accordingly (the only thing that changed is the versions.yml content, now correct instead of stale). Clean nf-core modules lint and nextflow lint. Generated by Claude Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QVaft273VdxeEKN4F76S5K
|
Thanks @Joon-Klaps |
PR checklist
nf-core modules lintpasses clean.Description
Collects one or more
subread/featurecountsraw output files (typically one per sample, for the same feature type/ORF caller) into a single long-format table (orf,chr,start,end,strand,length,sample,count,tpm), computing per-sample TPM and dropping zero-count rows.Ported from two independently-hardened, now-functionally-identical local modules — nf-core/magmap's and nf-core/metatdenovo's own
COLLECT_FEATURECOUNTS— as part of a cross-repo effort to share their summary-table-generation code (tracked in nf-core/magmap#237). Both pipelines' local copies were first made caller/pipeline-agnostic by moving their one pipeline-specific ID-reconciliation special-case (magmap's genome-accession join, metatdenovo's Transdecodercds.-prefix stripping) into small local post-processing steps — see nf-core/magmap#238 and nf-core/metatdenovo#483.tpmis rounded to 6 decimal places. This carries over a fix from nf-core/metatdenovo#483/#484: when the same underlying counts get independently re-computed elsewhere (e.g. after a caller-consolidation step using a different R backend —dtplyr/data.tablehere vs plaindplyrthere), the unrounded double can round its last significant digit differently for the same mathematical value, which made an exact-equality test downstream intermittently fail. Rounding here keeps output stable regardless of which backend a consumer uses to recompute the same value.Testing
Verified with real Docker runs against two new fixtures (two samples, mixed zero/non-zero counts, one shared multi-exon gene) — nf-core/test-datasets#2247.
tpmcorrectly sums to exactly1e6per sample after rounding, zero-count rows are dropped, sample names are correctly stripped of their.sorted.bamsuffix.Generated by Claude