We found something remarkable, and wondered if you have an explanation for it. Case 1) In short, we had some Ribo-seq data where the STAR reported that a high fraction of the reads is unmappable (~50% other), and a low number of multi-mapping reads (17%):

MULTI-MAPPING READS:
        Number of reads mapped to multiple loci |       24119780
             % of reads mapped to multiple loci |       16.96%
        Number of reads mapped to too many loci |       225856
             % of reads mapped to too many loci |       0.16%
                                  UNMAPPED READS:
  Number of reads unmapped: too many mismatches |       0
       % of reads unmapped: too many mismatches |       0.00%
            Number of reads unmapped: too short |       15408484
                 % of reads unmapped: too short |       10.83%
                Number of reads unmapped: other |       71667333
                     % of reads unmapped: other |       50.38%

Using the same version of STAR, and leaving out the parameter:

--twopassMode Basic

(which is set within Ribo-seq) gives the following numbers:

MULTI-MAPPING READS:
        Number of reads mapped to multiple loci |       99216594
             % of reads mapped to multiple loci |       69.75%
        Number of reads mapped to too many loci |       127419
             % of reads mapped to too many loci |       0.09%
                                  UNMAPPED READS:
  Number of reads unmapped: too many mismatches |       0
       % of reads unmapped: too many mismatches |       0.00%
            Number of reads unmapped: too short |       13022312
                 % of reads unmapped: too short |       9.15%
                Number of reads unmapped: other |       50600
                     % of reads unmapped: other |       0.04%

Now all of a sudden we have 70% multi-mapping reads, and no unmapped reads (other) anymore. Is this a known-phenomenon?

I'm wondering if what is happening here is that on the first pass in the two pass mode, we are getting a high number of false discovered splice junctions. In a 28nt long read at best we can get 14nt at either side of a junction to confidently call an intron. If these are allowed to span 100kb then the chances of getting false positives are very high. Mapping to this extended junction set makes it easy to get multimapped reads and hard to find a read that wont map somewhere I am always skeptical of anything that relies on Ribo-Seq data to call transcript isoforms.

you meant that practically all of 71 million reads get assigned splice junctions because of 14-14bp split. That sounds very drastic. Which reminds me, I should probably look at the reads that make it through to the BAM file to decide whether the look dramatically different between the two mapping modes

14-14 is the most confident we could be in the alignment. I'd imagine those aren't all that likely in themselves but other combinations would be more spurious and more likely to occur. Add in 3 mismatches or whatever we allow and its far more likely. But you are of course right to be skeptical because I kind of just made it up based on what I think would happen Even looking in the STAR log just above the multimapping stats should tell us more. Could you share the UNIQUE READS section for these two outputs?

Of course, here is without TwoPass:

Mapping speed, Million of reads per hour |       975.41

                          Number of input reads |       142247282
                      Average input read length |       27
                                    UNIQUE READS:
                   Uniquely mapped reads number |       29830357
                        Uniquely mapped reads % |       20.97%
                          Average mapped length |       26.87
                       Number of splices: Total |       3866566
            Number of splices: Annotated (sjdb) |       3755129
                       Number of splices: GT/AG |       3832307
                       Number of splices: GC/AG |       29697
                       Number of splices: AT/AC |       2927
               Number of splices: Non-canonical |       1635
                      Mismatch rate per base, % |       0.64%
                         Deletion rate per base |       0.01%
                        Deletion average length |       1.01
                        Insertion rate per base |       0.00%
                       Insertion average length |       1.04
                             MULTI-MAPPING READS:
        Number of reads mapped to multiple loci |       99216594
             % of reads mapped to multiple loci |       69.75%
        Number of reads mapped to too many loci |       127419
             % of reads mapped to too many loci |       0.09%
                                  UNMAPPED READS:
  Number of reads unmapped: too many mismatches |       0
       % of reads unmapped: too many mismatches |       0.00%
            Number of reads unmapped: too short |       13022312
                 % of reads unmapped: too short |       9.15%
                Number of reads unmapped: other |       50600
                     % of reads unmapped: other |       0.04%
                                  CHIMERIC READS:
                       Number of chimeric reads |       0
                            % of chimeric reads |       0.00%

So there is about 5x more novel junctions but I'm struggling to see if that adds up to the differences we are looking to explain...

I’ve got the BAM files now, will update later today (after some meetings)

I have no idea what is happening and am no expert for STAR, but here are are some puzzle pieces: Alex Dobin himself says:

unmapped-other are the reads for which the anchor seeds are not found.

The two main reasons for that are (i) poor read quality (ii) "contamination" with sequences not present in the reference genome.

So "no seed found" would mean that the seed is not in the index, which does not make sense if the same index was used for both runs, unless the length of the seed is changed (with the longer ones not found). Might be worthwhile trying different STAR versions to see if this is a regression ...

