看到一篇关于kallisto软件原理的文章，内容通过图文的形式详细地介绍了软件的工作原理，让读者能够直观地理解以便更好地使用软件，内容通俗易懂值得学习与分享。内容如下：
kallisto是一款基于 k-mer 方法从 RNA-Seq 数据中估算转录本表达丰度的软件。与之前使用的方法相比，在速度和内存使用方面有了显著改进，同时保持了类似的准确性。kallisto非常轻捷且方便， 直接以 TPM 的形式输出定量结果。但它是如何工作的呢？
传统的定量方法依赖于比对，即对 RNA-Seq 测序的reads与参考基因组先进行比对。reads被分配到基因组中的一个位置，基因或转录本的表达值通过计算重叠reads的数量而得到。
kallisto 的核心思想是依赖于伪比对，它不识别reads在转录本中的位置，只确定它们可能的转录本来源。因此，它避免了对每个read与参考基因组进行比对。实际上，它只使用转录组序列而不是整个基因组。
在分析之前，构建 kallisto 索引是必不可少的。kallisto 首先从转录组中找到的所有 k-mer 来构建一个 colored de Bruijn graph (T-DBG)。图中的每个节点对应一个 k-mer，并以不同的颜色区分它们属于哪些转录本。图中具有相同颜色的线性部分对应于重叠群或转录本。一旦构建了 T-DBG，kallisto 就会存储每个 k-mer 与其转录本来源以及在转录本中的位置之间的映射关系。这一步相同的基因组只需做一次，并且依赖于提供的注释文件。
如上示意图所示，对于给定的测序样本，kallisto 将每个read分解为k-mer，并用这些 k-mer 在 T-DBG 中找到一条覆盖路径。这条转录组图的覆盖路径，其中一条路径对应一个转录本，为每个 k-mer 生成 k-兼容类，即节点上潜在转录本来源的集合。通过取其 k-mer 的 k-兼容类的交集，可以得到一个read的潜在转录本来源。为了加快伪比对的速度，kallisto 会跳过冗余的 k-mer，因为相邻的 k-mer 往往属于同一个转录本。然后，kallisto 使用 expectation-maximization (EM) 算法优化以下 RNA-Seq 似然函数。
在上面函数中，F 是片段(fragment或read)的集合，T 是转录本的集合， 是转录本 t 的有效长度，
是一个兼容性矩阵，定义为如果片段 f 与 t 兼容则为 1，否则为 0。参数
是从转录本 t 中选择读段的概率。这些
是我们感兴趣的参数，因为它们代表转录本的丰度或相对表达量。
为了加快速度，兼容性矩阵被分解成等价类。等价类由所有与相同转录本子集兼容的读段组成。EM 算法应用于等价类，每个 将被优化以得到最大化的似然。
下面用一个小例子来说明 kallisto 涉及的核心步骤。例如，一个只有 3 个转录本的小基因组，假设 RNA-Seq 实验产生了 4 个reads，如下图所示。
首先，是用kallisto构建 T-DBG 图索引，如下图。所有转录本序列被分解成 k-mer (这里k=5) 以T-DBG构建图。这里仅为了说明，所以没有在图中展示所有的节点/k-mer。重点是每个不同的转录本会在图中形成一条不同的路径。需要注意的是，kallisto的定量方法属于非链特异性。
构建好索引，就可以分析测序样本的 4 个reads了。它们被分解成 k-mer (这里 k=5)，然后使用预先构建的索引确定每个 k-mer 的 k-兼容类。然后，计算每个read的 k-兼容类。例如，对于read1，其所有 k-mer 的 k-兼容类的交集表明它可能来自转录本1或转录本2。
当4个reads都完成后，就可以构建兼容性矩阵，这是需要最大化的 RNA-Seq 似然函数的一部分。在这个例子中，将 EM 算法应用于4个reads没有问题，但在现实中，将其应用于数百万个reads就会很慢。
因此，兼容性矩阵被折叠成等价类，并为每个类计算一个计数 (有多少reads由这个等价类代表)。EM 算法使用这个折叠的信息来最大化等价的 RNA-Seq 似然函数并优化
。
EM 算法会改变 以最大化 L(α)，即它会尝试每个转录本的不同概率值，以获得尽可能大的似然。这里用一个直观的方式来说明 EM 算法的作用，来直观的认识一下。回到在这个例子，将3个概率
(即从每个转录本中选择reads的概率) 初始化为相等的值0.33，得到的 L(α) 为 8.4e-05。
[0.33, 0.33, 0.33] ==> 8.4e-05
[0.45, 0.15, 0.4] ==> 0.00011
[0.98, 0.01, 0.01] ==> 0.00042
[0.9, 0.05, 0.05] ==> 0.00037
[0.05, 0.9, 0.05] ==> 1.5e-05
[0.05, 0.05, 0.9] ==> 2.2e-06
[0.15, 0.45, 0.4] == >3.3e-05
如果将丰度值从 [0.33, 0.33, 0.33] 改为 [0.98, 0.01, 0.01] (更有可能来自转录本1)，得到的 L(α) 为 4.2e-4，这比 8.4e-05 大。
在这个例子中，考虑到这4个reads，它们很可能来自转录本1，该转录本是3个转录本中丰度最高的，所以认为 [0.98, 0.01, 0.01] 是这些转录本的相对丰度。在 kallisto 中，对于每个转录本t， > 0.01（其中 N 是读段总数）的变化小于 1% 时，EM 算法停止。
参考资料：https://bioinfo.iric.ca/understanding-how-kallisto-works
