PacBio consensus quality

[phleroy]

Hello

I have sequenced a BAC clone with PacBio RSII
To make the assembly I used Facon through pbbioconda and for polishing I used quiver

To have an estimation of the consensus quality I re map the original bam reads file against the consensus

How to estimate a mean quality value, in other world a consensus Phred score for the base calls of the consensus ... :-)

Thank you in advance
Philippe

[SNPsaurus]

We polish with arrow and just list one of the outputs as fastq "-o sample_consensus.fastq" and it generates a fastq file with a consensus for each contig and the quality score. You might check if quiver has the same option, or switch to arrow (here's a blog about doing so https://dazzlerblog.wordpress.com/tag/arrow/ ).

[phleroy]

Thank you very much for this suggestion
We do have the possibility to obtain a fastq file with quiver with the option -o out.fastq as you have mentionned for arrow

The question, is then, how you recover the mean QV for the consensus?

into the fastq file we can see :
@000000F|quiver
ATCATTGTTACTACTAGAGGAAGAATCTTTCTTG ...
+
"RQQPQQQRQQQQQSRRQSTSSQRQRSSSRRRQQRSRRRQRSRQ ...

I guess the quality value for each consensus nucleotide is the second line ? but how to calculate it ?

Thank you again for any help
Philippe

[Magdoll]

You can convert the Phred QV scores to probabilities then sum over the probabilities over the entire sequence to get the expected number of errors.

You can use this Python script to calculate expected acc from a FASTQ files (though this is in a repo meant for PacBio transcriptome data, this script is generic):

https://github.com/Magdoll/cDNA_Cupcake/wiki/Sequence-Manipulation-Wiki#expacc

[phleroy]

Thank you so much, I will try this option as soon as possible and tell you :-)

[phleroy]

I tried the python script (calc_expected_accuracy_from_fastq.py) on our fastq consensus sequence which was obtained with quiver and obtained as expected the "expected_accurancy" which was : expected_accuracy=0.997

In a previous analysis I used two smrtlink python scripts to estimate the mean_QV
- summarize_coverage.py to obtain a alignment summary gff file
- polished_assembly.py to obtain the csv file which gives the a mean_qv of 48.65

I have the feeling that the two values estimate different metrics ? I am not a specialist of this area and I am curious to have any remarks or suggestion

Nevertheless, these two values : mean_qv and expected_accuracy should give an estimation of the quality of the consensus assembly. I just need to understand precisely what interpretation to have for each value

Thank you in advance
Philippe

[rhall]

If you assemble a set of reads, then use them to polish the assembly, there is no way to measure any truly meaningful consensus quality without an orthogonal datatype, or knowledge of ground truth. The expected accuracy from the fastq that results from polishing is highly dependent on the consensus algorithm and may not be a true indication of the quality of the consensus.