It does seem like my logic about updating the index with new splice junction was actually backwards.

I have now had a look at the alignments themselves, and the numbers are fairly dramatic:

Total reads	reads in 8 regions (5kb)
noTwoPass.   [MQ20]	   33,691,827	          4,465,411
TwoPassBasic [MQ20]	   35,862,325	          4,385,250
TwoPassBasic [MQ0]	   53,116,995	         12,672,285
noTwoPass    [MQ0]    102,632,078	         61,206,338

The same BAM files were imported only allowing uniquely mapping reads (MQ20), or allowing multimappers (MQ0) (primary alignments only). It’s clearly not that STAR cannot map to these regions generally, but it is specifically the TwoPass mode that excludes the reads from the file… I have also identified the regions where these reads get mapped to, where just 4 regions (of 5kb each) attract more than 50 Million reads. 3 out of 4 of these have a rRNA precursor or pseudogene annotated nearby.

So these are the samples that no rRNA was removed from? Those reads multimapping wouldnt be all that surprising. Combining that with Thomas' point, the two pass mode appears to have different parameters for the seed stage of the alignment and thats why they go from 'other' to multimapped

These are the samples where

sortMeRNA

has already removed 3% of reads (they have no_RNA in the name), but they clearly still contain > 50% rRNA reads. I think in this case the unique alignments are comparable in either mode (and if we want to continue only with those it will not matter much), but without two pass everything aligns (as multi-mapping), but with twoPass these reads are all of a sudden unmappable. If you ask me, me this looks like a bug and is an unwanted behaviour. It will look very different in a QC plot, potentially send you down on a contaminant hunt (which I am now on for a few days), and potentially spark a complaint with the sequencing facility/company why the samples are contaminated….

But surely you'd see that rRNA removal performed poorly and then assume those were the contaminants? I really dont think you want these reads/alignments in your downstream analysis anyway. Focus has to be on getting as many high-confidence rpfs in your final profiles as possible. Theres definitely room for some innovation on how we handle multimappers. Currently best advice if you want to ensure they are all actual RPFs is to stick with unique mappers but as Ira pointed out before, this precludes you from studying whole gene families. There as got to be a smarter way to include/exclude multimappers

I agree, we want to focus on unique mappers for the downstream things, but in my mind this behaviour of STAR is unexpected and potentially plain wrong. So in terms of the Ribo-seq (for which we considered getting a refund from the CRO) we are probably fine, but I feel like I should probably open an issue on the STAR repo….

I believe Alex Dobin has moved on to the Arc Institute so the rate of reply on the repo has dropped off.

I just saw 😞 944 open issues, hardly any replies….

So you feel this is a bug rather than just a feature of updating the index with the new splice junctions? My understanding now is that the search of the seed agains the suffix array returns too many hits because the first pass added these regions as splice junctions. Almost like a built in rRNA filter!

That’s a good point. I am currently trying ‘--outFilterMultimapNmax 1000’ (max number of multiple alignments allowed for a read: if exceeded, the read is considered unmapped) together with ‘--winAnchorMultimapNmax 2000’ to see if the few additional splice junctions exceed the multi-mapping threshold

I am still concerned about aligning to splice junctions that are identified using such short reads. I wonder if we could look into a pipeline that does a manual two pass. First pass with RNA-Seq to update the index and second pass with the Ribo-Seq to avail of the improved junction set https://groups.google.com/g/rna-star/c/rBQK-ujtSh8/m/viTPg1UFKl8J

Ha, it’s that! The injected junctions from these very short reads do indeed cause the alignments to exceed the multi-mapping filters! This is with

--twopass Basic --outFilterMultimapNmax 1000 --winAnchorMultimapNmax 2000

MULTI-MAPPING READS:
        Number of reads mapped to multiple loci |       96029346
             % of reads mapped to multiple loci |       67.51%
        Number of reads mapped to too many loci |       2553
             % of reads mapped to too many loci |       0.00%
                                  UNMAPPED READS:
  Number of reads unmapped: too many mismatches |       0
       % of reads unmapped: too many mismatches |       0.00%
            Number of reads unmapped: too short |       10176682
                 % of reads unmapped: too short |       7.15%
                Number of reads unmapped: other |       6188987
                     % of reads unmapped: other |       4.35%

OK I am glad this is not a true bug. Thanks all for your input here! (It might have helped to call the reads “Number of reads exceeding multi-map filters” rather than “unmapped”, but I guess that’s just semantics….)

Thanks for the update! I realize mapping is complicated, and there are many decisions to be made with lots of edge and corner cases, but it still is a bit unsatisfying ...

I tend to agree 🥹