Question: should the quantifications in the *.quant and *.prob output files agree?

[DNA-Dave]

I wanted to ask if, for any given transcript that was quantified in the *.quant file, if the sum of all the read probabilities that were assigned to that transcript across all reads in the *.prob file should match the num_reads field of that transcript in the *.quant file? My understanding is that the read probabilities file contain the final EM weights for each read to its respective transcripts. However, I'm comparing the two quantities and while sometimes they match closely, other times they are significant different from each other.

For context, I have human samples with PacBio Kinnex and long-read PacBio genome sequencing data, and I've phased the Kinnex reads to each haplotype, and would like for each transcript to calculate the relative contribution of each haplotype to the total transcript expression of the sample. Therefore, I figured I could use the read-to-haplotype mappings from my Kinnex read phasing and the *.prob file to phase the readcounts generated by the *.quant file output file.

Many thanks!

[rob-p]

Hi @DNA-Dave,

Yes; they should roughly match (can differ due to some numeric/rounding issues). The first thing to do is make sure you have the latest version of oarfish as there was a bug in earlier versions where the conditional rather than the posterior probability was written.

If the discrepancy still persists, we'd be happy to take a look at some test/sample data that reproduces the issue for you to see what's going on!

Best,
Rob

[DNA-Dave]

Hi Rob,

Thanks for the quick response. I've verified that I'm using the latest oarfish version (0.9.2) and reran one of my samples and got the same issue. I've uploaded a representative *.fastq.gz file containing analysis-ready Kinnex reads of the one retested sample and the GENCODE transcripts fasta used in my oarfish run that is having the issue. I've also uploaded the oarfish log file and all outputs. One example of an isoform that has significantly different counts between the *.prob and *.quant files is ENST00000491404.2, which has a 35.73 counts in the *.quant file but only 17.74 counts in the *.prob file. I've also uploaded a Jupyter notebook with the calculations for this transcript's quantifications. The uploaded files are in a Google Drive folder and I've sent your email an access link. The fastq is quite large, and I apologize for that.

Here's the specific oarfish command I ran to run quantification:

oarfish -j 11 --verbose --reads HG00320.fastq.gz -o HG00320 --write-assignment-probs=compressed --filter-group no-filters --model-coverage --annotated gencode.v44.transcripts.fa --seq-tech pac-bio-hifi 2> oarfish.log

Many thanks for looking into this issue!

-David

[rob-p]

Thanks, David! We’ll let you know what we find and keep you in the loop on a fix!

Best,
Rob

[DNA-Dave]

Thanks Rob! Let me know if you need anything else, and looking forward to seeing what you find!

Best,
David

[rob-p]

Hi @DNA-Dave,

So, after some further investigation, I have an answer. Thank you again for the extensive testing data. That made reproducing and tracking the issue very straightforward.

The short explanation is that, to avoid writing a lot of very small probabilities, when posterior probabilities are written to the .prob file, they are first compared against a DISPLAY_THRESH. This threshold is computed based on the normalized posterior probability, so that they sum to 1 for all reported alignments of a read. The value of DISPLAY_THRESH is 0.001, and the behavior of oarfish is to discard all transcripts with a normalized posterior probability lower than this value, and to re-compute the normalized posterior probabilities among what remains.

However, in the case of this transcript (and, I presume, the others that are divergent), a substantial part of the overall read mass comes from a very large number of very low posterior probability alignments. Here is the behavior of the values in the .prob file for different values of the DISPLAY_THRESH:

DISPLAY_THRESH	sum of probs
0.001	17.7415
0.0001	35.689
0.00001	35.7252
0.0000001	35.7276

So, as you can see, while it is a somewhat rare occurrence, it does sometimes happen that a considerable amount of the total posterior read mass for a transcript comes from the sum of many small individual alignment contributions.

Regardless, this discrepancy should clearly be explained in the docs, so thank you for helping us figure this out. The remaining question is how you think may be best to proceed here.

I could lower the DISPLAY_THRESH --- as you can see, by lowering only 1 more order of magnitude, we regain most of the read mass in this case (this transcript must have many alignments close to the threshold). Of course, it is then conceptually possible to have other transcripts close to this threshold, etc. Yet, the lower we make the threshold, the larger (by an order of magnitude) the number of contributing alignments must be to make up substantial read mass.

Another option would be to remove this entirely, though then we often get many trivially low posterior probabilities (some of which even display as 0) in the .prob file, taking up space.

A third option, that trades off between the two, is to make the default threshold lower, but then allow it to be even further adjusted via a command line parameter.

In any case, we'll be sure to better document this behavior and explain the possible effect of changing this value.

What do you think is the best way forward here?

Best,
Rob

P.S. Do you have any thoughts on this @zzare-umd or @mikelove?

[DNA-Dave]

Hi Rob,

Thanks so much for the quick investigation! I appreciate it. It makes sense that the two numbers are different due to the display threshold. For me personally, I think that keeping the threshold the same or slightly lowering it but adding the option for users to set it to arbitrarily low values or even remove the threshold entirely is the best. I can't imagine most people need this granular of data, but for my purpose of calculating allele-specific counts, I do need this very granular data and so I would either set the threshold very low or remove it and would happily take the hit on storage space. I can't think of too many other situations where you would need this granular data, however, which is why I don't think its necessary to lower the threshold too much by default. I hope my perspective is helpful.

Just thinking out loud, I imagine that reads with very low posterior probabilities become a larger problem under some situations:

(1) Data is lower quality/fragmented, as happens with ONT in particular

(2) High sequencing depth with lots of reads (more chance to get reads that map to multiple isoforms)

Thanks again, and looking forward to a updated version to rerun my samples on.

Best,
David

[rob-p]

Hi David,

Thanks for the feedback! I agree. I've pushed a new release (0.9.3) that implements what is described above. The new default threshold becomes 1e-6, and it is now settable via an optional command line parameter --display-thresh. You can pass this option a floating point value, or optionally none (which will set it to machine precision for 64-bit floats). The newest version was just merged in from bioconda as well and so should be available from bioconda soon. I'll also find out where to add an appropriate section to the docs for this. Thanks again for reporting this!

--Rob