前几天，一位小朋友问了我这个问题：为什么没有亲缘关系的样本，IBD显示他们是同卵双胞胎或者重复样本。

具体来说，使用PLINK的--genome参数计算后得到的PI_HAT（Proportion IBD）全是1。

如下所示：

image

鉴于这位小朋友的微信ID特别有味道（屎尿屁之类），而且赞赏过我文章，让我印象很深刻。

于是，我决定亲自操刀，让他发测试数据给我。

测试数据是ped/map格式。
map如下图所示，可以看到，RS号没有统一好（第二列）：

image

好人做到底，给人家统一一下RS号，统一好后如下所示：

image

神清气爽了！

做数据分析，清洗数据很重要，太脏的数据不仅影响工作效率，还影响结果。

比如本推文出现的问题，在没有看PLINK源代码的情况下，我们是不知道是根据位置和染色体信息，还是根据RS号信息计算IBD。

如果是根据位置和染色体信息，那么只需要确保这两个信息准确就行了；

但如果是根据RS号，没有统一好RS号的话，会丢失掉很多位点，影响结果。

数据清洗后，开始计算IBD：

plink --bfile file --indep-pairwise 50 5 0.2 --out file_indep #Pruning
plink --bfile file --extract file_indep.prune.in --genome --out file_indep.prune.in.ibd #计算亲缘关系
awk '$10>=0.95' file_indep.prune.in.ibd.genome #提取PI_HAT大于0.95的样本

清洗后的结果还是跟之前一样，本无亲缘关系的样本还是有亲缘关系，如下图红框所示：

image

到这一步至少确定了，PLINK计算IBD是根据位置和染色体信息，不需要统一RS号。

但我们还是无法找到问题所在。

想确认样本间是不是同卵双胞胎/重复样本，最万无一失的方法是计算样本间碱基的一致性（kappa值）。但我懒得写脚本。

于是我使用了一种偷懒的办法：把样本拷贝后更换ID变成新的样本，再计算亲缘关系。如下：

#拷贝样本
cp file.bed dd.bed
cp file.bim dd.bim
cp file.fam dd.fam

随后，修改样本ID。

原始file.fam的ID如下所示：

image

修改dd.fam样本ID变成新的样本：

image

实际上，dd.fam的63547_63547样本就是file.fam的63547；
同理，dd.fam的63559_63559样本是file.fam的63559；
更改ID是为了合并；

更改ID后，合并样本，计算IBD:

plink --bfile file --bmerge dd.bed dd.bim dd.fam --make-bed --out merge #合并样本
plink --bfile merge --indep-pairwise 50 5 0.2 --out merge_indep #Pruning
plink --bfile merge --extract merge_indep.prune.in --genome --out merge_indep.prune.in.ibd #计算亲缘关系
awk '$10>=0.95' merge_indep.prune.in.ibd.genome #提取PI_HAT大于0.95的样本

IBD的结果如下所示：

image

毫无意外的， 样本63547和样本63547_63547之间的PI_HAT为1，样本63559 和样本63559_63559之间的PI_HAT为1。他们本来就是同一个样本，被我们拷贝过来的。

此外，我们也可以观察到样本63547和样本70973的PI_HAT为1; 样本63559 和样本69111的PI_HAT为1，与重复样本（样本63547和样本63547_63547、样本63559 和样本63559_63559）的结果完全一致。

到这里就说明了，1）要么他们（样本63547和样本70973、样本63559 和样本69111）是同卵双胞胎/重复样本，2) 要么贴错样本ID，使得样本重复测序；

除此之外，我想不到还有什么理由，让已知没有亲缘关系的样本变成同卵双胞胎/重复样本。

各位有些相关经验的，欢迎找我讨论，我想听听别的声音。