Fixes to quality plots with fwd & rev reads by jgrg · Pull Request #2116 · samtools/samtools

jgrg · 2024-09-18T16:08:29Z

Fixes warning:

"38911_1#35_F0xB00-quals2.gp" line 932: warning: y2 axis range undefined or overflow, resetting to [0:0]

which resulted in missing ticks on quals2 plot.

Fixes quals-hm broken plot layout:

to:

and changes heatmap to a (hopefully) more interpretable one which avoids dark zones between blue and green and between green and red.

Fixes broken table HTML in gallery file and broken links when file names contain characters which need to be URI escaped.

…nd heatmap

…lipped

jgrg · 2024-09-19T09:12:51Z

Should add this commit too:
Reduce width of plot so that "Number of bases" heatmap label is not clipped

jkbonfield · 2024-09-23T10:08:59Z

Big improvement on layout obviously.

Overall I like the colours, but there's one slight issue that going through very light colours in two places (the very top and around aqua) can lead to confusion when you have a very dotty image instead of the more continuous gradations we see on NovaSeqs. See these examples.

The second picture, to my red/green colourblind eyes, I don't easily distinguish the two light shades.

The blue colours are Q10 are easier to see though. I wonder if we just need a slightly truncated colour-wheel?

jgrg · 2024-09-25T10:21:27Z

Yes, the rainbow heatmap is not ideal for this. This one was based on a MatLab heatmap in the gnuplot documentation. It's OK for smooth contours, like plotting mathematical functions, but where there are adjacent areas from different parts of the heatmap we need one where it's obvious which is larger. I'll have to a go at choosing one which is consistently dark > light alongside any colour change.

jkbonfield · 2024-09-25T10:31:55Z

Thanks James

jgrg · 2024-09-27T11:42:21Z

I think we want one of the perceptually uniform matplotlib colourmaps.
Choosing "plasma", but squashing the dark end a touch to make low values more visible, I get:

jkbonfield · 2024-09-30T08:48:57Z

It's a better gradation, but personally I think it's a bit too washed out at low frequencies.
Eg compare:

With

I agree with what you were attempting to do by moving away from a black background, but this feels too light to me. Thoughts on this change? (Just a quick hacked up tweak of the lowest two thresholds)

