四色连珠，离散矩反演与K-mer方法

休息时偶然想到这个问题：

假设有一串四色连珠，它们的总数固定且可知，但不知道各自多少个。我还需要哪些条件，才能在不知道正确答案的情况下反推其各自数量？

答案其实有很多，我可以询问谁比谁多多少，询问谁比谁少多少，只要能够联立各自独立（即“各自线性无关”）的 4 个式子组成四元一次方程组，我就可以得出各自的数量。——用矩阵的说法，这个四元一次方程组表示的矩阵是可逆的，即一定有唯一解。

说到矩阵，不得不提这个例子：

\[\begin{bmatrix} 1 & 1 & 1 & 1 \\ a & b & c & d \\ a^2 & b^2 & c^2 & d^2 \\ a^3 & b^3 & c^3 & d^3 \end{bmatrix} \begin{bmatrix} x\\y\\z\\w \end{bmatrix} = \begin{bmatrix} n\\S_1\\S_2\\S_3 \end{bmatrix}\]

它叫范德蒙德矩阵（Vandermonde Matrix），其特点是，它的四元一次方程组的形式，相当于联立它们的系数和、系数平方和、系数立方和，以及系数四次方和：

\[\begin{cases} x + y + z + w = n \\ ax + by + cz + dw = S_1 \\ a^2x + b^2y + c^2z + d^2w = S_2 \\ a^3x + b^3y + c^3z + d^3w = S_3 \end{cases}\]

其中 $x$，$y$，$z$，$w$ 为各自珠子的数量，$a$，$b$，$c$，$d$ 为各自珠子代表的系数。

你可以理解为，给不同颜色指定一个系数，然后将所有珠子的系数全取下来，对它们做上面这种运算，我就能知道它们各自出现了几次——前提是我需要知道系数有多少种。这在现实世界里不太好做，但在计算机世界，我们的CPU还是能轻松做到的。

等等？计算机？计算机世界里，我不仅能做加法减法，还可以做位平移嘛！别忘了位平移运算是最方便的，计算乘 2 或除 2 的方法。这样也能解决问题吗？

回到那个问题。我们对系数乘 1、乘 2、乘 4、乘 8，得到另一个方程组：

\[\begin{cases} ax + by + cz + dw = n \\ 2ax + 2by + 2cz + 2dw = S_1 \\ 4ax + 4by + 4cz + 4dw = S_2 \\ 8ax + 8by + 8cz + 8dw = S_3 \end{cases}\]

——不太行，它们其实表示的都是第一个方程，因此不可能有唯一解。

不过，如果它们的系数这样来设置，就能够实现我们的设想了：将系数设置为 1、16、256 和 4096（$2^0$、$2^4$、$2^8$、$2^{16}$）。虽然稍稍偏离最初的设想，但仍然保留了二进制的精髓。

……还是无法理解吗？好吧，我们将权值设得更大一些，设置为 1、100、10000 和 1000000（$10^0$、$10^2$、$10^4$、$10^6$），这其实是上面那个的十进制版本。此时你拿到结果之和，$4030201$，你立马就能知道，它们到底出现过几次！如果能拿到系数与颜色的对应关系，你甚至可以知道每种颜色的珠子出现了几次。

但……它与 K-mer 方法有哪些关系？

在生物信息学世界，DNA 里的四种碱基可以理解为一串字母序列（ATGC）、一个字符串字面量（\x41\x54\x47\x43），或者任意的 4 个数字（预先构建并对每个碱基执行 map）。而在序列分析技术中，K-mer 方法（用一个固定长度的滑动窗口读取窗口内的序列信息）经常在启发式序列搜索（BLAST）、基因组调查，以及 RNA-seq 分析中得到使用。

以基因组常用的 19-mer 为例，19 的二进制长 5 个比特，也就是说如果使用大权值的方法，它的和不会超过 20 位——一个 int32 就能搞定。但后果也是显而易见的，每一个碱基占用的比特数从 8 位上升到 16 位，大小直接翻倍。对于经过 gzip 压缩的序列，这样的序列翻倍或许能够接受。

回到那个 Vandermonde Matrix。要进一步压缩空间，或许仍然需要将它转换成 2-bit 序列——也就是分配 0、1、2、3（A、T、G、C）这四个权值。

（从某处借来一个 Python REPL）

>>> import numpy as np
>>> mer19=np.random.randint(low=0,high=4,size=19)

