feat: Reimplementation of umi_tools prepare-for-rsem

[MatthiasZepper]

Description

Unless I have overlooked something, fgumi currently has no equivalent tool to the umi_tools prepare-for-rsem command.

Since this command frequently takes a significant time e.g. in the nf-core rnaseq pipeline, a faster rewrite in Rust would be desirable. Since fgumi will prospectively be supplementing/replacing umi_tools in that pipeline, adding this functionality to fgumi rather than publishing a standalone tool seems reasonable.

I have an AI-generated draft for that feature sitting on a branch of my fork as preview. I have however not yet reviewed it extensively and am currently running comparisons of both tools with real-world data.

Please let me know if you are generally open to that addition, then I will spend more time and effort on human curation to get it into shape for a proper PR.

Related Issues

No response

Version

0.2.0

Impact

Inside the rnaseq pipeline, prepare-for-salmon/rsem is frequently the longest running task, taking up to 12h per sample. Far longer than alignment and also longer than umi_tools extract and umi_tools dedup. Replacing this script could thus significantly shorten runtimes of the pipeline and it would be great if fgumi would cater for that step as well.

Suggestions

AI-generated draft for that feature

Test input

I am still looking for suitable test data. Neither umi_tools itself nor the nf-core module seem to be using a particular test dataset that truly covers all edge cases.

Checklist

• Read the CONTRIBUTING guide
• Checked there are no existing Issues for this feature request
• Checked there are no Discussions that address this issue
• Checked you have installed at least the minimum required version of Rust
• Ensured that the latest release does not contain this feature
• Provided all information necessary to reproduce the current behavior

[nh13]

Thanks for raising this, @MatthiasZepper.

Before scoping a new tool, I want to push on whether fgumi already covers this. A few things worth checking:

1. RSEM's actual ordering requirement is query-grouped, not query-name sorted. From the rsem-calculate-expression docs: "RSEM requires all alignments of the same read group together. For paired-end reads, RSEM also requires the two mates of any alignment be adjacent." --sort-bam-by-read-name is an output option for determinism, not an input requirement.

2. fgumi dedup output already satisfies that. It is template-coordinate sorted (all records of a template adjacent) and preserves secondary/supplementary alignments via the tc-tag mechanism added by fgumi zipper. So fgumi zipper → fgumi sort --order template-coordinate → fgumi dedup → rsem-calculate-expression should work directly.

3. The umi_tools tool exists to repair unpaired secondaries (cases where an R1 alignment has no corresponding R2 alignment record), not to enforce "one R1 and one R2" per read. RSEM is fine with multi-mappers as long as each R1 alignment has its mate R2 alignment adjacent.

Could you check, on your real-data comparison, whether the post-fgumi dedup BAM actually has unpaired secondaries that rsem-calculate-expression rejects? If it runs cleanly, the right action here is documentation rather than a new tool. If it doesn't, knowing exactly which records are unpaired (and which aligner produced them) would let us scope a minimal fix.

[MatthiasZepper]

Thank you very much for your instant and thoughtful reply. I realize now this issue requires more attention and careful exploration. So far, I was simply annoyed about the long run times and pointed the AI agents to the existing Python script for porting it to Rust using the libraries and code style of fgumi.

I am neither familiar with the intricacies of the SAM format specification nor know which edge case prompted the initial creation of that script. But I can certainly dig out some info from old nf-core Slack posts and rnaseq pipeline issues. Allow me a few days to think that through and ideally identify/generate a suitable test dataset.

[nh13]

@MatthiasZepper if you have a dataset I can test on, or point me to the testdata on nf-core or elsewhere, I'd be glad to also do my own verification. I'd be glad to help, as I am now also on the annoyed bandwagon of it being slow!

BTW, I really really appreciate how you framed the original issue and the transparency about the state of your implementation, I am going to learn a few things from that.

[nh13]

@MatthiasZepper feel free to reopen if you have more questions!

[MatthiasZepper]

Sorry, I should have kept you up to date. The project in unfortunately on hold at the moment.

I discussed that performance bottleneck with my manager and for the time being, an internal web application that I am building has the higher priority, since I need input from colleagues for that one. You know that after Midsommar many people are on holiday here in Sweden, so that would complicate matters. But I will be working over the summer and we figured that would be the more appropriate time to revisit that issue. I was thinking to build a derivative of the RNA-seq pipeline to benchmark all three solutions side-by-side including quantification.

So please go ahead with the stable release of your tool as planned. Your explanation why it is not really needed sounded quite convincing.