基因组序列片段压缩与无性二倍体谱系研究

### 基因组序列片段压缩与无性二倍体谱系研究 在基因组研究中，数据的高效存储和分析是至关重要的。本文将介绍基因组序列片段压缩的相关方法，以及无性二倍体谱系重建的新问题和解决方案。 #### 基因组序列片段压缩 ##### 马尔可夫编码压缩质量值 在基因组数据中，质量值（Q - values）的压缩是一个关键问题。通过经验数据集D，我们可以估计马尔可夫链中状态转移的概率Pr(q → q′)： \[ Pr(q \to q') = \begin{cases} 0 & \text{if } q' = 0 \\ \text{reads with initial quality } q' \text{ fraction} & \text{if } q = 0 \\ \frac{\#\text{pairs } (q, q') \text{ in } D}{\#\text{occurrences of } q \text{ in } D} & \text{otherwise} \end{cases} \] 假设片段长度为L，马尔可夫分布的熵H(M)为： \[ H(M) = \frac{1}{L}H(P_1) + \frac{L - 1}{L} \sum_{q,q' \neq 0} Pr(q \to q') \log \left(\frac{1}{Pr(q \to q')}\right) \] 经验计算表明，熵约为3.3比特。为了匹配这个熵，我们使用自定义的马尔可夫编码（Markov - encoding）方案，每个转移q → q′使用 - log(Pr(q → q′))比特的霍夫曼编码。 | 压缩方式 | 每字符比特数 | 压缩比 | | ---- | ---- | ---- | | 原始文件 | 8 | 1 | | bzip2 | 3.8 | 2.11 | | Δ（未明确） | 4.25 | 1.88 | | 马尔可夫（霍夫曼） | 3.45 | 2.32 | 从表中可以看出，马尔可夫编码方案提供了2.32倍的压缩，每字符需要3.45比特。进一步使用bzip2压缩并没有改善这个结果。 ##### 质量值的有损压缩 虽然质量值的压缩仍然是一个有吸引力的研究问题，但即使技术不断改进，Q - 值压缩可能仍然是片段级压缩的瓶颈。我们考虑是否可以丢弃质量值，但目前许多下游应用，如比对、变异检测等，都将Q - 值视为推理的关键部分，因此不能简单丢弃。 我们提出了一个不同的问题：下游应用是否对Q - 值的小变化具有鲁棒性？如果是，那么“有损编码”可能对下游应用无关紧要。 我们定义了有损Q - 分数编码： 1. 首先，记不同质量值的数量为|Q| = Qmax - Qmin，用log₂(|Q|)比特编码。 2. Q - 分数计算本身就涉及精度损失，我们可以通过随机舍入（randomized rounding）来减少误差。随机舍入定义为： \[ r_{rand}(x) = \begin{cases} \lceil x \rceil & \text{with probability } x - \lfloor x \rfloor \\ \lfloor x \rfloor & \text{otherwise} \end{cases} \] 3. 对于参数b，有损Q - 分数编码为： \[ LQ - score_b(Q) = r_{rand}\left(\frac{Q - score \cdot 2^b}{|Q|}\right) \] 我们使用马尔可夫编码对2^b个不同的值进行编码。下游应用将看到Qmin + |Q| · LQ - scoreb(Q)而不是原始值Q。 我们测试了有损方案对Illumina的CASAVA应用的影响。在GAHUM的Chr 2的50M宽区域上运行CASAVA，原始Q - 值返回了17,021个符合内部阈值的变异集S。对于参数b ∈ {1, …, 5}，我们用有损编码后的分数重新运行CASAVA，得到变异集Sb。结果令人惊讶，即使b = 1（用1比特编码Q值），S中98.6%的变异调用是一致的；b = 3时，一致性提高到99.4%。而且，大部分不一致的SNP的等位基因质量接近阈值，85%的不一致SNP的等位基因质量 ≤ 10。 ##### 无损方案是否总是更好？ 我们考虑了Chr 2中38个有损（3比特）压