Bonaparte：贝叶斯网在灾难遇难者识别（DVI）中的应用

荷兰一家软件公司开发的Bonaparte系统已成功应用在多起灾难的遗体鉴别任务中。根据其官网给出的联系方式，联系到了其技术主管Willem Burgers，向他询问Bonaparte背后的技术原理。Willem Burgers很热情地将他的新书《Interactive Collaborative Information System》发给我。遂翻译出书中关于Bonaparte背后原理的章节。

1. 前言

大型灾难不仅出现在电影里，在现实生活中也时有发生，比如世贸大厦的恐怖袭击袭击、海啸、飞机失事等等。灾难发生后，搜索遇难者遗体并鉴定死者身份非常重要。现代DNA技术大大促进了灾难遇害者鉴别（Disaster Victim Identification, DVI）。灾难后收集到的遗体残骸往往难以鉴别身份，但法医们可以从少量的遗体残骸中采录到DNA图谱，从而为DVI提供可能。鉴定方法是将身份不明的遗体残骸与报告失踪人口的亲属进行DNA匹配，因为亲属的DNA图谱与遇难者DNA是相关的，一级家庭成员的DNA有50%是相同的。
在仅有一名遇难者的案件中，DVI对法医而言是一项简单直接的任务。但随着遇难者人数的增加，DVI就变得难以手工完成，必须借助自动化程序来实现。
贝叶斯网非常适合于对家族谱系中的DNA图谱进行统计相关性建模。贝叶斯网的一个优势是它非常灵活，能适应不同的家族谱系类型，且允许合并其它因素发挥作用，比如测量误差概率、数据缺失、更高级的基因标记等。

图1. 匹配的问题。根据受害者和失踪者亲属的DNA图谱，将身份不明的受害者(右图蓝色)与报告的失踪人员(左图红色)进行比对。方代表男性，圆代表女性，实线表示可以获得DNA图谱，虚线表示不可获得。

目前，有一款叫Bonaparte的荷兰软件可用于DVI。Bonaparte的计算引擎使用了自动生成的贝叶斯网和贝叶斯推理方法，能够基于DNA图谱结合家族谱系信息进行亲属关系分析。其支持在涉及数百名遇难者的大规模数据集上进行推理。
本文将介绍Bonaparte所采用的贝叶斯网方法。首先将简要介绍DNA图谱，模型中将计算DNA图谱在不同家族谱系中的似然率。

2. 两个假设的似然率之比

假定我们已经获得一个遇难者(Missing Person, MP)的家族谱系。在该谱系中，有些家族成员(Family members, FAM)提供了DNA样本。此外，还获得了一个待鉴别身份者(Unidentified Individual, UI)的DNA。我们的问题就是UI是否等于MP？为了判断该问题，我们假设有一个关于DNA图谱的概率模型P，其中以家族成员的DNA信息作为证据变量。同时，我们提出两个假设$H_p$和$H_d$：

• $H_p$：假设MP=UI，则MP是已观测到的，证据变量为 E = { D N A M P + D N A F A M } E=\{DNA_{MP}+DNA_{FAM}\} E={DNAMP​+DNAFAM​}
• $H_d$：假设UI是一个不相关的人U，则U是已观测到的，MP是未观测到的，证据变量为 E = { D N A U + D N A F A M } E=\{DNA_U+DNA_{FAM}\} E={DNAU​+DNAFAM​}

在概率模型P中，这两个假设的似然率之比为：

$$LR=\frac{P(E|H_p)}{P(E|H_d)}$$

若已知先验几率$P(H_p)/P(H_d)$，则后验几率$P(H_p|E)/P(H_d|E)$可通过先验几率与似然率之比相乘得到：

$$\frac{P(H_p|E)}{P(H_d|E)}=\frac{P(E|H_p)P(H_p)}{P(E|H_d)P(H_d)}$$

3. DNA图谱

