血浆MSI（bMSI）算法介绍

2021-05-14 本文已影响0人名本无名

前言

在这之前，我们介绍了一个在组织中常用的 MSI 检测方法，如 MSIsensor、MANTIS 等

今天，我们要介绍几种应用于血浆 cfDNA 样本的 MSI 状态预测算法

1. bMSISEA

首先，我们介绍由燃石公司开发算法 bMSISEA，算法流程如下

流程图

算法分为两部分：

构建 baseline
预测样本 MSI 状态

构建 baseline

1. 计算分布

该算法从 bam 文件出发，先使用 MSIsensor 算法识别出 MS 位点的 reads 分布情况，得到的文件如下

chr1 16248728 ACCTC 11[T] AAAGG
N: 0 0 0 0 0 0 0 0 2 12 313 38 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
T: 0 0 0 0 0 0 0 0 5 378 2097 133 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

其中，第一行表示的是位点信息，前三个值分别代表染色体号、起始位置和终止位置，两个碱基序列为该位点的两端侧翼序列，中间的 11[T] 表示该为 MS 点有 11 个 T 碱基

第二行为配对样本的 reads 分布，其中数值代表该位点长度分别为 1-100 个重复 T 的 reads 数目，例如 313 代表该位点上有 313 条 reads 为 11 个 T 重复

第三行为肿瘤样本的 reads 分布

2. 识别覆盖模式

针对上述的输出文件，可以着手构建 MSS 和 MSI-H 的 reads 分布模式

MSS 分布模式

使用 MSS（WBC）样本作为训练集，以分布的峰值（如，上述示例的峰值为 2097，对应为 11 个 repeat）作为起始，向两边扩展，使得在该区域内的 reads 占比超过 75%。同时该区间能够在 25% 的 MSS 样本中重现

也就是计算一段最小的 repeat 区间，使得落在该区间内的 reads 数占该位点总的 reads 数的比例超过 75%，并以该 repeat 长度区间作为 MSS 分布模式

MSI-H 分布模式

