Repository navigation
Port rMATS to the nf-core rmats/prep module and refactor RMATS_POST to the nf-core module structure - #285
Conversation
Replace the local RMATS_PREP module with the nf-core rmats/prep module. It
preps one BAM file at a time, so each BAM is read once whatever the number
of contrasts it takes part in, instead of every BAM of a contrast in one
go. rMATS matches the .rmats files to the --b1/--b2 lists by BAM file name
and never reads the BAMs again in the post step, which is what makes the
split possible: the RMATS subworkflow gathers the .rmats files of the
samples of each contrast by sample id and hands them to RMATS_POST with
the bam lists.
Refactor RMATS_POST to the nf-core module template: [ meta, rmats,
bam_list1, bam_list2 ], [ meta2, gtf ] and read_length inputs, per file
type emits, environment.yml, meta.yml, a stub and nf-tests (unpaired,
paired, single condition, stub). --cstat, --paired-stats and --novelSS
with --mil/--mel move to ext.args in modules.config, together with the
library options of the prep step, and the paired model is marked by the
prefix instead of a nested folder. --statoff stays in the module, since
it follows from the absence of a second bam list. The module also reports
the PAIRADISE version.
The .rmats files are published to rmats/prep/{sample}.rmats and the
results of a contrast to rmats/{contrast}{_paired}/, with the log next to
the folder. The intermediate tmp/ folder of the post step is no longer
published. The pipeline snapshots change accordingly; the hashes of the
rMATS result tables are unchanged.
Generated by Claude Opus 5
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019SiaVrAVHBtxZcGgoiwQun
|
erikrikarddaniel
left a comment
There was a problem hiding this comment.
Claude's comments:
I checked the part the pipeline snapshots can't cover. The test data is paired-end and unstranded, which are rMATS's defaults for -t and --libType. So unchanged hashes alone can't show that post is safe without -t, --libType, --variable-read-length and --allow-clipping. I ran rmats.py 4.3.0 from the module's container on two of the chrX test BAMs, with prep set to stranded paired-end and, separately, to single-end. In both cases post's output was byte-identical with and without those four options (about 1500 events). The prep settings do change the result (1499 vs 1520 events), so the options sit on the right step.
One question inline, about duplicate sample ids with --source genome_bam.
A small naming thought, not blocking: modules/local/rmats/post would sit next to modules/nf-core/rmats/prep and match the nf-core name if RMATS_POST ever goes upstream.
| // The post step of a contrast takes the `.rmats` files of every sample in its bam | ||
| // lists. The bam lists carry the BAM file names, which is what rMATS matches the | ||
| // `.rmats` files by, so the samples of a contrast are looked up by sample id here | ||
| ch_rmats_by_sample = RMATS_PREP.out.rmats |
There was a problem hiding this comment.
Can a --source genome_bam samplesheet have two rows with the same sample? The comment at line 33 mentions technical replicates sharing a sample id, and assets/schema_input_genome_bam.json has no uniqueEntries. Unlike the fastq source, the genome_bam branch of PIPELINE_INITIALISATION doesn't groupTuple either.
If it can, each BAM gets its own RMATS_PREP task with the same prefix, so both write ${meta.id}.rmats. Here both then match the same id, and RMATS_POST stages them into rmats_tmp/*, which fails. I checked that part with a minimal process:
Process `P` input file name collision -- There are multiple input files for each of the following file names: rmats_tmp/S1.rmats
The two published rmats/prep/{sample}.rmats would also overwrite each other. The old all-in-one prep took both BAMs in one list, so it didn't hit this. I didn't run the pipeline with such a samplesheet.
If duplicates aren't meant to be supported, "uniqueEntries": ["sample"] in the schema would say so up front. If they are, naming the prep output after the BAM instead of the sample would keep them apart, e.g. ext.prefix = { genome_bam.baseName } on RMATS_PREP. The bam lists already carry the BAM names that rMATS matches on.
There was a problem hiding this comment.
Good catch, and the answer is that duplicates were never viable for this source, so I went with the first option. With --source genome_bam every BAM goes through BAM_SORT_STATS_SAMTOOLS first, whose prefix is ${meta.id}_sorted (conf/modules.config), so two rows of the same sample already produced two S1_sorted.bam: the old all-in-one prep staged both into one task (the same name collision you saw) and CREATE_BAMLIST would have listed the name twice, which rMATS rejects as duplicate input bam files. That also rules out naming the prep output after the BAM, since by then the BAM is already named after the sample.
e91fe6f adds "uniqueEntries": ["sample"] to the genome_bam, transcriptome_bam and salmon_results schemas, the three sources without a merge step, and says so in docs/usage.md. Checked with a fifth row repeating ERR188383 on the test_genome_bam samplesheet:
The following errors have been detected in dup.csv:
-> Entry 5: Detected duplicate entries: [sample:ERR188383]
The three source tests still pass with unique names.
Only the fastq source merges the rows of a sample. With --source
genome_bam two rows of the same sample already collided in
BAM_SORT_STATS_SAMTOOLS, which names the sorted BAM after the sample id,
and from there in CREATE_BAMLIST and rMATS's own duplicate BAM check. The
per sample RMATS_PREP now adds a `${meta.id}.rmats` collision to that
list, so the schemas of the genome_bam, transcriptome_bam and
salmon_results sources declare `uniqueEntries` on `sample` and the
samplesheet fails validation up front instead. Documented in usage.md.
Raised in review of nf-core#285.
Generated by Claude Opus 5
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019SiaVrAVHBtxZcGgoiwQun
|
Thanks for checking the prep/post option split against stranded and single-end data, that was the part the unstranded paired test set could not tell. On the name: agreed that Generated by Claude Opus 5 |
|
thanks @erikrikarddaniel ! |
Summary
Port the rMATS modules to the nf-core structure. There is an nf-core
rmats/prepmodule but normats/post, so:RMATS_PREPis replaced by the nf-corermats/prepmodule, installed withnf-core modules installand left untouched. It preps one BAM file at a time, so each BAM is read once whatever the number of contrasts it takes part in, instead of every BAM of a contrast in one go. TheRMATSsubworkflow gathers the.rmatsfiles of the samples of each contrast by sample id and hands them to the post step with the bam lists.RMATS_POSTis refactored to the nf-core module template:[ meta, rmats, bam_list1, bam_list2 ],[ meta2, gtf ]andread_lengthinputs, per file type emits (mats,from_gtf,raw_input,summary,log),environment.yml,meta.yml, a stub and nf-tests. It also reports thePAIRADISEversion.--cstat,--paired-statsand--novelSSwith--mil/--melmove toext.argsinmodules.config, together with the library options of the prep step (-t,--libType,--variable-read-length,--allow-clipping).--statoffstays in the module, since it follows from the absence of a second bam list.Why the split is safe
rmats.py --task postnever reads the BAM files. It walks--tmpfor.rmatsfiles and matches each one to the--b1/--b2lists by the BAM file name recorded on its first line (split_sg_files_by_baminrmatspipeline.pyx), ignoring.rmatsfiles of BAMs that are not in the lists.CREATE_BAMLISTand the nf-core prep module both write the staged file name, so they agree. The hashes of the rMATS result tables in the pipeline snapshots are unchanged by the port.Output layout
rmats/prep/{sample}.rmatsandrmats/prep/{sample}_read_outcomes_by_bam.txt(wasrmats/{contrast}/rmats_temp/*with timestamped names, plus armats_prep.log)rmats/{contrast}{_paired}/*andrmats/{contrast}{_paired}.log(wasrmats/{contrast}/rmats_post{_paired}/*and.log)tmp/folder of the post step is no longer publisheddocs/output.md,tests/.nftignoreand thestable_pathignore lists of the pipeline tests are updated accordingly.One thing to keep an eye on
The nf-core
rmats/prepmodule isprocess_single(1 CPU, 6 GB), where the old all-in-oneRMATS_PREPwasprocess_high. One BAM per task needs far less than all of them at once, but a large BAM may still want more memory; that is awithName: RMATS_PREPresource override inmodules.configif it turns out to be needed.Testing
nf-test test modules/local/rmats_post --profile +docker: 4/4 pass (unpaired, paired, single condition, stub), snapshot stable on rerundefault,genome_bamandmultiple_runsregenerated; the rMATS result table hashes are the same as beforenf-test test tests/suite: 7/7 passnf-core pipelines lint: 0 failures;nf-core modules lint --local: only the generic warnings every local module here gets;prekcleanGenerated by Claude Opus 5
PR checklist
nf-core pipelines lint).nf-test test tests/ --profile +docker).docs/output.mdis updated.CHANGELOG.mdis updated.README.mdis updated (including new tool citations and authors/contributors).🤖 Generated with Claude Code
https://claude.ai/code/session_019SiaVrAVHBtxZcGgoiwQun