Starting with case 2. I assume you have experience...
# riboseq
j
Starting with case 2. I assume you have experienced better success with other runs using SortMeRNA? Is the read length distribution of the samples that just high 3% rRNA any different to a typical sample? My guess is that the shorter read lengths suffer with
window_size
set to the default of 18. Perhaps we drop the window size to 14 or 15 to increase the number of windows per read. Perhaps we should test RiboDetector also? https://doi.org/10.1093/nar/gkac112
j
That's a really good point
f
Setting the seed length to
-L 14
makes exactly no difference
j
Interesting! Not even to the E-values?
f
default:
Copy code
Parameters summary: 
    Reference file: rDNA.fa
        Seed length = 18
        Pass 1 = 18, Pass 2 = 9, Pass 3 = 3
        Gumbel lambda = 0.575284
        Gumbel K = 0.307525
        Minimal SW score based on E-value = 53
    Number of seeds = 2
    Edges = 4
    SW match = 2
    SW mismatch = -3
    SW gap open penalty = 5
    SW gap extend penalty = 2
    SW ambiguous nucleotide = -3
    SQ tags are not output
    Number of alignment processing threads = 32
    Reads file: file.non_rRNA.fastq.gz
    Total reads = 142247282

 Results:
    Total reads passing E-value threshold = 21172086 (14.88)
    Total reads failing E-value threshold = 121075196 (85.12)
    Minimum read length = 20
    Maximum read length = 44
    Mean read length    = 27

 Coverage by database:
    rDNA.fa             14.85

 Wed Jun 11 11:02:42 2025
`-L 14`:
Copy code
Parameters summary: 
    Reference file: rDNA.fa
        Seed length = 14
        Pass 1 = 18, Pass 2 = 9, Pass 3 = 3
        Gumbel lambda = 0.575284
        Gumbel K = 0.307525
        Minimal SW score based on E-value = 53
    Number of seeds = 2
    Edges = 4
    SW match = 2
    SW mismatch = -3
    SW gap open penalty = 5
    SW gap extend penalty = 2
    SW ambiguous nucleotide = -3
    SQ tags are not output
    Number of alignment processing threads = 32
    Reads file: file.non_rRNA.fastq.gz
    Total reads = 142247282

 Results:
    Total reads passing E-value threshold = 21172086 (14.88)
    Total reads failing E-value threshold = 121075196 (85.12)
    Minimum read length = 20
    Maximum read length = 44
    Mean read length    = 27

 Coverage by database:
    rDNA.fa             14.86
(this is the same rDNA.fa file that produces 51% alignments using Bowtie2)
j
So its that the minimum score is too high. Basically we would need a 27nt long read to align perfectly to reach the minimum score (27 x 2 = 54). Looking at the code it doesnt seem like seed length is used in the E-value score at all (which i thought it would based on how BLAST works). I think I can see why the model breaks down if the read length distribution has high variance or is skewed (but this doesnt appear to be the case here). Regarding why it could have issues with shorter reads is that the algorithm expects the length of the high scoring pairs to be shorter than the query sequence length. It should be possible to get multiple high scoring pairs per sequence for the Karlin-Altschul algorithm to work. This is pure ChatGPT but it does kind of track with what I remember from studying these algorithms. I got Claude to make an Artefact. https://claude.ai/public/artifacts/e3e9eb24-d0d7-494c-96b9-f791f128cef5