ID/DIV plots don't shown sequences for certain edge cases

[ressy]

Hi Chaim,

I've found some edge cases where no tile is drawn in the ID/DIV plots even though there's a sequence there.

If I use this example id-div.tab file:

sequence_id v_gene    germ_div mab 
1           IGHV1-AFS 0        80  
2           IGHV1-AFS 5        85  
3           IGHV1-AFS 10       90  
4           IGHV1-AFS 15       95  
5           IGHV1-AFS 20       100

I then get tiles drawn for sequences 2 and 4, but not 1, 3, or 5:

1 and 5 are apparently just getting removed by ggplot's X/Y scaling. The docs say:

Note that setting limits on positional scales will remove data outside of the limits. If the purpose is to zoom, use the limit argument in the coordinate system (see coord_cartesian()).

(I guess they define the limits as an open interval but just don't say so?) If I nudge the limits a little or switch to an equivalent coord_cartesian in plot_all, it then shows those two points:

Sequence 3 is due to a weirder issue. I think it's because the MASS::kde2d call the plot function uses splits the contribution of that one point perfectly between two bins in X and two bins in Y, so it smears it across four bins. Then the color scaling in scale_fill_gradientn misses it, since the maximum value around that region is just a tad too low to apparently count as one sequence's worth of stuff. If I tweak the .../2 on this line I can make the tile get drawn, but I'm not sure about the details of these variables:

b     <- (sum(g$z) / length(data$germ_div))/2

If I instead do .../5 it shows up:

..but that throws off the definition of what counts as one sequence, so I'm not sure what the right approach would be.

Since I'm not familiar with these kernel density estimations here's a quick demo of what I'm trying to describe, with one observation on a 10x10 grid positioned either "just wrong" (split evenly across four bins) or "just right" (centered on one bin):

g <- MASS::kde2d(0, 0, h = 1, n = 10, lims = c(-1, 1, -1, 1))
pheatmap::pheatmap(g$z, cluster_cols = F, cluster_rows = F)
g <- MASS::kde2d(1/9, 1/9, h = 1, n = 10, lims = c(-1, 1, -1, 1))
pheatmap::pheatmap(g$z, cluster_cols = F, cluster_rows = F)

In practical terms this is less of a problem than that clipping issue above, though, since the individual point is still there and can be selected in the interactive plot (it gets drawn when the plot adds the individual sequences for the selected island). In contrast I don't think I can access sequences 1 and 5 during island selection unless the limits of the plot area are expanded a little.

[scharch]

Good to know about X/Y scaling, thanks - I will make that change.
I am not surprised about the problems with the KDE, I've never really been happy with it. The scale for sequence count is something I hacked together manually, and I can't swear to how accurate it is (though it should at least give a decent feel for relative counts). Probably the solution is to add a layer with the sparsest points...

note to future self: https://stackoverflow.com/questions/13094827/how-to-reproduce-smoothscatters-outlier-plotting-in-ggplot

[ressy]

Thanks!

I had a thought about the KDE thing-- since you start with counts and then loop back around to an estimate of counts, why not bin them directly and plot the binned counts? I played around with using cut on both axes and it made some promising-looking results:

  xbreaks <- seq(0, 60, length.out = 100)
  ybreaks <- seq(40, 100, length.out = 100)
  densf <- as.data.frame(table(
    cut(data$germ_div, breaks = xbreaks, include.lowest = TRUE),
    cut(data[[native]], breaks = ybreaks, include.lowest = TRUE)))
  densf$x <- xbreaks[densf$Var1]
  densf$y <- ybreaks[densf$Var2]
  densf$z <- densf$Freq

  b <- 1
  t <- b * 10^4
  if ( max(densf$z) > t ) {
	       t <- b * 10^ceiling( log10( max(densf$z) ) )
	}
	r     <- 10^seq(log10(b), log10(t), 1)

Before:

After:

It makes the contours more jagged so maybe you do still want some of that smoothing effect though.