首先，使用 MSI-H 的组织样本，由于存在样本肿瘤纯度的问题，假设肿瘤纯度为 $\color{red}u$ ，那么样本中属于肿瘤的 reads 数目 $\color{red}Y_{tumor}$ 满足如下等式
$\bbox[#8dd3c7, 15px]{ Y_{tumor} = Y_s - (1 - u) * X_s }$

其中， $\color{red}Y_s$ 和 $\color{red}X_s$ 为肿瘤和配对正常样本中的 reads 总数。

需要使用肿瘤纯度 $\color{red}u$ 来标准化肿瘤样本的 reads 总数，如何计算纯度值呢？

他们提出一个假设：肿瘤样本中落在 MSS 模式区间内的 reads 都是来自其配对正常样本的。那么就有如下公式：
$\bbox[#8dd3c7, 15px]{ Y_{mss} = (1 - u) * X }$

其中， $\color{red}{Y_{mss}}$ 为肿瘤样本中落在 MSS 模式区间内的 reads 数，X 为其配对正常样本中该区间的 reads 数。并以此估算出肿瘤纯度 $\color{red}u$ ，并使用该值去标准化 MSI-H 肿瘤样本的 reads 数，并用于后续计算。

MSI-H 的分布模式也是一段区间，满足 MSS 样本在该区间内的 reads 占比少于该位点总 reads 的 0.2%，而 MSI-H 样本在该区间的 reads 占比超过该位点总 reads 数目的 50% 以上。

目标就是为了使区间内 reads 占比在 MSS 和 MSI-H 样本中具有很高的区分度。

将满足上述要求的位点作为候选的标记位点，他们找到了 8 个位点

3. 构建 baseline

富集分析

对于每个标记位点，使用累积超几何分布，在所有 MSS 样本中对落在 MSI-H 模式区间与 MSS 模式区间上的 reads 数目进行富集分析

$\bbox[#8dd3c7, 5px]{ P(X = k) = \frac{\binom{K}{k}\binom{N-k}{n-k}}{\binom{N}{n}} }$

其中 N 为所有样本中落在 MSS 和 MSI-H 模式区间内的所有 reads 总和，K 为所有样本中 MSI-H 模式区间内的 reads 总和，相应的 N-K 为 MSS 模式区间内的 reads 总数，n、k 分别为其中某一样本中模式区间内的 reads 总数和 MSI-H 模式总 reads

在所有 MSS 样本中，使用上式计算出 P 值之后，进行负对数转换
$\bbox[#8dd3c7, 15px]{ H_s = -log(P_s(X > k_s)) }$
并将该值作为该位点的富集性指数

然后计算该位点在所有样本中富集性指数的均值和方差： $\color{red}mean_{s}$ 、 $\color{red}SD_{s}$

最后的 baseline 包含

位点位置信息
MSS 模式区间
MSI-H 模式区间
N
K
$\color{red}mean_s$ 、 $\color{red}SD_s$

阈值

对于 MSS 样本，利用 baseline 信息及超几何分布，计算出所有位点的富集性指数，并计算所有位点的富集性指数 Zscore 之和
$\bbox[#8dd3c7, 5px]{ MSscore = \sum_{s \in markers}\frac{H_s - mean_{s}}{SD{s}} }$

然后，计算 MSscore 的均值 mean 和方差 SD，并将阈值定义为
$\bbox[#8dd3c7, 15px]{ cutoff= mean + 3 * SD }$

将大于该阈值的位点定义为不稳定。

预测

预测方法与训练阈值的过程类似，也是计算出样本中，每个位点的富集性指数，然后计算所有位点的富集性指数 Zscore 之和。

最后，比较 MSscore 与 cutoff 的大小，大于阈值的样本认为是 MSI-H

2. 思路迪方法

我们介绍的第二个方法，是思路迪公司 2020 年发表论文 Wang, Zhenghang, et al.

该方法的大致流程为

流程图

该方法也是从 bam 文件开始，先计算出每个位点的不同 repeat（重复）长度对应的 reads 数目分布，具体怎么算的没有说明，应该和 MSIsensor 差不多。

可以用上面的数据分布带入来理解

算法原理

训练集包括 20 个组织 MSI-H(tMSI-H) 和 100 个组织 MSS(tMSS) 样本。

首先，他进行了一步位点过滤，未在文章中看到具体的方法，提到是根据对训练样本的不稳定敏感性排序，并取 top100

猜测应该是根据每个位点在 MSI-H 样本和 MSS 样本中的 reads 分布情况，进行了某种检验： $\color{red}\chi^2$ 检验或 T 检验，并按照显著性或统计量对位点进行排序，并取 top100 的位点进行后续分析。

如果不是，未尝不可以如此尝试

1. 计算阈值

对于每个目标位点，分别计算其

在 MSS 样本中每个 repeat 所对应的平均 reads 占比
在 MSI-H 样本中每个 repeat 所对应的平均 reads 占比

然后为每个位点计算累积平均占比，并将累积占比差异最大的 repeat 作为阈值 $\color{red}C_i$

我们举个例子来说明，假设某一位点的分布如下

其中，repeat_ 列为重复单元的长度，mss 和 msi 分别为位点在 MSS 和 MSI-H 样本中对应 repeat 的 reads 总数

mat <- data.frame(
  repeat_ = 9:13,
  mss = c(2, 12, 313, 38, 5),
  msi = c(5, 378, 2097, 133, 8)
)

> mat
  repeat_ mss  msi
1       9   2    5
2      10  12  378
3      11 313 2097
4      12  38  133
5      13   5    8

将 reads 数目转换为占比

> mat$mss <- mat$mss / sum(mat$mss)
> mat$msi <- mat$msi / sum(mat$msi)
> mat
  repeat_         mss         msi
1       9 0.005405405 0.001907669
2      10 0.032432432 0.144219763
3      11 0.845945946 0.800076307
4      12 0.102702703 0.050743991
5      13 0.013513514 0.003052270

再计算每个 repeat 的累积占比

> mat$mss <- cumsum(mat$mss)
> mat$msi <- cumsum(mat$msi)
> mat
  repeat_         mss         msi
1       9 0.005405405 0.001907669
2      10 0.037837838 0.146127432
3      11 0.883783784 0.946203739
4      12 0.986486486 0.996947730
5      13 1.000000000 1.000000000

识别出累积占比差值最大的 repeat

> mat[which.max(abs(mat$mss - mat$msi)), 1]
[1] 10

最后求得 $\color{red}C_i = 10$

定义：落在阈值之前的 reads 为不稳定 reads，即 repeat 长度小于 $\color{red}{C_i}$ 的 repeat 所对应的 reads。相对应，落在阈值之后的 reads 为稳定 reads

上述例子中，在 10 之前的所有 reads 都是不稳定的

2. 位点不稳定性

确定了每个位点的 repeat 长度阈值 $\color{red}{C_i}$ 之后。对于每个测试样本，使用二项分布模型计算每个位点的不稳定性：
$\bbox[#8dd3c7, 10px]{ P(X = n_i) = C^{n_i}_{N_i}p_i^{n_i}(1-p)^{N_i - n_i} }$

其中 $\color{red}{i}$ 表示某一测试位点， $\color{red}{p_i}$ 表示 MSS 样本中阈值 $\color{red}{C_i}$ 处的累积占比（如，上述例子 0.0378）。 $\color{red}{n_i}$ 表示不稳定的 reads 数， $\color{red}{N_i}$ 表示该位点的总 reads 数。

如果累积概率 $\color{red}{P(X \geq n_i)} \leq 0.001$ ，则认为位点是不稳定的。

如果样本中不稳定位点的比例超过 20%，则认为该样本为 MSI-H

3. msisensor-ct

最后介绍 2021 年求臻医学发表的 bMSI 检测方法 msisensor-ct

1. 数据描述

他们的研究中，真实的配对血浆样本（cfDNA 和 WBC）只有 39 例，其中有 5 例的 MSI 状态是经过 IHC 验证的，其他样本的状态是使用 MSIsensor 推断的。

在使用 MSIsensor 时，进行了过滤：

对于 MSIsensor 输出的 reads 分布文件，即最开头的例子文件，去除正常样本中 repeat 长度与参考等位不一致且 reads 占比超过 10% 的位点
将 cfDNA 和 WBC 配对样本中剩下的位点分别进行抽样，抽样的平均深度为 200X，然后使用 $\color{red}\chi^2$ 检验对肿瘤和配对样本进行检验，如果 p < 0.05 则位点不稳定，如果样本中不稳定位点超过 10% 则认为样本是不稳定的。
最后分出来 4 个 MSI 和 35 个 MSS 样本

疑问：

用这种方式来判断样本的状态作为标准，然后再使用其他方法来构建血浆样本预测模型

是不是舍近求远呢？如果不相信这种方法预测的结果，需要开发其他算法，那为何又使用这个方法的预测结果作为标准呢？

最后，样本的构成如下：

样本构成

TCGA 的 1565 个 WES 样本作为训练集，然后使用这些样本进行 cfDNA 样本的模拟，模拟的 cfDNA 纯度为 0.1%，并用作算法的测试集。

另外的 159 个 EGA 和 NGDC 模拟样本作为独立验证集

2. 算法原理

算法流程

2.1 位点过滤

筛选的位点需要满足：

测序深度 $\color{red}\geq 20X$
在 $\color{red}\geq 10\%$ 的样本中出现
不稳定分布比例 $\color{red}\geq20\%$

最后保留了 25861 个位点

2.2 位点模型

对于每个位点，分别获取其在所有肿瘤样本中的分布，随机将其分为 80% 的训练集和 20% 的验证集。

由于 cfDNA 含量很低，所以在肿瘤样本的 reads 分布中，与参考基因组相同的等位频率很容易影响低频变异，所以将该等位从分布中删除

对于某一位点，将其在所有样本中的分布合并成一个矩阵，称为 site-set。位点在样本中的状态由 MSIsensor 算法确定，即如果 $\color{red}\chi^2$ 检验肿瘤和配对正常的分布有显著差异（p < 0.05），则位点的状态为不稳定。

然后分别应用如下机器学习模型进行训练：

XGBoost
AdaBoost
GBDT
SVM
LR

构建二分类模型

2.3 状态预测

利用上一步训练出的模型，在模拟的测试集中进行预测。并计算预测结果的 AUC 值，通过限制 AUC 的值来进一步筛选位点，最后得到了 1476 个位点。

每个位点都有自己的模型，对于一个测试样本，每个位点模型都有一个预测状态，通过计算预测为不稳定的位点的比例，来判断样本的 MSI 状态。