本节将对基于DNA的亲缘关系检测原理进行简要介绍。人类的DNA分布于细胞核中的染色体上，一个正常的人类细胞有46条（23对）染色体，每对染色体中，一份遗传自父亲，另一份遗传自母亲。其中有22对常染色体，它们有相同的长度，通常包含相同的基因（DNA功能单元）；另外还有一对性染色体。男性为XY染色体，女性有两条X染色体。
任意两人的DNA有99%以上是相同的，因此大多数DNA对于身份鉴定是无用的。然而，染色体上有一些特定的点位，在这些点位上的DNA存在个体间差异，这种差异被称为基因标记。在基因学上，这些特定点位被称为基因座(复数形式loci,单数形式locus)。
在法医研究中，短串联重复序列(Short Tandem Repeat, STR)基因座是用得最多的。STR也称微卫星DNA(microsatellite DNA)，通常是基因组中由1~6个碱基单元组成的一段DNA重复序列，由于核心单位重复数目在个体间呈高度变异性并且数量丰富，构成了STR基因座的遗传多态性。一般认为人类基因组平均每15kb就存在一个STR基因座。其多态性成为法医物证检验个人识别和亲子鉴定的丰富来源。不同人体基因组卫星DNA重复单位的数目是可变的，因此，形成了极其复杂的等位基因片段长度多态性。美国FBI利用13对STR基因座加性染色体位点制作一个“联合DNA索引系统(CODIS)”的全国性DNA数据库，应用到犯罪鉴定与法医学鉴定。

图2. 联合DNA索引系统

STR基因座本质上是一类变异，表现为两个或两个以上碱基对重复发生，如下：

$$(CATG{)}_3=CATGCATGCATG$$

碱基对的重复次数$x$（本例中$x=3$）是存在个体差异的。有时也存在不完整重复，比如$CATGCATGCATGCA$，此时的重复次数$x=3.2$。对于法医鉴定中所使用的STR基因座，$x$的值都能精确建档。
这些STR基因座的集合构成了DNA图谱。由于染色体是成对存在的，所以一个图谱由成对的STR基因座组成。比如CODIS所使用的13对基因座可记作：

$$\bar{x}={(}^1{x}^1{,}^1{x}^2),{(}^2{x}^1{,}^2{x}^2),\cdots ,{(}^{13}{x}^1{,}^{13}{x}^2)$$

其中${}^{\mu }{x}^s$表示在基因座$\mu$处的重复次数。由于染色体成对存在，因此每一个基因座都有一对等位基因序列${}^{\mu }{x}^1$和${}^{\mu }{x}^2$，其中一个继承自父亲，另一个继承自母亲。然而，目前的DNA分析技术还无法鉴别等位基因来自父母哪一方。因此${(}^{\mu }{x}^1{,}^{\mu }{x}^2)$与${(}^{\mu }{x}^2{,}^{\mu }{x}^1)$无法区分。为了数学表达的一致性，我们在记录时遵守${}^{\mu }{x}^1\le$${}^{\mu }{x}^2$的规范。

图3. 一个最基本的亲缘分析示例，包含父亲、母亲和子女。方框代表男性，圆圈代表女性。右图为对应的贝叶斯网络，灰色节点为观察到的节点，上标p表示来自父亲的等位基因，上标m表示来自母亲的等位基因。

染色体遗传自父母，父母双方从每对染色体中等概率地遗传其中一条给孩子。且在染色体遗传的过程中，等位基因发生突变的可能性只有约为0.1%。
最后，在DNA分析中，有时会出现某个点位的等位基因丢失，如观察值为${(}^{\mu }{x}^1,?)$，其中?表示一个不确定项。

4. 基于贝叶斯网的亲缘分析

在本小节，我们将描述如何为DNA图谱建立贝叶斯网络。首先，各STR基因座的遗传以及我们对它们的观察都是独立的，因此对每一个基因座，我们可以建立一个独立的概率模型$P_{\mu }$。在下文，我们只考虑对一个基因座进行建模，因此为了记录简便，我们将基因座的位置$\mu$省略。

4.1 等位基因遗传模型

在一个族谱中，除了家族的创始者，其他任何一个个体$i$都有一对父母，父亲记为$f(i)$，母亲记为$m(i)$。由于交配后，等位基因会发生交换重组，因此DNA图谱和亲属的等位基因之间的统计相关性可以借助于该族谱进行计算。对于一个给定的基因座，个体$i$拥有一个来自父亲的等位基因$x_i^f$和一个来自母亲的等位基因$x_i^m$，这对等位基因记作$x_i=(x_i^f,x_i^m)$。我们有时也用$x_i^s$来记这对等位基因，其中上标$s$的取值为 { f , m } \{f,m\} {f,m}，$x_i$有$N$种取值，$N$表示基因座上的等位基因数。

