Hi all, we were just reviewing a recent run of the workflow with Bowtie2-based rRNA removal (which worked well). For a test sample, 46 million sequences survived the adapter/quality filtering and rRNA-filtering, and were fed into STAR. Around 36M sequences aligned either uniquely or multiply. The UMI-deduplication using UMI-tools reports however: 12M alignments left over 82M alignments removed This is obviously a lot more alignments that had been fed in, so I assume it also uses secondary alignments. At least in my mind, these stats are not really what we would like to have reported, as we would like to go from

raw reads

to

trimmed_reads

to

non_rRNA_reads

to

aligned_reads

to

deduplicated_reads

, as a QC metric for the experiment/library/CRO. Could we instruct UMI-tools to only consider primary alignments, or filter to primary alignments before (

samtools view -F 0x900)

? Alternatively, should we limit STAR to output only a single alignment, such as with:

--outFilterMultimapNmax 1

? Any opinions? @Jonathan Manning

Interesting observation. I wonder if this might be related to CGATOxford/UMI-tools#454, which describes UMI-tools dedup duplicating mate records when it encounters secondary alignments. If that's what's happening here, it could explain why the stats show more alignments processed than reads fed in. It's also possible UMI-tools is just counting all alignments (including secondaries) in its stats, which would inflate the numbers even without the duplication bug. Either way, your suggested approaches seem reasonable: 1. ***Filter to primary alignments before dedup*** (

samtools view -F 0x900

) - stats would then reflect actual read counts, though you'd lose multi-mapping info in the deduplicated BAM 2. ***Limit STAR output*** (

--outFilterMultimapNmax 1

) - cleaner, but discards multi-mappers entirely which might not be ideal if you want Salmon to use them for quantification I'd probably lean toward option 1 for the QC metric flow you're describing, but I haven't tested either approach with this specific issue. Worth experimenting?

For the reported stats, we like to report read( pair)s, but UMI tools apparently reports alignments. The relationship between reads to alignments is one-to-many (for mapped reads). So whenever you report something, you have to distill out the actual reads from the alignments. With mate pairs that is even more complex, because each mate may or may not have one or more alignments. The way to count mapped read pairs that is guaranteed to be correct is

samtools view -F 4 "$bam" | cut -f 1 | sort -u | wc -l

. With the sorting this is of course highly inefficient, so ideally there are filters that eliminate the sorting step, like using "primary alignments", but that depends on whether the aligner uses that flag correctly, and you still have to sort out how to get from aligned mates to a number of read pairs (some pairs may only have one mate aligned and the other unmapped). The dirty way is

-f 64

to only count the "first in pair", but that is inaccurate of only the "second in pair" (its mate) is mapped. I would not throw out multimappers at the BAM creation (

--outFilterMultimapNmax 1

), it makes wrong pretenses about the accuracy of the alignment. STAR does by default create a few more multimappers than required, by counting alignments that differ by only one mismatch as the same (e.g. a read with one perfect alignment and one or more other alignments with one mismatch would be counted as multimapped). This can be switched off with

--outFilterMultimapScoreRange 0

, but since it only affects a few reads, it does nothing to keep the statistics straight.

Thanks for your thoughts on this. In the case of Ribo-seq we are sort of lucky as it tends to be single-end only, so a lot less complex than for paired-end reads. How is this situation handled in say the RNA-seq workflow with paired-end reads?

How is this situation handled in say the RNA-seq workflow with paired-end reads?

You'd have to look at the #CE8SSJV3N source code for details, but with RNAseq there are fewer steps to lose reads. I think for the most part it is sufficient to look at the trimming losses, the STAR log (STAR has to handle every single read pair anyway, so it can keep track of where it goes without extra cost) and then at (sum of the) the reads counted (by Salmon or RSEM).

Maybe we could add a flag

--primary_only

for the UMI-deduplication (possible set to ON by default) ...