我们生成了一个序列，即 mer19，它代表 19-mer 序列上四种碱基（ATGC）的位置。我使用了赋值符号，因此在 REPL 里看不到生成结果。

接下来，尝试对这个数列进行操作：

>>> sum(mer19 ** 0)
np.int64(19)
>>> sum(mer19 ** 1)
np.int64(28)
>>> sum(mer19 ** 2)
np.int64(66)
>>> sum(mer19 ** 3)
np.int64(172)

此时的方程组为

\[\begin{cases} x + y + z + w = 19 \\ y + 2z + 3w = 28 \\ y + 4z + 9w = 66 \\ y + 8z + 27w = 172 \end{cases}\]

解得

\[\boxed{ \begin{aligned}\begin{cases} x = 5 \\ y = 5 \\ z = 4 \\ w = 5 \end{cases} \end{aligned}}\]

我还真尝试过解这个方程组，但已经相当费力，而且还解错了。感慨接触编程近 5 载，投身代码海洋一年有余，数学能力反而倒退好多……

不想吃这等苦头的话，我们还有第二种方法。

>>> V = np.vander(np.array([0,1,2,3]), N=4, increasing=True).transpose()
>>> V
array([[ 1,  1,  1,  1],
       [ 0,  1,  2,  3],
       [ 0,  1,  4,  9],
       [ 0,  1,  8, 27]])
>>> s = np.array([19, 28, 66, 172])
>>> # 解 Vx=s
>>> np.linalg.solve(V,s)
array([5., 5., 4., 5.])

可以看到，结果是浮点数。你当然可以手动改变数据类型，但 SymPy 其实更加适合这种计算。

>>> import sympy as sp
>>> V = sp.Matrix([
...     [1, 1, 1, 1],
...     [0, 1, 2, 3],
...     [0, 1, 4, 9],
...     [0, 1, 8, 27]
... ])
>>> s = sp.Matrix([19,28,66,172])
>>> V.LUsolve(s)
Matrix([
[5],
[5],
[4],
[5]])

如果只想做 GC 比，其实只需要套用一下计算公式：

\[GC\%=\frac{-7S_1+9S_2-2S_3}{6k}\]

即得

\[GC\%=47.37\%\]

……终于有点尊严感了。

补充1

在 \(\{0, 1, 2, 3\}\) 映射下，各项计数分别为

\[\boxed{ \begin{aligned} x_0 &= k - \frac{11}{6}S_1 + S_2 - \frac{1}{6}S_3 \\ x_1 &= 3S_1 - \frac{5}{2}S_2 + \frac{1}{2}S_3 \\ x_2 &= -\frac{3}{2}S_1 + 2S_2 - \frac{1}{2}S_3 \\ x_3 &= \frac{1}{3}S_1 - \frac{1}{2}S_2 + \frac{1}{6}S_3 \end{aligned}}\]

如果只想算 GC 比，那么将 G 和 C 映射到 1 和 2 其实是最爽的：

\[\boxed{GC\% = \frac{GC}{k} = \frac{3S_1 - S_2}{2k}}\]

因此，如果你预先将序列按 {'A':0, 'G':1, 'C':2, 'T':3} 映射，你只需要维护一个 K-mer 滑窗与计数器，计算整个数字序列的和与平方和，我就可以演算出所有 K-mer 的 GC 比；更进一步地，我可以通过平方和计算所有 K-mer 中各碱基的出现频率，通过立方和对结果进行核实验证。

补充2

手动解算过程：

对于方程组

\[(1) &\quad x + y + z + w = 19 \\ (2) &\quad y + 2z + 3w = 28 \\ (3) &\quad y + 4z + 9w = 66 \\ (4) &\quad y + 8z + 27w = 172\]

\((3) - (2)\)：

\[\begin{aligned} (5)\quad2z + 6w = 38 \end{aligned}\]

\((4) - (3)\)：

\[\begin{aligned} (6)\quad4z + 18w = 106 \end{aligned}\]

\((6)·2 - (5)\)：

\[6w = 30\\ w = 5\]

将 \(w\) 回代至 \((5)\)，得

\[z = 4\]

将 \(z\)、\(w\) 回代至 \((2)\)，得

\[y = 5\]

将 \(y\)、\(z\)、\(w\) 回代至 \((1)\)，得

\[x = 5\]

即得

\[\boxed{ \begin{aligned}\begin{cases} x = 5 \\ y = 5 \\ z = 4 \\ w = 5 \end{cases} \end{aligned}}\]