Add RNA-seq Variant Discovery (RNAvar) workflow

[emregulben]

Summary

Adds the RNA-seq Variant Discovery (RNAvar) workflow.

Technical Details

Implementation based on nf-core/rnavar logic:

• Alignment: STAR (2-pass mode).
• Preprocessing: GATK4 MarkDuplicates, SplitNCigarReads, and BaseRecalibrator.
• Variant Calling: GATK4 HaplotypeCaller using RNA-seq parameters.
• Annotation: SnpEff (hg38).

Validation

• Tests: Planemo-validated suite verifying workflow outputs and technical parity with reference results.

FOR CONTRIBUTOR:

• I have read the Adding workflows guidelines
• License permits unrestricted use (educational + commercial)
• Please also take note of the reviewer guidelines below to facilitate a smooth review process.

FOR REVIEWERS:

• .dockstore.yml: file is present and aligned with creator metadata in workflow. ORCID identifiers are strongly encouraged in creator metadata. The .dockstore.yml file is required to run tests
• Workflow is sufficiently generic to be used with lab data and does not hardcode sample names, reference data and can be run without reading an accompanying tutorial.
• In workflow: annotation field contains short description of what the workflow does. Should start with This workflow does/runs/performs … xyz … to generate/analyze/etc …
• In workflow: workflow inputs and outputs have human readable names (spaces are fine, no underscore, dash only where spelling dictates it), no abbreviation unless it is generally understood. Altering input or output labels requires adjusting these labels in the the workflow-tests.yml file as well
• In workflow: name field should be human readable (spaces are fine, no underscore, dash only where spelling dictates it), no abbreviation unless generally understood
• Workflow folder: prefer dash (-) over underscore (_), prefer all lowercase. Folder becomes repository in iwc-workflows organization and is included in TRS id
• Readme explains what workflow does, what are valid inputs and what outputs users can expect. If a tutorial or other resources exist they can be linked. If a similar workflow exists in IWC readme should explain differences with existing workflow and when one might prefer one workflow over another
• Changelog contains appropriate entries
• Large files (> 100 KB) are uploaded to zenodo and location urls are used in test file

[Copilot]

Pull request overview

This PR adds a new Galaxy workflow under workflows/variant-calling/rna-variant-discovery/ implementing an RNA-seq variant discovery pipeline (STAR → GATK preprocessing/calling → filtering → SnpEff annotation) along with Dockstore metadata, documentation, and a Planemo test suite with bundled test data.

Changes:

• Added the RNAvar.ga Galaxy workflow for RNA-seq variant calling and annotation, with MultiQC reporting.
• Added a Planemo test definition (RNAvar-tests.yml) and corresponding small test datasets.
• Added workflow packaging/documentation files (.dockstore.yml, README.md, CHANGELOG.md).

Reviewed changes

Copilot reviewed 9 out of 11 changed files in this pull request and generated 3 comments.

Show a summary per file

File	Description
workflows/variant-calling/rna-variant-discovery/RNAvar.ga	New Galaxy workflow definition for RNA-seq variant discovery and reporting
workflows/variant-calling/rna-variant-discovery/RNAvar-tests.yml	Planemo test job + output assertions for key workflow outputs
workflows/variant-calling/rna-variant-discovery/.dockstore.yml	Dockstore registration metadata pointing to the workflow and tests
workflows/variant-calling/rna-variant-discovery/README.md	High-level explanation of workflow steps and expected behavior
workflows/variant-calling/rna-variant-discovery/CHANGELOG.md	Initial release entry
workflows/variant-calling/rna-variant-discovery/test-data/test_rnaseq_1.fastq.gz	Test input R1 reads for Planemo validation
workflows/variant-calling/rna-variant-discovery/test-data/test_rnaseq_2.fastq.gz	Test input R2 reads for Planemo validation
workflows/variant-calling/rna-variant-discovery/test-data/genome.fasta	Small reference FASTA used by tests
workflows/variant-calling/rna-variant-discovery/test-data/genome.gff3	Small annotation file used by tests
workflows/variant-calling/rna-variant-discovery/test-data/dbsnp_146.hg38.vcf	Known variants VCF used for BQSR / known-sites
workflows/variant-calling/rna-variant-discovery/test-data/mills_and_1000G.indels.vcf	Known indels VCF used for BQSR / known-sites

[Copilot]

SnpEff is configured with a hardcoded named genome database (hg38). Since the workflow also accepts a user-provided reference FASTA and annotation, this can silently produce incorrect annotations (or fail) when the inputs are not hg38. Consider either (a) building a SnpEff database from the provided reference+annotation (as other variant-calling workflows here do) or (b) exposing the SnpEff genome version as an explicit workflow parameter and documenting that the workflow is hg38-specific.

[lldelisle]

If I may add to copilot review, I think would be good to describe in the README each of the expected inputs and the expected outputs.
Also add the README in the workflow so users can display it as 'help for this workflow'.
I am not familiar with the field, do you usually provide a fasta or do you use built-in genome indexes like hg38 etc..?

[wm75]

It makes limited sense to hard-code the SnpEff database and, thereby, restrict analysis to just one organism and genome version. The database choice should instead be turned into a workflow parameter.

I would populate the user choices from the "Locally installed SnpEff databases" instead of letting the user download databases on demand.

Please let us know if you need help with this.

[emregulben]

Done here: 5146320

[wm75]

Here, the README should clearly state that the aim is to provide a faithful reimplementation of https://nf-co.re/rnavar and that currently there are still limitations.

