排查疑似污染样本的方法

“变异信息的统计对质控来说是很重要的一步，可以通过此步排除疑似污染样本。”

我们对测序结果进行分析的时候，在实验和下机数据质控层面可能没有发现异常，但进行后续分析时，可能就会出现问题，比如与参考基因组比对后发现某样本各项指标异常，或变异检测注释后统计位点信息时发现某样本异常。今天要和大家分享的就是后面一种情况。

01

—

分析方法

以WES全外显子组测序人的样本为例，如何根据变异注释结果分析出样本是疑似污染样本呢？方法如下：

snp数量

indel数量

het ratio

patho/likey patho的数量

rediscovery rate

index拆分

实验

12例样本采用WES测序技术，捕获panel为Aglient 64M。按照上述的分析标准，结果分别如下：

表 1. 疑似污染样本与11例样本的变异注释统计结果

注：J表示疑似污染样本sample J，11 sample mean表示其他11例样本每项指标的均值。NFE：Non-Finnish European

通过上面的结果可以看出，样本J在突变数量、杂合率及疑似致病位点数量上，远高于其他11例样本的均值。其中，patho/likely patho是根据ACMG分级的结果。而在正常人群突变频率比较时发现，样本J在东亚人数据库的比例远低于其他11例样本，而位点在非洲人群的比例高于其他11例样本。

排除Index拆分错误及实验方面的原因，因此判定样本J为疑似污染样本。其中所有12例样本为单Index，按照不容错拆分。

重点讲下rediscovery rate是怎么计算的：

rediscovery rate位点重现率是基于1000g、exac、gnomAD、berry400k这4个数据库进行分析的。由于ANNOVAR注释的时候只用了1000g_all和1000g_eas，所以统计结果里只有这2个数据库。而exac和gnomAD是有多个人种数据库的。

千人基因组1000g是全基因组测序，平均测序深度7.4X，样本数量2504。其实每例样本也做了WES，平均测序深度65X。但为了便于计算，位点在数据库重现率计算上，1000g和berry400k直接统计位点在数据库频率大于0的位点。而exac数据库是基于WES测序的，gnomAD整合了exac、1000g和ESP，所以需要先统计有多少位点在exac的exome variants calling interval list里出现过，再统计有多少位点在exac里有记录，即频率大于0。

先将annovar注释结果转为vcf，此步自己写程序处理即可。转成vcf的目的是下一步，应用intersectBed进行与exac捕获Panel overlap位点的查找，因为intersectbed的输入文件只能是vcf或bed。

命令行：

为什么不用exac exome panel与我们Aglient 64M取交集，统计落在交集里的位点呢？因为目的是先看有多少位点在exac panel范围里，再看有多少位点是在exac数据库有频率记录的位点。如果看exac没有覆盖到的位点，可以用两者取交集的方法去做。

补充知识：

在计算rediscovery rate时，分子是在数据库有频率记录的位点，分母是该样本有多少个位点。

突变频率 frequency = allele count / allele number

如果在数据库中频率为0，是说明该位点的allele count为0，但这个位点是存在的。

bed文件格式：

column 1: chromosome

column 2: start position (0-base)

column 3: end position (1-base)

02

—

解决办法

如果确定这个样本是疑似污染的样本，那么下一步该如何分析呢？

方法一：

剔除疑似污染样本，继续后续分析。

方法二：

重新送样，新样本，相同表征的不同患者。

通过排查疑似污染样本，学习了变异质量控制方法。希望大家有所收获。

大家有什么疑问，欢迎随时交流，后台留言就可以哟

我是yyt，分享基因知识、研究心得、职场感悟、碎碎念等

由于2018年3月注册的公号没有留言功能，欢迎后台留言，我会及时回复哟

加油~↖(^ω^)↗~