diff --git a/misc/plot-bamstats b/misc/plot-bamstats
index 3343a7fa..f4e0acf8 100755
--- a/misc/plot-bamstats
+++ b/misc/plot-bamstats
@@ -825,8 +825,8 @@ sub plot_qualities
             # Dark end squished to make low values more distinguishable.
             set palette model RGB
             set palette defined ( \\
-                 0 '#0c0887' \\
-              ,  5 '#4b03a1' \\
+                 0 '#0c0847' \\
+              ,  5 '#4b0381' \\
               , 10 '#7d03a8' \\
               , 20 '#a82296' \\
               , 30 '#cb4679' \\

jkbonfield · 2024-09-30T09:01:10Z

Or alternatively moving the 0 entry to 1 and making a new 0 (even black) really highlights the difference going from zero or negligble amounts to significant.

Have a play with a few data sets I think

jkbonfield · 2024-09-30T09:05:36Z

Ultima Genomics demonstrates other inadequacies in our plotting and data gathering too. (Don't even try it on ONT!)

With a new black at level 0 and dark blue at level 1 shows the low-level background.

jgrg · 2024-09-30T09:17:03Z

Or alternatively moving the 0 entry to 1 and making a new 0 (even black) really highlights the difference going from zero or negligble amounts to significant.

Yes, that sounds like a good idea.

jkbonfield · 2024-09-30T09:22:17Z

I suggest you play with some single ended runs too. I see this change also broke things there, removing the xtics and colourwheel. A stab at fixing it is below (with other colour hacks as I was just experimenting):

diff --git a/misc/plot-bamstats b/misc/plot-bamstats
index 3343a7fa..d6bd12a6 100755
--- a/misc/plot-bamstats
+++ b/misc/plot-bamstats
@@ -715,7 +715,8 @@ sub plot_qualities
     # Average, mean and quality percentiles per cycle, forward and reverse reads in separate plots
     $args = get_defaults($opts,"$$opts{prefix}quals2.png",wh=>$is_paired ? '700,500' : '600,400');
     $fh = $$args{fh};
-    my $pos_size = $is_paired ? " set rmargin 0; set lmargin 0; set tmargin 0; set bmargin 0; set origin 0.1,0.1; set size 0.4,0.8" : '';
+#    my $pos_size = $is_paired ? " set rmargin 0; set lmargin 0; set tmargin 0; set bmargin 0; set origin 0.1,0.1; set size 0.4,0.8" : '';
+    my $pos_size = " set rmargin 0; set lmargin 0; set tmargin 0; set bmargin 0; set origin 0.1,0.1; set size 0.4,0.8";
     print $fh qq[
             $$args{terminal}
             set output "$$args{img}"
@@ -812,7 +813,12 @@ sub plot_qualities
     $fh = $$args{fh};
     my $max = defined $lmax && $lmax > $fmax ? $lmax : $fmax;
     my $ytics = generate_ticks([ map { $_->[0] } @ffq ], $is_paired);
-    $pos_size = $is_paired ? " set origin 0.05,0.5\n set size 0.85,0.5" : '';
+    my $colorbox = "set colorbox vertical user origin first ($nquals+1),0 size screen 0.025,0.812";
+    $pos_size = $is_paired
+       ? " set origin 0.05,0.5\n set size 0.85,0.5" :
+         " set origin 0.05,0.05\n set size 0.85,1.0\n$colorbox" ;
+    my $set_xtics = $is_paired ? "unset" : "set";
+
     print $fh qq[
             $$args{terminal}
             set output "$$args{img}"
@@ -825,8 +831,9 @@ sub plot_qualities
             # Dark end squished to make low values more distinguishable.
             set palette model RGB
             set palette defined ( \\
-                 0 '#0c0887' \\
-              ,  5 '#4b03a1' \\
+                 0 '#000000' \\
+              ,  1 '#0c0847' \\
+              ,  5 '#4b0381' \\
               , 10 '#7d03a8' \\
               , 20 '#a82296' \\
               , 30 '#cb4679' \\
@@ -850,7 +857,7 @@ sub plot_qualities
             set ylabel "Cycle (fwd reads)" offset character -1,0
             unset ytics
             set ytics ($ytics)
-            unset xtics
+            $set_xtics xtics
             set title "$$args{title}" noenhanced
             splot '-' matrix with image
         ];
@@ -873,7 +880,7 @@ sub plot_qualities
                 set ytics ($ytics)
                 set xrange  [0:$nquals]
                 set xtics
-                set colorbox vertical user origin first ($nquals+1),0 size screen 0.025,0.812
+                $colorbox
                 set cblabel "Number of bases"
                 splot '-' matrix with image
             ];

jkbonfield · 2024-09-30T10:26:08Z

The yticks could be improved when we start looking at longer reads too. Eg:

     my ($cycles, $skip_even) = @_;
 
     my @marks;
+    my $div = 10;
+    while (scalar(@$cycles)/$div > 20) {
+       $div *= 2;
+    }
+    $div = nice_num($div);
     foreach my $v (@$cycles) {
-        if ($v % 10 == 0) {
+        if ($v % $div == 0) {
             if ($skip_even) {
                 my $label = $v % 20 ? $v : '';
                 push @marks, qq{"$label" $v};

Although frankly samtools stats itself is pretty broken once you get to pacbio lengths and longer as it's inappropriate for it to record data for every cycle and it just explodes in memory and file size.

A Revio example though:

jkbonfield · 2024-10-24T10:46:56Z

Any update on this James?
Thanks.

jgrg · 2024-10-29T10:22:47Z

Sorry, I haven't had time to look at this. I was hoping to submit a quick fix and move on, and not end up rewriting it, in which case I'd probably choose a plotting library other than gnuplot! I'd also want to add a diverse set of example input, but I only have examples of NovaSeqX stats. We don't use bamstats for PacBio data.

jkbonfield · 2024-10-29T11:35:53Z

Your improvements are clearly significantly improving things, so thanks.

If it helps I'm happy to pick up from where you got to and add my changes above too in a second commit. It's obviously better than it languishing due to insufficient data or time for you to progress on it.

Fix single-ended run plots. Reduce rate of tick labels for long-reads. Adjust colour wheel to bring out more fidelity with very low qualities so we have better distinction between all zero and some low qual. Replaces / fixes samtools#2116

Fix single-ended run plots. Reduce rate of tick labels for long-reads. Adjust colour wheel to bring out more fidelity with very low qualities so we have better distinction between all zero and some low qual. Replaces / fixes samtools#2116 (cherry picked from commit aa8cb59)

jgrg added 3 commits September 18, 2024 16:41

Fix generation of HTML gallery file

ba58089

Fix y-axis ticks on quals2 plot and improve fwd/rev quals-hm layout a…

31c5a6f

…nd heatmap

Reduce width of plot so that "Number of bases" heatmap label is not c…

05d8072

…lipped

daviesrob assigned jkbonfield Sep 19, 2024

Merge branch 'samtools:develop' into develop

197d8af

Use "plasma" heatmap

b54b88a

jkbonfield mentioned this pull request Nov 18, 2024

Improve plot-bamstats qual-hm #2143

Merged

whitwham closed this in aa8cb59 Dec 11, 2024

Uh oh!

Conversation

jgrg commented Sep 18, 2024

Uh oh!

jgrg commented Sep 19, 2024

Uh oh!

jkbonfield commented Sep 23, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jgrg commented Sep 25, 2024

Uh oh!

jkbonfield commented Sep 25, 2024

Uh oh!

jgrg commented Sep 27, 2024

Uh oh!

jkbonfield commented Sep 30, 2024

Uh oh!

jkbonfield commented Sep 30, 2024

Uh oh!

jkbonfield commented Sep 30, 2024

Uh oh!

jgrg commented Sep 30, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jkbonfield commented Sep 30, 2024

Uh oh!

jkbonfield commented Sep 30, 2024

Uh oh!

jkbonfield commented Oct 24, 2024

Uh oh!

jgrg commented Oct 29, 2024

Uh oh!

jkbonfield commented Oct 29, 2024

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

jkbonfield commented Sep 23, 2024 •

edited

Loading

jgrg commented Sep 30, 2024 •

edited

Loading