[wm75]

After presenting the current workflow logic, it wouldbe good to have a "Limitations" section, where you discuss which features of the nf-core ppeline are currently not yet implemented.

Things hat I can immediately think of here are

• Sample Sheet support
• support for SnpEff alternatives for variant annotation
• optional UMI-Tools-based preprocessing
• no parallelization of the GATK4 SplitNCigarReads step (is this something that should be implemented at the tool level similar to what's done in the freebayes wrapper?
• base recalibration and variant filtering are currently non-optional when they are optional in th nf-core pipeline

Basically, compare https://nf-co.re/rnavar/1.2.3/docs/usage/ to what's implemented here.
Of course, if you still want to fix any of these limitations, you're welcome to do so.

[emregulben]

Added the limitations section: 5146320

[wm75]

Be more precise here and state the origin of that test data.

[wm75]

The RNAStar Galaxy wrapper should already generate coordinate-sorted BAM output, so is this step necessary?

Please check and/or clarify.

[emregulben]

I verified STAR already handles the sorting. Removed the step: 5146320

[wm75]

Try to understand if and how this step changes the data. Is there a way to use it in a meaningful way either now without having implemented sample sheet funtionality or maybe in the future?

[emregulben]

This step is required for the GATK4 tools used later in the workflow. Tools like HaplotypeCaller and BaseRecalibrator are designed to look for 'Read Group' metadata in the BAM header. Without these tags, the tools throw an error and won't process the file (as noted in the GATK documentation on Read Groups).

By including this step now with default values, we ensure the workflow is robust and GATK-compliant. It also sets us up nicely for the future. If we implement sample sheet support later, this tool is exactly where that metadata would be injected into the pipeline.

[wm75]

Please order the different MultiQC input reports in a logical way proceeding from raw read metrics through alignment, variant calling and annotation.

[emregulben]

5146320 ✔️

[wm75]

The primary advantage of building the STAR index on te fly (which otherwise is just very wasteful) is that you can control the --sjdbOverhang parameter (see https://nf-co.re/rnavar/1.2.3/docs/usage/#alignment-options).
This means that you need to expose the sjdbOverhang wrapper parameter as a workflow parameter with a helpful annotation.

[lldelisle]

For me sjdbOverhang is always exposed if you use a prebuilt genome reference without a prebuit gene model.

[wm75]

That's also a good reason, yes.
Do you think we should go this way @lldelisle ? Let the user choose from cached, prebuilt references, but bring their own gene model? Could be a very good compromise.

[lldelisle]

This is always the way I do, this is also the way it is done in the RNA-seq IWC workflow.
I think it should be cached prebuilt reference + own gtf.
Therefore I would add a step toolshed.g2.bx.psu.edu/repos/iuc/snpeff/snpEff_build_gb/5.4+galaxy0 to build the snpEff database matching the same reference fasta and the same gtf.
Using the cache option (reuse jobs with identical parameters) will allow to not really rebuild the database if it was already built with the same gtf.

[wm75]

This is always the way I do, this is also the way it is done in the RNA-seq IWC workflow.
I think it should be cached prebuilt reference + own gtf.

Cool, sounds excellent to me, thanks!

In that case, however, how can we make sure that the reference fasta matches the one used for building the STAR index?

Also,

Therefore I would add a step toolshed.g2.bx.psu.edu/repos/iuc/snpeff/snpEff_build_gb/5.4+galaxy0 to build the snpEff database matching the same reference fasta and the same gtf.

I'm not so convinced this is a good idea, building SnpEff annotation files from gff is a fragile art for eukaryotic genomes and the pre-built cached ones might contain fewer errors.

So, essentially, the question to answer is: what's the least error-prone way to have matching genomes at the STAR, GATK and SnpEff steps.

Potential differences between SnpEff and custom gff annotations are another thing, but wouldn't be catastrophic for the workflow. Users who have a carefully crafted SnpEff db based specifically for their gff could also ask for server-side installation?
For exotic cases where there's no reasonably close existing SnpEff db, we could offer an optional branch of the WF to build the db from the reference and the gff, I just wouldn't make this the default.

Of course, I'm only making all of this up while typing so feel free to disagree with me :-)

[emregulben]

I’ve exposed sjdbOverhang as a mandatory workflow parameter. Regarding the SnpEff database discussion, I’ve kept the database selection as a string input to use pre-built/cached local versions, avoiding the potential fragility of building custom databases from GFFs during the run. Let me know if the current version is fine. 5146320

[wm75]

These filters are user-configurable in the nf-core pipeline so should become workflow parameters, too.

[emregulben]

Done: 5146320 ✅

[lldelisle]

You need to force to choose among the prebuit SnpEff databases:

With attempt restriction based on connections

[lldelisle]

Here is my proposed variation: https://usegalaxy.eu/u/delislel/w/rna-seq-variant-discovery-ld-version
What I changed:

• use as input fastqs: PE collection, I think it is easier if you want to process multiple samples at the time and even if you don't this ensure that all your output files will have the name of your sample.
• use built-in genome instead of fasta
• add an output to STAR (genecounts) for me this is very useful to see this statistics in the MultiQC to check you put the right GTF (you should have a high proportion of reads mapping to genes if the GTF is correct).
• I also rearranged the inputs to they are in the order they are used in the pipeline
• Finally, I put 'hidden' all intermediate steps

Tell me what you think and feel free to change it.