Hey rna-seq developers! A few months ago I published a GPU accelerated STAR module on nf-core built on Parabricks rnafq2bam. Now I want to add it to the rnaseq pipeline to make it easier to deploy (See Github issue). However, I'm having trouble figuring out the best way to insert this module. Rnaseq right now runs ALIGN_STAR then BAM_SORT_STATS_SAMTOOLS then BAM_SORT_STATS_SAMTOOLS, however Parabricks runs STAR and MarkDuplicates back to back with no way to split it up. So I would have to create a new path for this workflow, running BAM_SORT_STATS_SAMTOOLS afterwards. My question is: 1. Is this new path okay to create? Is there a better way I should handle it? 2. Will this eventually require a change to the subway map?

For me yes to both. And I think Jon expects me to update the map anyway

Thank you for your quick response! I will get to work on this change. If anyone else has any thoughts on this or other GPU integrations into rnaseq, please let me know!

Maybe wait for @Jonathan Manning. I am overly enthusiastic about pretty much everything

Yeah, my concern is that we already have the accelerated STAR @Rike Hanssen added via sentieon, I'm concerned things will get super confusing for users

I see your concern about potential user confusion. I would have to make it clear when users are selecting CPU acceleration (Sentieon) or GPU acceleration (Parabricks).

I think it should be as simple as adding an option to the align_star subworkflow in the same way as the sentieon stuff, and then use a conditional to skip the mark duplicates step later on (so modify this line). The tube map is by necessity a simplification anyway, I don't think this really needs a change there. But I'm allergic to them anyway, so it would be Maxime!

Thank you, Jonathan. I understand what you're saying and that would be a minimally invasive change.

Maxime bought the perpetual responsibility for the tube map when he animated it!

Fair, I wouldn't touch that either

One small remark, though regarding this:

however Parabricks runs STAR and MarkDuplicates back to back with no way to split it up.

In that case, you should probably add a conditional samtools run to it to remove those duplicate markings again in case

with_umi

is true.

Thank you for the note, Matthias. There is functionality to skip the Parabricks MarkDuplicates step altogether and just run STAR. Would it accomplish the same goal if I tell Parabricks to skip MarkDuplicates when

with_umi==true

?

Yes, exactly. In an NGS experiment, duplications may arise for biological and technical reasons. For Illumina instruments and chemistries, five different sources are of relevance: • True biological duplication: The extent of true biological duplication varies greatly depending on the experiment and sequencing depth. Due to the myriad of possible combinations, biological duplication in Hi-C experiments is entirely absent, whereas similar fragments in cDNA libraries used for RNA-seq may well originate from highly expressed genes. In WGS, low-level duplication may be due to repeats, yet virtually no sequences should occur more than once in any lane in PCR-free WGS libraries. • PCR and enrichment duplicates: Even for PCR-amplified whole genome/exome sequencing libraries, most sequences will only appear once in the lane, if PCR amplification was started with a high enough input. Only when starting with a very low input, duplicate reads in a sequencing lane are much more likely and may be a problem in a WGS/WES library. In contrast, PCR duplicates are particularly relevant for amplicon sequencing methods like 16S/18S or STARR-seq. The fewer unique regions are present, the more random biases in a PCR will affect biological interpretation by skewed over-amplification of some regions. • Optical duplicates in Illumina data may occur on HiSeq 2500/MiSeq/NextSeq machines. They emerge from large clusters being called as two separate clusters, or from local re-clustering of the original library molecule once it is released from polymerase copying. Since duplicate reads are characterized by close proximity to each other on the same tile, detecting such optical duplicates bioinformatically is straightforward. • Exclusion Amplification duplicates are restricted to patterned flowcells, used e.g. in the HiSeq 4000/X or NovaSeq 6000 systems. Similarly to the non-patterned flowcells, local re-clustering of the original library molecule may occur once it is released from polymerase copying. • Tile-edge duplicates are limited to the border zone of distinct tiles. In particular NextSeq machines suffer from this type of duplicates, where they account for the vast majority of all duplicates (>80%). In a NextSeq, the six CCD cameras record the nine tiles sequentially. Poor calibration and motor imprecision can result clusters close to the border being recorded twice in adjacent tiles. In RNA-seq, duplicated cDNA fragments may reflect true biological duplication arising form highly expressed genes. Therefore, physically removing duplicate reads from the alignments is discouraged, unless the same UMI in the adapter confirms technical duplication. By default, and without having UMI information available, the RNA-seq pipeline therefore runs picard’s

MarkDuplicates