It would be nice to see in the BAM all the locations a certain read maps to. And if you make it optional, it won't help with the stats, because you can't rely it. Since STAR sets (AFAIK) the "primary alignment" flag correctly, using that for counting should be fine. (I don't have enough experience with UMIs to say if there is any other log output, etc., that might help.)

Coming back to this:

It would be nice to see in the BAM all the locations a certain read maps to.

If we went with filtering for primary alignments only before the deduplication, one could still go back to the original BAM file to look at all different mapping locations if one wanted to - so that would be fine?

And if you make it optional, it won’t help with the stats, because you can’t rely it.

I would suggest that the default is to set this to

true

by default so the stats check out. Switching it off would then have to come with a warning that says that including secondary alignments in the UMI-deduplication step would lead to stats that are somewhat disconnected to the number of input reads (“proceed at your own risk”).

I still think it would the cleanest approach to just disregard non-primary reads when counting (for the statistics), rather than removing them from the BAM. What would be the problem/disadvantage with this?

I suppose the problem with this is that UMI tools doesn’t seem to support this… This is Claude’s wisdom: ” no, UMI-tools doesn’t have a “primary only” option; you get the correct reporting by pre-filtering to primary (and optionally non-supplementary) alignments with samtools view -F 256 -F 2048 and then running UMI-tools on that BAM (or stream).”

I was thinking of using

samtools view -c

(with reasonable filters) to do the counting, independently of UMI tools. This is not nearly as efficient as taking values from a log file, but if the log doesn't show what we want, it's at least correct. There are no other outputs from UMI tools from which the correct values can be derived, right?

This is a UMI-tools log file, which is being picked up by MultiQC:

# UMI-tools version: 1.1.5
# output generated by dedup -I BulkD83_S8.genome.sorted.bam -S BulkD83_S8.umi_dedup.genome.sorted.bam -L BulkD83_S8.umi_dedup.genome.sorted.log --output-
stats BulkD83_S8.umi_dedup.genome.sorted --method=directional --random-seed=100
# job started at Thu Jan 22 11:20:19 2026 on <http://ip-10-3-1-113.lab.altoslabs.com|ip-10-3-1-113.lab.altoslabs.com> -- 8be2ef46-f571-4fec-a4d8-b5e157b0049c
# pid: 132, system: Linux 6.1.141-155.222.amzn2023.x86_64 #1 SMP PREEMPT_DYNAMIC Tue Jun 17 10:29:47 UTC 2025 x86_64
# assigned_tag                            : None
# cell_tag                                : None
# cell_tag_delim                          : None
# cell_tag_split                          : -
# chimeric_pairs                          : use
# chrom                                   : None
# compresslevel                           : 6
# detection_method                        : None
# filter_umi                              : None
# gene_tag                                : None
# gene_transcript_map                     : None
# get_umi_method                          : read_id
# ignore_tlen                             : False
# ignore_umi                              : False
# in_sam                                  : False
# log2stderr                              : False
# loglevel                                : 1
# mapping_quality                         : 0
# method                                  : directional
# no_sort_output                          : False
# out_sam                                 : False
# output_unmapped                         : False
# paired                                  : False
# per_cell                                : False
# per_contig                              : False
# per_gene                                : False
# random_seed                             : 100
# read_length                             : False
# short_help                              : None
# skip_regex                              : ^(__|Unassigned)
# soft_clip_threshold                     : 4
# spliced                                 : False
# stats                                   : BulkD83_S8.umi_dedup.genome.sorted
# stderr                                  : <_io.TextIOWrapper name='<stderr>' mode='w' encoding='utf-8'>
# stdin                                   : <_io.TextIOWrapper name='BulkD83_S8.genome.sorted.bam' mode='r' encoding='UTF-8'>
# stdlog                                  : <_io.TextIOWrapper name='BulkD83_S8.umi_dedup.genome.sorted.log' mode='a' encoding='UTF-8'>
# stdout                                  : <_io.TextIOWrapper name='BulkD83_S8.umi_dedup.genome.sorted.bam' mode='w' encoding='UTF-8'>
# subset                                  : None
# threshold                               : 1
# timeit_file                             : None
# timeit_header                           : None
# timeit_name                             : all
# tmpdir                                  : None
# umi_sep                                 : _
# umi_tag                                 : RX
# umi_tag_delim                           : None
# umi_tag_split                           : None
# umi_whitelist                           : None
# umi_whitelist_paired                    : None
# unmapped_reads                          : discard
# unpaired_reads                          : use
# whole_contig                            : False
2026-01-22 11:20:19,121 INFO command: dedup -I BulkD83_S8.genome.sorted.bam -S BulkD83_S8.umi_dedup.genome.sorted.bam -L BulkD83_S8.umi_dedup.genome.sort
ed.log --output-stats BulkD83_S8.umi_dedup.genome.sorted --method=directional --random-seed=100
2026-01-22 11:22:04,051 INFO total_umis 73103229
2026-01-22 11:22:04,051 INFO #umis 26262
2026-01-22 11:22:07,764 INFO Written out 100000 reads
2026-01-22 11:22:21,083 INFO Parsed 1000000 input reads
...
2026-01-22 12:03:29,107 INFO Reads: Input Reads: 73103229
2026-01-22 12:03:29,107 INFO Number of reads out: 12133067
2026-01-22 12:03:29,107 INFO Total number of positions deduplicated: 3940867
2026-01-22 12:03:29,107 INFO Mean number of unique UMIs per position: 4.49
2026-01-22 12:03:29,107 INFO Max. number of unique UMIs per position: 10063
# job finished in 2590 seconds at Thu Jan 22 12:03:29 2026 -- 2562.59 33.76  0.00  0.00 -- 8be2ef46-f571-4fec-a4d8-b5e157b0049c

It doesn’t look like there are other, better stats in there we could use… Here are the numbers from

samtools

, filtering out secondary alignments:

samtools view -c MAT-LIBR-19610_BulkD83_S8.genome.sorted.bam
73,103,229
samtools view -F 256 -c MAT-LIBR-19610_BulkD83_S8.genome.sorted.bam
26,062,970
samtools view -c BulkD83_S8.umi_dedup.genome.sorted.bam
12,133,067
samtools view -F 256 -c BulkD83_S8.umi_dedup.genome.sorted.bam
7,884,975

To get a feel for the UMI-aware de-duplication in the sample, and put in relation to the other numbers in the run, i.e.

raw reads

>

trimmed/length excluded reads

>

reads after rRNA removal

>

reads that can be aligned to genome

>

UMI-deduplicated alignments

it seems we’ve got to calculate 7,884,975 / 26,062,970 *100 = 30.2 % are left over sequences. Or in other words, the duplication rate of the sample was ~70%. If the values reported are 12,133,067 / 73,103,229 *100 = 16.6 % left over sequences, the implied duplication rate is not only ~84%, but also the numbers do not relate to any of the other steps mentioned above. In my mind, the numbers need to be sensible throughout the report to allow for a quality assessment of the library, so we should either: • 1) filter out secondary alignments before running UMI-tools (that would fix it outright), or • 2) run UMI-tools as is, but choose not to feed its output log to the MultiQC process, and then run

