rob-p
commented
3 years ago Hi @reganhayward,
Thank you for the detailed report. It's interesting that this happens when running with STAR but not when running with selective alignment. However, salmon will attempt to solve the optimization problem with the alignments it is given, regardless of if those come from STAR or from it's built-in selective alignment. While I would generally expect these to be similar, the alignment algorithms are different; see e.g. the differences between SA/SAF & STAR here.
Nonetheless, it is possible that for a small subset of transcripts, the probabilistic allocations are so ambiguous, that you get large swings in the resulting quantification estimates based on tiny variations in where the optimization starts (which is, itself, stochastic due to the asynchronous nature of salmon's online inference phase).
One way we can test this hypothesis is as follows. You can run salmon with --numGibbsSamples 100 and -d. This will tell salmon to perform posterior Gibbs sampling (--numGibbsSamples 100) and to dump the range-factorized equivalence classes used for offline quantification (-d). The Gibbs sampling files will contain the traces for the transcripts in question over the various iterations of the sampling procedure. Transcripts where there is a tremendous amount of ambiguity will tend to have highly anti-correlated posterior samples, and similarly, if you were to consider the abundance output of these transcripts as a group, there would be a large reduction in inferential relative variance. In fact, we wrote a whole paper on this topic. Consider this example from that paper:
The x-axis is samples from the Gibbs chains, and the y-values denote the estimated number of reads assigned to both transcripts in each sample. The green line at the top is what you get if you sum the abundances of these two transcripts. The main point is that the inferential relative variance (adjusted ratio of the variance over the mean) is much smaller for the sum of these transcripts than for either individually. This is strong evidence that they are inherently uncertain given the read evidence and alignments used for quantification. The tool described in that paper, called terminus, is a tool for automatically finding such groups of transcripts.
Anyway, once you have the Gibbs samples in hand, we can walk you though how to do some assessment of these transcripts (tagging @hiraksarkar here since he's most likely to have access to scripts that will let us look at the posterior samples from individual transcripts). Similarly, if you can provide the quantification directory, we can help examine this too. If this is the case, that the posterior distributions are highly anti-correlated, it is likely that the ambiguity you are seeing is simply inherent given the alignments salmon is being provided. If you have the quantification folder resulting from the same sample using selective alignment, we could compare and contrast the two.
At that point, there are a few options depending on how deeply you want to dive. You could try to see how STAR and selective alignment are mapping differently to these transcripts. One potential difference is that STAR is a lot more happy to softclip reads, which selective alignment won't do by default (you can test the effect with the --softclipOverhangs to allow selective alignment to softclip reads that hang off the transcript end or --softclip to allow softclips anywhere). Note that selective alignment may still be a bit more conservative than STAR about softclips simply because of the nature of the scoring function it uses. This might give you a sense if one of these alignment methodologies is more consistent with your expectations in this case. Another option is to consider doing a grouping with terminus. This will reduce the set of "genes" that you can call as DE, because it will be happy to group together transcripts from different genes. However, it should help considerably in eliminating DE from highly-uncertain point estimates. Finally, you might consider performing DE with swish (cc @mikelove as he might have some input here) rather than DESeq2 (though we've typically been using swish at the transcript level rather than the gene level). Unlike DESeq2, swish will explicitly take into account the inferential uncertainty in the abundance estimates, using the Gibbs samples produced by salmon. This will allow it to avoid spurious DE calls that might otherwise occur when you have highly uncertain transcripts that, by chance, end up being assigned very different abundances in different samples / over different runs.
Sorry for the information dump, but I wanted to lay out what might be going on, how to assess it, and what some potential solutions might be. If you dive in to start investigating this, feel free to reach out in this issue along the way if you get stuck or have follow-up questions.
mikelove
reganhayward
hiraksarkar
We are using Salmon to help allocate multimapped reads in dual RNA-seq data. What this means is we combine a host and pathogen genome (bacterial genome + host transcriptome) and also use a combined annotation file as well. Basically, we can then find out which reads uniquely map to each gene/transcript from both organisms.
The issue we are having is when using Salmon in alignment-based mode (using STAR), we are observing quite substantial variations in how ambiguous (multimapping) reads are being assigned.
We use the same quant command/params and the same input files - and we run the same data twice in different output locations. The correlations between quant files are generally really good ~ 99.9x (spearman).
However, the number of allocated reads to a small proportion of transcripts varies quite a lot.
For example:
We initially noticed this when examining differentially expressed genes downstream, which were a bit unusual.
These two transcripts in particular have a high number of ambiguous counts (generally >1000) when examining the
ambig_info.tsvfile (the number of ambig reads does remain the same though)I tried additional params and runs including: Accounting for biases:
--seqBias --gcBias --posBiasIncluding bootstrap information:--numBootstraps 100 --bootstrapReprojectBoth params reduced the variation between runs slightly, but not enough to reduce the disparity.
We don't see this difference when running in selective alignment mode
We are using
v.1.3.0andsingle endreads.