, to mark but not remove duplicates purely based on sequence similarity of the fragment alone. Unique Molecular Identifiers are short, artificial sequences used to create unique adapter-insert combinations in a sample library, which facilitates the accurate identification of read duplicates. Like in a HiC library, the the myriad of possible combinations ensures that the combination of fragment and UMIs is surely unique in the library. If the same fragment is encountered multiple times with the identical UMI, it can be discarded as a certain duplicate. While MarkDuplicates assumes a fragment may be a possible duplicate, umi-tools dedup or umi-collapse can tell for sure. Therefore, running those tools instead of MarkDuplicates is preferable, if UMIs are available. This will generate a filtered BAM file after the physical removal of duplicates. The

--umitools_grouping_method

parameter affects how similar, but non-identical UMIs are treated.

directional

, the default setting, is most accurate, but computationally demanding.

percentile

or

unique

can be considered if many samples are to be processed.

Thank you, Matthias for explaining it to me! I will make these changes to my branch and submit the PR soon. Will post here when it's ready for anyone who wants to review.

I see it here @Gary Burnett: https://github.com/nf-core/rnaseq/pull/1685 Are you done?

I'm actually not done yet. That was a premature PR. But I can ping you when it is ready!

Sorry to break it to you, but to be it seems that also the latest-everything tests do actually not work. Because you can see lots of errors littered in the passed tests like:

ERROR ~ Error executing process > 'DESEQ2_QC'

Caused by:
Process requirement exceeds available CPUs -- req: 6; avail: 4

ERROR ~ Error executing process > 'STAR_GENOMEGENERATE_IGENOMES (genome.fasta)' 
 Caused by: 
 Process requirement exceeds available CPUs -- req: 8; avail: 4

Depending on whether these modules can run with fewer cores, you would either need to adapt your test configs or talk to the <!subteam^S08LAB76NBU> about different self-hosted runners for tools like this. In the failing tests, there is additionally this error, which may be related to the failing subtasks not publishing the expected output? But that is just an assumption, since I am myself far from a Nextflow expert…

> ERROR ~ Workflow `PREPARE_GENOME` declares 27 input channels but 26 were given 
 > 
 > -- Check script '/home/runner/_work/rnaseq/rnaseq/.nf-test-750785c194709d4a30dbb88c8b3d614c.nf' at line: 72 or see '/home/runner/_work/rnaseq/rnaseq/~/tests/750785c194709d4a30dbb88c8b3d614c/meta/nextflow.log' file for more details 
Assertion failed: 
 
assert workflow.success 
 | | 
 workflow false 
 
java.lang.RuntimeException: Workflow has no output channels. workflow.out can not be used. 
 FAILED (9.544s)

oh yep, the latest tests run for information, but always pass so as not to block PRs

Thank you for letting me know about the errors in those latest-everything tests. Parabricks will require slightly more resources to run the test data so I will start a conversation with the infra team about runners. I ran into this issue when I initially submitted this module for review and would love to find a more permanent solution.

Would scaling test data down help?

Let me run some more tests and get back to you. I will need to submit another PR to update that module anyways (had to change the outputs to better integrate into the rnaseq pipeline).

We can totally upload some special test data to rnaseq test-data branch if useful

I found some unexpected (to me) behavior in the

prepare_genome

tests. When I run the Default Options test in Prepare Genome (using

nf-test test --profile debug,test,docker,gpu --verbose --update-snapshot --tag "default options" subworkflows/local/prepare_genome/tests/main.nf.test

) on the latest

dev

branch, the subworkflow stops at

BBMAP_BBSPLIT

and then says the test finished successfully. But when I look at subworkflow code, there are several steps I believe it should be running after that (Ex. STAR GenomeGenerate). This is preventing me from testing the Parabricks star genome_generate, but maybe I am missing a parameter to ensure it keeps going.

You're right, doesn't look like that test file was ever written correctly. I'll PR a fix

Thank you for the fast response! I will review this change now.

Looks like I need to work on the tests, will probably finish it next week, but you see the way it's goign

Yeah this is helpful to get a peek into what the changes will be. Thank you.

sorted! just needs approval

That's awesome! I'm looking through right now.

Hey @Gary Burnett can you link me to one of the errors?

Sure thing, here is one of them

For GPU runners, can you add a GPU CI file like this https://github.com/nf-core/methylseq/blob/dev/.github/workflows/nf-test-gpu.yml. Be sure To also tag the tests with cpu and gpu so they get picked up correctly

Thank you @mashehu, I have the GPU runners up and running. The methylseq example you shared was extremely helpful. My new GPU module is passing the the subworkflow CI tests (yay I can almost taste victory), but hitting SIGKILL errors in the full pipeline tests, no matter how much I limit Parabricks memory usage (using the

--memory-limit

flag). It leads me to believe that the full pipeline is running more concurrent tasks. I've played around with

maxForks

to try to limit that, but I'm not sure if that's the right idea or if this is a common occurrence when scaling up to the full pipeline tests.

@Gary Burnett I'm going to merge a couple of big PRs I've been working on this afternoon that will generate conflicts for you. But I can resolve them for you.