在族谱中，来自创始者的等位基因有两条，它们彼此独立同分布，我们假定这个分布是给定的，为在全体人口中的分布，记为$P(a)$。当然，在更高级的模型中，我们可以进一步考虑这对来自创始者的等位基因存在相关性，且不同分布，比如创始者来自某个特定的子人口集，从而可提高建模的精度。但是，做这种相关性假设后，将会导致模型消耗巨大的内存和计算资源。因此，本专题我们仅在独立同分布的假设下进行考虑。

若个体$i$在族谱中有父母，则他的等位基因的分布满足如下关系：

$$P(x_i|x_{f(i)},x_{m(i)})=P(x_i^f|x_{f(i)})P(x_i^m|x_{m(i)})$$

其中：

$$\begin{gathered} P(x_i^f|x_{f(i)})=\frac{1}{2}\sum _{s=f,m}P(x_i^f|x_{f(i)}^s) \\ P(x_i^m|x_{m(i)})=\frac{1}{2}\sum _{s=f,m}P(x_i^m|x_{m(i)}^s) \end{gathered}$$

该公式表明，个体从其父母处各获得两条等位基因中的一条是等概率事件。其中条件概率$P(x_i^f|x_{f(i)}^s)$和$P(x_i^m|x_{m(i)}^s)$是由基因突变模型$P(a|b)$给出的，表明遗传自父辈的等位基因是$b$，在子女处变异为$a$的概率。对于STR基因座中等位基因的突变机制尚不明确，但有证据表明从父亲到子女的突变率是从母亲到子女突变率的10倍。我们假定在族谱中的所有个体的性别是已知的，但为了记录简便，我们忽略父母的性别因素对基因传递过程中突变率的影响，也忽略基因位于不同基因座位置时的突变率差异。

基因突变率最简单的模型是不考虑突变。这种假设有一定合理性，因为突变的发生率很低，我们可以针对个案单独考虑突变模型。这一假设可以极大提高概率推理的算法效率。然而，一旦存在哪怕一个微小突变，都会导致基因匹配百分百被拒绝，即使其它所有基因座都能完美匹配。基因突变模型对于模型的鲁棒性很重要。我们可以设计一个简化的基因突变模型，比如平均分布模型：

$$P(a|b)=\{\begin{array}{ll}1-\mu & a=b\\ \mu /(N-1)& a\ne b\end{array}$$

其中$\mu$是基因突变率（不要与前文的基因座序号$\mu$混淆)。该突变率$\mu$包含了基因座位置、性别等所有相关影响因素。

相比于忽略基因突变的模型，该基因突变模型对算法性能影响不大。且由于应用了平均分布的突变模型，基因的概率分布也趋于平坦。有人会说这与基因突变产生的生物多样性现象一致。然而，回忆我们对族谱中创始者等位基因分布的假设，如果族谱中存在多个不同代的创始者，他们的等位基因被赋予了相同的分布，这与此处的基因突变模型存在差异，因此随着族谱中增加更多关联祖先，会导致似然率发生微小偏差。

4.2 观察模型

在族谱上，我们可观察到个体$i$的某对等位基因，记为${\bar{x}}_i$，省略下标则不指定个体。我们无法观察到等位基因来自父母哪一方，因此$x^f=a,x^m=b$与$f^f=b,x^m=a$是等价观察，记作$\bar{x}=(a,b)$，其中约定$a\le b$。当发生基因缺失时，我们用？来代替缺失的基因，如$\bar{x}=(x,?)$。我们用$L$来对基因缺失事件进行建模：$L=1$表示缺失了一个等位基因，如$\bar{x}=(x,?)$；$L=0$表示不缺失，如$\bar{x}=(a,b)$；当$L=2$时，我们实际上没有观察到任何信息，因此不对这种情况建模。从而，可对变量建立如下观察模型：

当$L=0$时：

