Filter genomes with too many differences

[ArtPoon]

Filter genomes with too many differences

[ArtPoon]

require(MASS)
fit <- rlm(count ~ coldate, data=clock)
plot(clock$coldate, clock$count, cex=0.5, ylim=c(0, 50))
abline(fit, col='red')

> summary(fit)

Call: rlm(formula = count ~ coldate, data = clock)
Residuals:
      Min        1Q    Median        3Q       Max 
-12.37036  -1.54694  -0.03015   1.58198 207.25696 

Coefficients:
            Value      Std. Error t value   
(Intercept) -1046.4271     5.2548  -199.1381
coldate         0.0574     0.0003   200.6884

Residual standard error: 2.325 on 79798 degrees of freedom
  (81740 observations deleted due to missingness)
> fit$coef[2] / 30000 * 365
     coldate 
0.0006986298 

This is pretty close to the clock rate estimate of about 8e-4.

[ArtPoon]

require(MASS)
fit <- rlm(count ~ days, data=clock)
hist(residuals(fit), breaks=100, ylim=c(0, 500))
rug(residuals(fit))
abline(v=20, lty=2)

20 seems like a reasonable cutoff for declaring outliers (TODO: look into a more sophisticated method later)

plot(clock$days, clock$count, cex=0.7, pch=16,
     col=ifelse(abs(residuals(fit)) > 20, 'red', 'blue'))
abline(fit, col='red')

[ArtPoon]

This is better:

slope <- 8e-4 * 29900 / 365  # clock rate (subs/site/year -> subs/genome/day)
origin <- as.Date('2019-12-01')  # based on the literature, e.g., Lu et al. 2000, Cell 181:997
x <- seq(origin, as.Date('2020-08-01'), length.out=10)
y <- as.integer(x-origin)*slope

plot(clock$coldate, clock$count, ylim=c(0, 50), pch=16, cex=0.5, col=rgb(0,0,0,0.2))
lines(x, qpois(0.5, y), col='red')  # median of Poisson distribution with time-dependent mean
lines(x, qpois(0.005, y), col='red', lty=2)
lines(x, qpois(0.995, y), col='red', lty=2)

[ArtPoon]

Implemented by filter_outliers():

    mu = rate*29900/365.  # convert to subs/genome/day
    t0 = parse_date(origin)
    for qname, diffs, missing in iter:
        coldate = qname.split('|')[-1]
        if coldate.count('-') != 2:
            continue
        dt = (parse_date(coldate) - t0).days
        max_diffs = poisson.ppf(q=1-cutoff, mu=mu*dt)
        ndiffs = len(diffs)
        if ndiffs > max_diffs:
            # reject genome with too many differences given date
            continue
        yield qname, diffs, missing