@Sebastian Uhrig @Jack Tierney @Naiara I've just about finished with the feature @Felix Krueger asked for around truncating the paired rna-seq reads before alignment/ quantification, I wonder if you might like to review? https://github.com/nf-core/riboseq/pull/135 I must say that the idea seems odd to me, but I'll take the expertise of those more native to this field.

I absolutely agree that it is an unusual practice. I have seen this type of approach in some papers, but generally I worry that "in silico" digestion/trimming generally lowers the accuracy/quality/unique mapability of alignment for the rna-seq modality and possibly introduces new biases. So I am glad this would be optional, not the default.

Obviously I would love your review too @Ira Iosub, didn't mean to exclude you!

this is interesting to hear Ira, I heard something very different 🙂 It certainly is the case that increased read length leads to increased mappability, especially more uniquely alignable fragments, and the same is true for switching from a single-end to a paired-end library strategy. As a consequence, aligning say 2x75bp RNA-seq reads will certainly result in much more, and potentially different alignments to a corresponding Ribo-seq sample with a mean read length of e.g. 28bp. To compare like with like, this step seems to make sense to me, even though I agree that it is unfortunate that it results in a net loss in accuracy/unique mappability - but is the Ribo-seq that is the rate limiting step here. In order not to miss out on the information contained within the (potentially much longer and paired-end) RNA-seq sample, we would suggest that the RNA-seq sample is processed using the RNA-seq workflow in addition to being used for DT-analysis here. Does that make sense? @Naiara, any comments?

I remain skeptical that we gain much by making quantifications on one side of the equation less accurate in the name of consistency, but as we know, this is not my field! I'll be interested to compare the impact on results with and without this setting.

my gut feeling is that the results won’t be “better”, but probably less wrong. I agree it would be nice to see a comparison! (maybe the answer is that it doesn’t matter much which way you go after all?)

My 2 cents on this: I have little doubt that the read length differences between Ribo-seq and RNA-seq introduce systematic biases. Either the downstream tools need to model these biases (quite hard, since it's a results of mappability), or the bias needs to be mitigated where it happens - at the mapping stage. In any case, I think it would be a benefit to have the trimmed counts available, because it enables both mitigation strategies (modeling/mapping). Trimmed counts do not have to be the default. But without them, we are blind to the biases that unmatched read lengths introduce. I'm trying to think of an analysis that can elucidate whether trimmed or untrimmed RNA-seq reads give us more accurate results (without having ground truth data) ... 🤔

For me we never really compare e.g. treatment A riboseq with treatment A rnaseq (or at least that's not the main point of interest). More pertinent is the difference of responses (treatment A riboseq vs treatment B riboseq) vs (treatment A rnaseq vs treatment B rnaseq). If I were to guess, this will just make things 'messier', so any signal we see without it will be less evident. I'll do some trials tomorrow with the pipeline's full size test data at least.

I agree with Jonathan's final point that the key comparison in TE analysis is change in Ribo vs change in RNA, not absolute Ribo/RNA ratios. For differential analysis, systematic length-dependent biases should cancel out in the contrast Regardless, degrading the mappability of the RNA-Seq wouldn't be the right approach in my opinion. Most abundance metrics already account for fragment length differences in their effective length calculations. Perhaps to address any bias originating from differences is number of mappable bases, the correct approach is to calculate mappability-aware effective lengths. Essentially counting only the positions where fragments can be uniquely assigned given each modality's read length

(fyi we are already using length-adjusted abundance values in the matrices we pass to analysis)

Yes but I don't believe the account for mappability differences which is what is in question here, right?

I don't know how we would calculate anything directly for mappability tbh. But my main point is still that I don't understand why we'd need to try, since we're not comparing rnaseq quantification directly with riboseq. The reads are so different in terms of sequence content etc, the profile of potential mapping locations, I'm not sure anything would close that gap in any case.

I believe the point is that these difficult to map genes (where 30nt wont align but 150 will) are systematically missed in these analyses. They will always be given a lower TE then the actually have. I 100% agree that this doesn't mean we will be falsely calling them up or down regulated just that we are likely mischaracterising the TE of some genes

Gotcha, thanks.

I don't think my mappability proposal is important to implement in this pipeline but it could be interesting from a research perspective. I have used GenMap for something in a similar vein before.

Why don't you make an issue so we don't forget though?

Sorry I am late to this party. My 2 cents on this here as well.... I think having this as an optional feature is the right move. To add to Jonathan's point about the comparisons: When calculating Translation Efficiency (Ribo-seq abundance normalized by RNA-seq abundance), if we estimate that efficiency using different 'effective transcriptomes' (30nt mappability vs 150nt mappability), the relationship itself becomes skewed. We aren't just comparing the fold changes of two separate modalities; we are statistically interacting them. By trimming, we ensure that the RNA reads in the denominator represent the same 'mappable universe' as the Ribo reads in the numerator. I agree with Jack though that for probabilistic alignment tools , which model effective length differences mathematically, trimming might be redundant. However, for standard alignment-based counting (e.g., STAR), the software doesn't correct for this mappability difference. Thus, trimming, as far as I understand it, becomes the 'only' way to reduce this bias and prevent false negatives in regions that look 'silenced' simply because Ribo-seq reads couldn't map uniquely there. I think that, this is why workflows for specific TE tools like Xtail or Riborex often imply or recommend matching read properties to ensure valid inputs.

I wonder a little bit if we could bypass some of this by doing straight Salmon quantification (rather than going via alignments and then doing Salmon as we do now) for the purposes of the TE stuff. We'd have to set the kmer size a bit lower than normal which would introduce a measure of the same ambiguity, but it feels nicer than that hard trimming. I could enable that as another option at least. This would be in addition to the alignment-based mode, since we need the alignments elsewhere.

That at least reminds me of this approach sounds like this approach: https://www.cs.cmu.edu/~ckingsf/software/ribomap/

Sounds like they re-implemented Salmon/ Kallisto!

Its from the same last author as Salmon

I've tweaked the docs based on @Naiara’s explanation above. I'll merge today based on @Jack Tierney’s approval unless someone would like me to hold off.

since it's optional, no reason to hold off from my side