$$P(\bar{x}|(a,b),L=0)=\{\begin{array}{l}1if\bar{x}=(a,b)\\ 0otherwise\end{array}$$

当$L=1,a\ne b$时：

$$\{\begin{array}{l}P(\bar{x}=(a,?)|(a,b),L=1)=\frac{1}{2}\\ P(\bar{x}=(b,?)|(a,b),L=1)=\frac{1}{2}\end{array}$$

当$L=1,a=b$时：

$$P(\bar{x}=(a,?)|(a,a),L=1)=1$$

5. 推断

将族谱上创始者的等位基因先验概率模型、基因遗传模型和观察模型相乘，我们就可以得到一个贝叶斯网，其中包含了等位基因$x$和族谱上已观察到的等位基因$\bar{x}$。假定族谱上包含的个体集为$I=1,2,\cdots ,K$，其中的创始者组成子集$F$，且假定各个体的等位基因缺失情况$L_j$已知，则该贝叶斯网的联合概率分布如下：

P ( { x ˉ , x } I ) = ∏ j P ( x ˉ j ∣ x j , L j ) ∏ i ∈ I \ F P ( x i ∣ x f ( i ) , x m ( i ) ) ∏ i ∈ F P ( x i ) P(\{\bar x,x\}_I)=\prod_jP(\bar x_j|x_j,L_j)\prod_{i\in I\backslash F}P(x_i|x_{f(i)},x_{m(i)})\prod_{i\in F}P(x_i) P({xˉ,x}I​)=j∏​P(xˉj​∣xj​,Lj​)i∈I\F∏​P(xi​∣xf(i)​,xm(i)​)i∈F∏​P(xi​)

在该模型下，对于一组DNA图谱，我们可以计算其似然率。族谱上部分个体的等位基因被采集和观察到，这些个体组成的子集记为$O$，其似然率为边缘概率分布记为 P ( { x ˉ } O ) P(\{\bar x\}_O) P({xˉ}O​)：

P ( { x ˉ } O ) = ∑ x 1 ⋯ ∑ x K ∏ j ∈ O P ( x ˉ j ∣ x j , L j ) ∏ i ∈ I \ F P ( x i ∣ x f ( i ) , x m ( i ) ) ∏ i ∈ F P ( x i ) P(\{\bar x\}_O)=\sum_{x_1}\cdots\sum_{x_K}\prod_{j\in O}P(\bar x_j|x_j,L_j)\prod_{i\in I\backslash F}P(x_i|x_{f(i)},x_{m(i)})\prod_{i\in F}P(x_i) P({xˉ}O​)=x1​∑​⋯xK​∑​j∈O∏​P(xˉj​∣xj​,Lj​)i∈I\F∏​P(xi​∣xf(i)​,xm(i)​)i∈F∏​P(xi​)

在实际应用中，由于等位基因的取值空间非常庞大，这将会使得“联合树"算法变得不可行。幸运的是，我们可以通过一个叫做“值约简"的方法大大降低算法的复杂度。假设我们所观察到的等位基因值为等位基因可能取值空间的一个子集，记为$A$，$||A||=M$，包含M个不同的取值点。我们将等位基因的所有其它可能取值都约简为$z$，从而将等位基因的取值可能数约简为$M+1$。即对于 a ∈ A , L ∈ { 0 , 1 } ) a\in A,L\in\{0,1\}) a∈A,L∈{0,1})，且$b_1,b_2,b_3,b_4\notin A$，有如下等式成立：

$$\begin{array}{rl}P(a|b_1)& =P(a|b_2)\\ P(\bar{x}|(a,b_1),L)& =P(\bar{x}|(a,b_2),L)\\ P(\bar{x}|(b_1,a),L)& =P(\bar{x}|(b_2,a),L)\\ P(\bar{x}|(b_1,b_2),L)& =P(\bar{x}|(b_3,b_4),L)\end{array}$$

从而，我们可以用$P(a|z)$代替$P(a|b)$，用$P(\bar{x}|(a,z))$代替$P(\bar{x}|(a,b))$，等等。其中$z$的条件概率可通过下式计算：

$$P(z|x)=1-\sum _{a\in A}P(a|x)$$

其中$x\in A\cup z$.

最终，我们的推断过程为：首先对数据进行“值约简"，然后对每一个基因座调用“联合树"算法计算似然率和似然几率。

6. 波拿巴软件

波拿巴软件能处理大规模基因匹配问题，该软件为CS架构，其计算核心和内部数据库运行在服务器上，所有匹配结果也存储在内部服务器上。服务器通过XML和https接口与外界相连。用户可通过web浏览器登录系统，因此用户本地计算机上无需安装任何软件。大家可访问www.dnadvi.nl来试用它的演示版本。

图3. Bonaparte基本架构

7. 总结

波拿巴软件是基于贝叶斯网的遇难者身份鉴别系统。在该系统中，贝叶斯网被用来对族谱中不同个体的DNA图谱建立统计关系模型。通过贝叶斯推理，可以为法医提供接受假设与拒绝假设这两个假设的后验几率之比。贝叶斯网中各变量之间的概率关系是建立在遗传学第一定律“基因分离定律"上的。该系统的一个特征是能从数据中自动、实时地迭代生成模型。