Bracken Bayesian abundance re-estimation — morie_bracken

A k-mer classifier assigns each read to the lowest node whose k-mers identify it uniquely, so reads from species sharing genome stop at a common ancestor. Reading species counts straight off the classifier therefore understates every species by a DIFFERENT amount, which is a distorted ranking rather than a noisy one.

Usage

morie_bracken_abundance(
  kraken_output,
  kmer_distribution,
  max_iter = 1000L,
  tol = 1e-12
)

Arguments

kraken_output: Reads assigned to each classification node.
kmer_distribution: P(node | species); columns sum to 1.
max_iter, tol: Expectation-maximisation controls.

Value

A list with estimate, fractions, reads_reassigned, iterations, converged, identifiable, log_likelihood.

Details

Bracken redistributes them by expectation-maximisation over \(P(j \mid i) \propto P(i \mid j)\theta_j\). What it cannot do is invent information: species with identical k-mer columns leave the likelihood flat along the direction trading one for the other, and identifiable reports that rather than returning an arbitrary point on a ridge.

References

Lu J, Breitwieser FP, Thielen P, Salzberg SL (2017) PeerJ Comput Sci 3:e104, doi:10.7717/peerj-cs.104 .

Examples

P <- rbind(c(0.6, 0, 0), c(0, 0.5, 0), c(0, 0, 0.4), c(0.4, 0.5,
    0), c(0, 0, 0.6))
reads <- c(1200, 900, 700, 1500, 1100)
morie_bracken_abundance(reads, P)
#> $estimate
#> [1] 1913.85 1686.15 1800.00
#> 
#> $fractions
#> [1] 0.3544167 0.3122500 0.3333333
#> 
#> $naive_species_reads
#> [1] 1200  900 1800
#> 
#> $reads_reassigned
#> [1] 1500
#> 
#> $reads_unassignable
#> [1] 0
#> 
#> $total_reads
#> [1] 5400
#> 
#> $iterations
#> [1] 28
#> 
#> $converged
#> [1] TRUE
#> 
#> $log_likelihood
#> [1] -8526.47
#> 
#> $identifiable
#> [1] TRUE
#> 
#> $max_column_cosine
#> [1] 0.3922323
#> 
#> $n_nodes
#> [1] 5
#> 
#> $n
#> [1] 3
#> 
#> $warnings
#> character(0)
#> 
#> $method
#> [1] "Bracken Bayesian abundance re-estimation"
#>