https://nf-co.re logo
Join Slack
Powered by
I am now getting terminal errors for `NFCORE_RIBOS...
# riboseq
f
I am now getting terminal errors for 
NFCORE_RIBOSEQ:RIBOSEQ:BAM_DEDUP_STATS_SAMTOOLS_UMITOOLS_TRANSCRIPTOME:UMITOOLS_DEDUP:
Copy code
PYTHONHASHSEED=0 umi_tools \
    dedup \
    -I Tm.Aligned.toTranscriptome.out.bam \
    -S Tm.umi_dedup.transcriptome.sorted.bam \
    -L Tm.umi_dedup.transcriptome.sorted.log \
     \
     \
    --method='directional' --random-seed=100
Copy code
Traceback (most recent call last):
  File "/usr/local/bin/umi_tools", line 11, in
    sys.exit(main())
  File "/usr/local/lib/python3.9/site-packages/umi_tools/umi_tools.py", line 61, in main
    module.main(sys.argv)
  File "/usr/local/lib/python3.9/site-packages/umi_tools/dedup.py", line 290, in main
    inreads = infile.fetch()
  File "pysam/libcalignmentfile.pyx", line 1100, in pysam.libcalignmentfile.AlignmentFile.fetch
ValueError: fetch called on bamfile without index
j
OK, maybe the quick fix didn't completely do the trick This PR is about to merge in rnaseq: https://github.com/nf-core/rnaseq/pull/1369 After that my plan is to factor all that UMI handling logic into a subworkflow we can deploy in both workflows. I suggest that we wait and see if the issue goes away when I do that.
f
Sounds like a good plan. Do you have a rough timeline in mind for this? (we are getting new Ribo-seq data in very soon…. 😨 )
j
It's the thing I want to do now. I'm a bit disrupted for the next day or so due to work travel, but should be done within a week or so 🤞
f
🤞
i
@Felix Krueger is this the first time you get this error? was it with just this dataset? Regardless, it seems a pretty clear issue - basically UMI-tools needs the bam file to be indexed, so prolly the bai wasn't staged correctly / the wrong bai is staged. Worthwhile for you to go into the right 
work
directory and see if a bai was staged at all. If there's no bai, not sure what the reason is yet. If there is a bai, there may be a mismatch between the bam and bai, it may be related to these lines in the riboseq workflow:
Copy code
ch_transcriptome_bam       = FASTQ_ALIGN_STAR.out.orig_bam_transcript
ch_transcriptome_bai       = FASTQ_ALIGN_STAR.out.bai_transcript
at the transcriptome level seems like these 2 are mismatched because bai is done only only on the sorted bam, whereas 
ch_transcriptome_bam
is defined as the original transcriptome bam 😛 . I pointed in PR #75 in my comment that I am not sure this is right.. However, when 
ch_transcriptome_bam
is re-defined as the sorted bam, other errors appear.
just my 2 cents on what might be going on, best we wait for @Jonathan Manning to do his magic
f
The following files are getting staged:
so there is no file 
Tm.umi_dedup.transcriptome.sorted.bai
i
Technically umi tools didn’t work, so the 3rd file in the list, the dedup bam should be empty, is it?
Also indeed the 2 first files are mismatched, one is the original alignment, whereas the bai corresponds to the sorted bam, which is why umi tools thinks the bam has no index ;)
f
it is 0 byte indeed
👍 1
i
Oki, it might because of reason I stated above that it fails ( but I have no way to check rn)
f
I think if Jon is about to implement this as a subworkflow for RNA-seq, and then redeploy here it would be the best solution
👍 1
i
Yeah I have zero plans of fixing this although I could; no point in duplicating work:)
💯 1
j
Just so you know I'm working on it 🙂 https://github.com/nf-core/rnaseq/pull/1466 https://github.com/nf-core/rnaseq/pull/1467 The UMI testing functionality in RNAseq wasn't being tested yet, so I needed to add that baseline before I mucked about, and it has proven fiddly. But well on the way now.
f
Thanks Jonathan, that’s much appreciated!
j
UMI handling now factored out to nf-core/modules, and we've tested it carefully in rnaseq. https://github.com/nf-core/rnaseq/pull/1466 Work started on riboseq here: https://github.com/nf-core/riboseq/pull/79
@Felix Krueger @Ira Iosub I still have work to do on documentation and tests, but if you'd like to test the UMI functionality, the 
umi_handling_from_subworkflow
branch on riboseq should work. If you could confirm while I'm tidying up the rest of it, that would be very helpful.
f
It succeeded! That's awesome!
thankyou 1
i
Thx @Jonathan Manning all good with my tests too! Very useful to have UMI testing also for rnaseq given how widely used that workflow is!
👍 1
f
@Jonathan Manning given that you are in our time zone and active right now, I am still having trouble doing stuff with Ribo-seq as it won’t render at all on our platform enterprise. Was there a release earlier this week that we could deploy? I am currently working from a ‘resume’ run from October where I selected the ‘work_with_umi_subworkflow’ branch (which is quite suboptipmal…). Thanks for any update you may have!
j
Hey Felix, let me come back to you later today on the Platform release, we just need to consult on something with engineering. The PR hasn't been merged and released yet because I've had some issues 'faking' UMIs in the test data to add UMI testing. Are you OK working with the branch for now? I can always release without the test stuff for the moment and do a patch release with it later if you'd prefer.
i
Would it be simpler to add new data with UMIs than faking UMIs on the existing data? And have a test_umi profile?
f
I am personally happy either way
j
It was a little bit involved generating the original test data- we have to use some intelligent reduction to a specific region, I was trying to involve re-doing it. We'd also need to select a new test dataset- do you have good candidates? Maybe I can just release as-is, for the moment. Someone else may be able to help out on the tests as well.
I'm off for 2 weeks after this Friday, so I'm a little strapped for time right now.
f
Just a comment on the side which may or may not be relevant for Ribo-seq folks: I have added a new option to Trim Galore to add demultiplexing for in-line barcodes of Ribo-seq libraries that have a strategy like this:
Copy code
ribosome protected sequence - UUUUUUUU - BBBBBB - AGATC.....................
U = UMI
B = barcode
AGATC... = adapter
It’s currently still in the dev branch, but it worked like a charm (see more here: https://github.com/FelixKrueger/TrimGalore/issues/199)
i
The riboseq-flow test data has UMIs let me get back to you so I double check where it’s from
I will also be away so I can’t really add the new test data myself soon
j
OK, I'll at least tidy up the docs and try and release tomorrow, we can worry about the testing next year
👍 1
keanu thanks 1
The UMI subworkflow is internally tested anyway, so I'm not too worried
f
is there any good report from UMI-tools that will yell us what it did?
or is that disabled for time/compute reasons?
j
There are logs you should see, beyond that I'm not sure- I just fixed up the existing logic across the two pipelines
f
There doesn’t seem to be a section for it in the MultiQC report
i
UMI tools doesn’t work with multiQC .. I forgot the reason why.. it’s a known issue or it was the last time I checked
f
so one would have to run a workflow twice, once with and once without UMI tools to see if there is a difference?
i
Logs should be informative
f
ok will try to find some
Update: I can find a UMI-extract.log, but this is largely uninformative (summarising resolved options, and outputting this:
Copy code
...
2024-12-19 08:26:58,359 INFO Parsed 22100000 reads
2024-12-19 08:26:59,389 INFO Parsed 22200000 reads
2024-12-19 08:27:00,356 INFO Input Reads: 22293894
2024-12-19 08:27:00,356 INFO Reads output: 22293894
# job finished in 238 seconds at Thu Dec 19 08:27:00 2024 -- 238.48  2.18  0.00  0.00 -- a50088f0-10f4-43c7-9486-a72d4309b677
There don’t seem to be any logs regarding the deduplication step (which is arguably the one that matters), or am I looking in the wrong place?
i
hold on i'll have s look
j
let me check where they should be published
f
thanks!
i
Copy code
if (params.with_umi && ['star'].contains(params.aligner)) {
        process {
            withName: '.*:BAM_DEDUP_STATS_SAMTOOLS_UMITOOLS_GENOME:UMITOOLS_DEDUP' {
                ext.args   = { [
                    meta.single_end                 ? '' : '--unpaired-reads=discard --chimeric-pairs=discard',
                    params.umitools_grouping_method ? "--method='${params.umitools_grouping_method}'" : '',
                    params.umitools_umi_separator   ? "--umi-separator='${params.umitools_umi_separator}'" : ''
                ].join(' ').trim() }
                ext.prefix = { "${meta.id}.umi_dedup.sorted" }
                publishDir = [
                    [
                        path: { "${params.outdir}/alignment/${params.aligner}/umitools" },
                        mode: params.publish_dir_mode,
                        pattern: '*.tsv'
                    ],
                    [
                        path: { "${params.outdir}/alignment/${params.aligner}/deduplicated" },
                        mode: params.publish_dir_mode,
                        pattern: '*.bam',
                        saveAs: { params.save_align_intermeds || params.with_umi || params.save_umi_intermeds ? it : null }
                    ]
                ]
            }
so the module itself is emitting the logs but it doesnt look like they are published
this can be fixed in this config for the future; but meanwhile you can look into the work dir where UMI dedup was ran, the log should be there
f
true. There adding an entry for 
pattern: '*.log'
should to the trick? I assume this is also not working for the RNA-seq workflow then?
Copy code
2024-12-19 09:45:44,321 INFO Parsed 122000000 input reads
2024-12-19 09:46:05,069 INFO Reads: Input Reads: 122615517
2024-12-19 09:46:05,069 INFO Number of reads out: 5538104
2024-12-19 09:46:05,069 INFO Total number of positions deduplicated: 1602398
2024-12-19 09:46:05,069 INFO Mean number of unique UMIs per position: 5.09
2024-12-19 09:46:05,069 INFO Max. number of unique UMIs per position: 1649
This is from the file:
🎉 1
j
Let me double-check, I'm sure they were getting published, because I was having to exclude them from workflow-level tests
🙏 1
i
what i see is this:
Copy code
if (params.with_umi) {
        process {
            withName: 'NFCORE_RIBOSEQ:RIBOSEQ:SAMTOOLS_SORT' {
                ext.args   = '-n'
                ext.prefix = { "${meta.id}.umi_dedup.transcriptome" }
                publishDir = [
                    path: { "${params.outdir}/${params.aligner}" },
                    mode: params.publish_dir_mode,
                    pattern: '*.bam',
                    saveAs: { params.save_align_intermeds || params.save_umi_intermeds ? it : null }
                ]
            }

            withName: 'NFCORE_RIBOSEQ:RIBOSEQ:UMITOOLS_PREPAREFORSALMON' {
                ext.prefix = { "${meta.id}.umi_dedup.transcriptome.filtered" }
                publishDir = [
                    [
                        path: { "${params.outdir}/${params.aligner}/umitools/log" },
                        mode: params.publish_dir_mode,
                        pattern: '*.log'
                    ],
                    [
                        path: { "${params.outdir}/${params.aligner}" },
                        mode: params.publish_dir_mode,
                        pattern: '*.bam',
                        saveAs: { params.save_align_intermeds || params.save_umi_intermeds ? it : null }
                    ]
                ]
            }
j
Yes, but there's usually default publishing rules
i
ok so it wasnt being published, tx for the change!
j
Yeah, I was testing with umicollapse, whose logs were being published 😉
😅 1
f
excellent. Is this limited to Ribo-seq then, or also an issue for RNA-seq?
j
Good point- it'll be in rnaseq too, since I copied the configs. I was just about to release there, will fix before I do
💯 1
4 Views
Open in Slack
PreviousNext
Powered by