samtools view -F 256 -c

on both the original sorted and UMI-deduplicated result files, and regenerate a UMI-tools like report section (as custom content?) To me, option 1 makes a lot more sense. Since UMI-tools requires a physical BAM file (and no pipes), we would need to add a step like this:

# Filter to primary (and optionally non-supplementary) alignments
samtools view -F 256 -F 2048 -b input.bam -o primary_only.bam

and use this as the input for UMI-tools. We can then choose whether to publish or not the primary-only BAM file. If this would be the default option, but can be skipped e.g. via

--skip_primary_alignment

filtering we should have all scenarios covered, or am I mistaken?

My take on this: The diverse methods of counting are confusing. We should harmonize this to make the numbers interpretable as @Felix Krueger says. At the same time, I find that we should do this with minimal changes to the outputs, i.e., leave the alignments as they are. The multi-mapping reads may serve a purpose in some use cases. I think it would be best to manually count the number of reads (not alignments) after UMI filtering, similar to how @Thomas Danhorn suggested. If we use

samtools collate

, we can avoid the sorting for better performance. The only question to me that remains is how to represent this in the MultiQC plot? Shall we omit the UMI counts and only report the manual counts? Or shall we keep the UMI counts and report both (potentially causing even more confusion)?

Just flowing another idea: we could also be leave the BAM files in there as they are, and add another process that will filter BAM files temporarily (to primary only) but without publishing the resulting for (or only with toggle button

--umi_dedup_primary_only

), and then use that (temporary) file for the UMI-module of MultiQC. While somewhat wasteful, this would probably super easy to implement as we could still make use of the UMI module of MultiQC, avoid having to remove and replace it by a manual kludge, and make all parties happy?

Yep, sorry if didn't explain very well, but that's what I was suggesting in my option 1) 🙂

some tests are still failing…

Yeah, I just have some snapshot updates to do. But I'm mostly interested if we have consensus on the approach

Can I quickly recap: In the new behaviour, UMI-tools will now by default run with

--umitools_dedup_primary_only true

so that the deduplication stats are connected in the QC flow from:

raw reads

>

trimmed

>

non-rRNA

>

aligned

>

deduplicated

; the deduplicated BAM file contains de-duplicated primary alignments only? If someone has an interested in the looking at all multimapping de-duplicated alignments and the associated disconnected (and confusing) filtering stats, they can set a switch to

--umitools_dedup_primary_only false

, which would be the current behaviour? As far as I am concerned, I am completely happy with this

Yep, that's the idea

I don't have concerns with your solution just havent had a moment look at it

yeah I also would like to have a chance to look at it, but more towards the end of this week

All tests passing now anyway

reviewed, thanks for implementing @Jonathan Manning!

Thanks! Will respond and merge