RNA-seq quantification: resolving multi-mapped reads with EM

Blocks N Transcripts n Reads M Read position Within-transcript repeats Seed

Simulation: N distinct 300-bp blocks; each transcript is 2–5 blocks, so transcripts share blocks and a read from a shared block is multi-mapped. Transcript fractions p_t are uniform random; a read picks a transcript in proportion to p_t·l_t and a position along it (uniform, or exponential from one end with mean 300 bp). All reads with the same hit pattern δ (hits per transcript) collapse into one line count: δ₁,…,δ_n. EM: start by giving each transcript 1/k of every read that hits k transcripts, so q_t⁽⁰⁾ = (Σ_i δ_it/Σ_sδ_is) / R (exact already when no read is multi-mapped), then repeat q_t^(k+1) = (1/R) Σ_i δ_itq_t^(k)/l_t ⁄ Σ_sδ_isq_s^(k)/l_s (E-step: split each read among the transcripts it hits by per-base abundance; M-step: average). The transcript fraction is p_t ∝ q_t/l_t and TPM = 10⁶·p_t. Edit the pattern lines to run EM on your own counts. Keys: → step, ← back, f run all, r reset, s re-simulate.

1 Simulated transcripts

2 Reads by hit pattern (editable)

3 Estimates at the current iteration

4 E-step with the current q: how the reads of each pattern are split among the transcripts they hit

5 Convergence of the transcript fraction p_t

estimated p_t by iteration simulated truth current iteration (click the chart to jump)

RNA-seq quantification: resolving multi-mapped reads with EM

1 Simulated transcripts Show transcripts

2 Reads by hit pattern (editable)

3 Estimates at the current iteration

4 E-step with the current q: how the reads of each pattern are split among the transcripts they hit

5 Convergence of the transcript fraction pt

1 Simulated transcripts

5 Convergence of the transcript fraction p_t