Is there a minimum read count or number of genes for Salmon bias correction?

[rob-p]

Question from Biostars:

I am trying to apply Salmon to a very small (artificial in-silico) genome for testing of an in-development pipeline. I am wondering if there is a limit to bias correction options (particularly --gcbias, --posbias but also --seqbias) in how few reads they need in order to be expected to operate decently? For example, is 100,000 reads too few for them to work? One million reads? I notice that seqbias is documented to use the first million reads: is that a minimum for proper functioning?

Similarly, are there restrictions on the number of distinct genes and/or transcripts needed for these corrections to be meaningful? Would you expect them to operate adequately if only a few dozen transcripts were expressed?

On a related note, is there any output from Salmon about the size of the observed biases and amount of 'correction' applied? I'm interested in values that could be compared between different samples that would indicate how much of a bias is present or how much Salmon was able to do to compensate.

These are stranded, paired-end bulk RNA-seq data if that is relevant.

[rob-p]

Original answer from Biostars by @mikelove:

I think 1 million is a good idea, because you want to be able to estimate parameters based on rare observations as well, e.g. low and high GC content fragments. I don't know if we did sensitivity analysis on how few fragments are needed.

The fragments should be from a random collection of genes/transcripts, this is pretty important or else the bias will be very skewed/wrong.

You can compare the effective transcript length (which has the bias baked in) and the nominal transcript lengths - mean fragment length (L_txp - mu_FLD). The more deviation, the more bias.

[rob-p]

Some additional thoughts on Mike's answer (which I agree with).

• On the use of 1 million

The use of 1 million is certainly a "rule of thumb". We have not done systematic testing across a wide range of values on how many fragments you need to observe to have a reasonable estimation of these biases. This is because most bulk samples are on the order of 10s of millions of reads (or more), and these are more than sufficient to learn the model. Nonetheless, I think that 1 million is a reasonable rule of thumb, and if the samples are randomly drawn, we can usually begin to see some stabilization form in the model at this point.

• On the randomness of the genes / transcripts being sampled

As Mike points out, it is critical that the sample of fragments be unbiased with respect to what is actually being expressed in the sample. If you draw a million fragments, but they are all from some particularly small subset of of expressed transcripts, then this can bias the model such that it will not properly represent GC or sequence specific biases (expressed) transcriptome wide.

• Measuring the bias

Mike's suggestion is absolutely the simplest and most straightforward here. The bias terms are folded into the effective length of the transcript, and you can measure the overall bias by looking at how the reported effective length differs from the "standard" effective length. In particular, you might look at the ratio of these values effective-length_reported / effective-length_standard. If this value is large (much > 1) then we are biased toward oversampling from this transcript and if this value is small (much < 1) then we are biased against sampling from this transcript. The individual parameters of the GC bias, positional bias, and sequence specific bias models are written out by salmon as well, though they are serialized in a binary format briefly described in the documentation, and we don't have an off-the-shelf tool to parse these models and provide a simple score from them.