• 测序基础
  • DNA 基础
  • PE 测序
• Barcode
  • 单端 barcode
  • 双端 barcode
  • 单读段 barcode
• 链偏 （Strand Bias）
  • 定义
• 概率模型
  • 马尔可夫链

测序基础

DNA 基础

• DNA 由两条反向互补的链组成，每条链都有一个 5’ 端和一个 3’ 端。
• 标准的 DNA 序列方向是从 5’ 到 3’。

PE 测序

• PE (Pair Ended) 测序可以提供片段两端的序列信息。
• 这种测序方式可以帮助我们了解片段的方向和距离信息。

假设测序的 DNA 片段是：

5' - ATGCACTGACTGAGCTA - 3'
3' - TACGTGACTGACTCGAT - 5'

测序结果为： 读段 1：5' - ATGCACTGA - 3' 读段 2：5' - TAGCTCAGT - 3'

###

Barcode

目的：在测序中，我们经常一次性测序来自不同样品的多个 DNA 片段。为了在测序完成后能区分各个样品，我们会在每个片段上加上一个独特的 DNA 序列，称为 barcode。测序后，我们可以根据这些 barcode 来拆分、分类数据。

单端 barcode

只在片段的一个端（如 5’ 端）加入 barcode。示例：

加入 barcode GGTT 后的片段：

5' - GGTTATGCACTGACTGAGCTA - 3'
3' - AACCTACGTGACTGACTCGAT - 5'

• 测序读段 1：5' - GGTTATGCA - 3'

双端 barcode

在片段的两端都加入 barcode。示例：

5’ 端的 barcode 为 GGTT，3’ 端的 barcode 为 AACC。

5' - GGTTATGCACTGACTGAGCTAAACC - 3'
3' - TTGGTACGTGACTGACTCGATGGTT - 5'

• 测序读段 1：5' - GGTTATGCA - 3'
• 测序读段 2：5' - GGTTATCGA - 3'

单读段 barcode

在 paired-end 测序中，只在一个特定的读段（如 Read1）添加 barcode。示例：

在 Read1 上添加 barcode GGTT：

5' - GGTTATGCACTGACTGAGCTA - 3'
3' - AACCTACGTGACTGACTCGAT - 5'

• 测序读段 1：5' - GGTTATGCA - 3'

链偏 （Strand Bias）

定义

概率模型

马尔可夫链

通过 Markov Chain 我们可以精确计算出，在观测到 \( n \) 个碱基时，出现 \( m \) 个特定类别（ACGT）的碱基 的概率。下面，我们以碱基类别 A 举例。

我们用一个 \( n+1 \) 维的向量¹，表示这个概率分布 \( \mathbf{x} \in \mathbb{R}^{n+1} \)。

1. 没有观察任何碱基的时候，系统状态是 \( \mathbf{x}_0 = [1, 0, …]^\top \)， 表示 A 必然出现 0 次。
2. 观察 1 个碱基后，我们有 \( p \) 的概率观察到 A，系统状态是 \( \mathbf{x}_1 = [q, p, 0, …]^\top \)； 其中，\( q = 1 - p \)，表示碱基没有出现的概率。
3. 当我们观察第 2 个碱基的时候，我们有以概率 \( p \) 多观察到一个 A，系统状态是 \( \mathbf{x}_2 = [qq, pq, pp, 0, …]^\top \)；

我们可以把这个过程抽象为一个 Markov Chain，其转移矩阵为，

从初始状态 \( \mathbf{x}_0 \) 开始，我们通过观测 \( n \) 次结果，得到最终的状态 \( \mathbf{x}_n \)，并且

我们可以用下面的代码，来计算最终状态，并且通过蒙特卡洛模拟，验证结果的准确性，

import numpy as np


def exact(p, n):
    M = np.eye(n + 1) * (1 - p)
    np.fill_diagonal(M[1:], p)
    M[-1, -1] = 1
    s = np.zeros((n + 1, 1))
    s[0] = 1
    for _ in range(n):
        s = M @ s
    return s


def monte_carlo(p, n, n_sample=int(1e5)):
    trajs = np.zeros(n_sample, dtype=int)
    for _ in range(n):
        selected = np.random.permutation(n_sample)[:int(p * n_sample)]
        trajs[selected] += 1
    count = np.bincount(trajs, minlength=n+1)
    return count / count.sum()


if __name__ == "__main__":
    from time import time
    import matplotlib.pyplot as plt

    p, n = 0.2, 10

    t0 = time()
    solution = exact(p, n)
    print(f"Exact calculation: {(time() - t0) * 1000:.2f} ms")

    t0 = time()
    simulations = []
    for _ in range(50):
        simulations.append(monte_carlo(p, n))
    print(f"Monte Carlo Simulation: {(time() - t0) * 1000:.2f} ms")

    plt.errorbar(
        x=np.arange(n+1),
        y=np.mean(simulations, axis=0),
        yerr=np.std(simulations, axis=0, ddof=1),
        label="MC simulation",
        mfc='none', color='tomato', marker='o', ls='none', capsize=2
    )
    plt.plot(solution, color='teal', zorder=2, marker='x', label="Analytical")
    plt.xlabel("Occurance", fontsize=16)
    plt.ylabel("PDF", fontsize=16)
    plt.yscale('log')
    plt.legend(ncol=2, fontsize=16, handlelength=1)
    plt.tight_layout()
    plt.show()

运行代码后，得到的结果如下

Exact calculation:      0.07 ms
Monte Carlo Simulation: 1170.33 ms

计算解析解的速度是蒙特卡洛模拟的 16714 倍²！