前几天，一位小朋友问了我这个问题：为什么没有亲缘关系的样本，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，使得样本重复测序；
除此之外，我想不到还有什么理由，让已知没有亲缘关系的样本变成同卵双胞胎/重复样本。
各位有些相关经验的，欢迎找我讨论，我想听听别的声音。
网友评论