Thank you for the heads up!

For now we can just document it as non-ARM. The longer-term fix is to get an arm Conda build done. We already have an old conda build and matching container (2.6.1d, via Wave) to support the ancient STAR indices in iGenomes, but it's in the

seqera

channel rather than

bioconda

, which isn't ideal.

Thank you, Jonathan. I would like to do that in my next PR and I found some some nf-core docs that I can get started with. Then all I'm left to wrestle with is this sample in the samplesheet that's running OOM in the full pipeline test. I can't think of a way to get past this besides: 1. Creating a new samplesheet without this sample that I would only use for the Parabricks tests 2. Increasing the memory resources limit in

conf/test.conf

from 15GB to 16GB which works locally and fits on the ci runner instance, but I'm not sure if removing that extra GB will cause other problems Should I take either of these approaches? I'm open to other ideas as well.

One option might be to use the prokaryotic test sheet I just added for https://github.com/nf-core/rnaseq/pull/1688? I bet those would reduce the memory requirements. Sorry I didn't resolve your conflicts today, other things came up, but can do tomorrow

Oooh I'll take a look at this new test sheet. And don't worry about not getting to the conflicts, I have plenty to keep me busy until then. I can take a look at it tomorrow morning, as well.

Long short term hoping to just drop this

@Gary Burnett I have: 1. switched the test strategy to the one I was describing that didn't need all the

cpu

tags etc 2. reconciled conflicts with dev 3. made a PR against your PR with some suggestions on better integration with existing STAR configs: https://github.com/gburnett-nvidia/rnaseq/pull/1 That includes use of the prokaryotic data for the pipeline tests, which still OOMs, so hasn't helped. I may revert that shortly.

I talked to the core Parabricks engineering team about this OOM issue and they confirmed that it's the STAR index directory that is causing such a large memory footprint even when the input fastq are small. I did some local testing and was able to get it running by bumping the memory requirements from 15GB to 30GB, but this would mean that the CI runner would have to get bumped from a

g4dn.xlarge

to a

g4dn.2xlarge

. We did this when I originally published the rna_fq2bam module as a one-time solution. Is it possible to do it again here? This is more of an infrastructure question, I'm happy to bring to another channel if it's more appropriate somewhere else.

If it is needed to make it run, go for it. Is star used in all gpu tests?

STAR index directory that is causing such a large memory footprint even when the input fastq are small.

Unless I am missing something the size of the index should be determined by the reference genome and not by the fastq data that is aligned. Switching from the 3.3 Billion bp human genome to a smaller reference genome would probably mitigate the OOM issue on the runner?

test("parabricks star") {
        setup {
            run("PARABRICKS_STARGENOMEGENERATE") {
                script "../../../../modules/nf-core/parabricks/stargenomegenerate/main.nf"
                process {
                    """
                    input[0] = Channel.of([
                        [ id:'test_fasta' ],
                        [ file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true) ]
                    ])
                    input[1] = Channel.of([
                        [ id:'test_gtf' ],
                        [ file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.gtf', checkIfExists: true) ]
                    ])
                    """
                }
            }
        }

I actually don't think that this is the issue, that file IS a small subset of chr22 IIRC. fq2bam itself seems to have some high baseline CPU mem requirement.

I actually don’t think that this is the issue, that file IS a small subset of chr22 IIRC. Fair. Saturday night at 1 a.m. I didn’t bother to check anymore. I just skimmed over the PR and noticed

homo_sapiens/genome/genome.fasta

. I assumed that the file name would somehow reflect if it was downsampled.

No worries! @Gary Burnett, when you're here, I think I mostly got it with the new machine type as above, after some tinkering about on how config is applied so that we don't get multiple fq2bam runs happening together.

Hey Jonathan! I'm just starting my day over here and catching up on the Slack threads. Thank you so much for helping me push this over the edge. I'm surprised the

--low-memory

flag didn't fix that T4 16 GB issue. But if you have something that's working, I'll go take a look now.

Yep, think it's good, and working, you way want to refine if I've done something silly. The only remaining failure was a snapshot issue because the Log.out etc weren't being prefixed to match the snapshot exclusion patterns in the pipeline. I've just fixed that in the pipeline to test it, but you'll need to update the module in nf-core/modules

I will go do that now and get back to you shortly!

Could you also do the topics migration please while you're at it, for the version handling?

Yes! I'll follow the instructions I found here

Thanks! The rest of the workflow is mostly converted, so you can check there for examples too

Huge thanks to @Jonathan Manning @mashehu and everyone who helped on this thread. The PR is now closed and merged into dev 🎉

Well done Gary, turned out to be a somewhat thorny one because of the other STAR behaviours it touched on

It was definitely a delicate operation, but I learned a lot about rnaseq and nf-core dev in general. Looking forward to doing more GPU integrations across the platform.