2022-11-09 19:38:20 星期三
目的
平时需要核查二代 DNA 数据检出的胚系变异是否可信，需要反复的打开 igv ，核查变异位置及相关条件，才能判断变异真假程度。此过程重复聒噪，当 igv 加载缓慢时尤为费时，而且人的判断可能带有主观意识，不能保证每次判断结果完全一致，特别面对不太确定变异。
思路
搜集评定变异可信度所需的证据项，给与每个证据项符合逻辑规律的加减分机制，最终得分越高，变异可信度越高，并划分合理的等级。（仿照 ACMG 的 5 等级）
证据项
支持性证据
- 等位基因深度（AD），测到变异的 reads 越多，变异真实性越高。因此该项赋分应该是正数，应符合一个单调递增函数，且 AD 较小时（AD<=8) 加分比重应该较小，即斜率与 AD 成正比。
- 等位基因频率（AF），测到变异的 reads 占总体比重越高，变异真实性越高。赋分应是正数，需符合单调递增函数，且 AF 越小时加分比重越小，即斜率与 AF 成正比。
- 测穿 reads 的对数，测穿变异的 reads 对数越多，变异真实性越高。赋值应该是正数，需符合单调递增函数，且其数量越高，加分比比重越小，即斜率与测穿 reads 的对数成反比。
反对性证据
- 总的测序深度（DP），该项赋分应该为负数，在 DP 越小，变异越不可靠。DP 应符合一个单调递减函数，且 DP 越小时，比重越高（斜率越大）。
- 是否是 reads 的唯一比对位置？唯一比对区域就不用考虑是比对造成的假变异，非唯一比对位置应该减分。
- 是否是复杂序列区域？变异处于简单序列区域（如 A 碱基重复区域）可能是测序错误导致的假变异，应该减分。
- 是否处于高质量 reads 上？低 MQ reads 上的变异较可能是假变异，应该减分。
代码示例
点击查看代码
import sys
import math
import logging
logging.basicConfig(
level = logging.DEBUG,
format = "[%(asctime)s - %(levelname)-4s - %(module)s : line %(lineno)s] %(message)s",
datefmt = "%Y-%m-%d %H:%M:%S"
)
def DP_score(DP: float) -> float:
""" DP 深度得分 """
if DP <= 6:
return round(-((DP-7) ** 2), 2)
return 0
def AD_score(AD: int) -> float:
""" AD 等位基因深度的得分 """
if AD <= 8:
return round(AD**2/2)
return 40
def AF_score(AF: float) -> float:
""" AF 等位基因频率比率得分, 越大得分越多 """
if AF < 1:
return round(AF**2 * 50, 2)
return 50
def Is_score(Is_reads: int) -> int:
""" 测穿 reads 对数量得分 """
if Is_reads <= 5:
return math.sqrt(Is_reads) * 17
return 45
def cal_level(score: int) -> str:
""" 根据得分划分等级 """
if score < 0:
return 'E'
elif 0 <= score < 40:
return 'D'
elif 40 <= score < 80:
return 'C'
elif 80 <= score < 120:
return 'B'
return 'A'
def call_score(
AD: int,
AF: float,
Is_reads: int,
Uniq_region: bool,
complex_region: bool,
high_MQ_reads: bool,
) -> float:
"""
Is_reads : 测穿的 reads 对数
Uniq_region : 是否为唯一比对区域
complex_region : 是否为高复杂度区域
high_MQ_reads : 是否在高 MQ 的 reads 上，变异所在 reads 的平均 MQ
不低于变异位置所有 reads 的平均 MQ 的 50%
"""
# 数据检查
if AF < 0 or AF > 1:
logging.error(f'AF 不在 (0,1] 范围内！')
sys.exit(1)
if Is_reads*2 > AD:
logging.error(f'测穿变异的 reads 数量比 AD 要多:\
\n\tAD:{AD} AF:{AF} 测穿 reads 对数:{Is_reads}')
sys.exit(1)
# 总深度 DP
score = 0
score += DP_score(AD/AF)
score += AD_score(AD)
score += AF_score(AF)
score += Is_score(Is_reads)
if not Uniq_region:
score -= 30
if not complex_region:
score -= 30
if not high_MQ_reads:
score -= 30
# 三相不满足，多减，两项不满足，多减
if not Uniq_region and not high_MQ_reads and not complex_region:
score -= 20
elif (not Uniq_region and not high_MQ_reads) or \
(not high_MQ_reads and not complex_region) or \
(not Uniq_region and not complex_region):
score -= 10
if Uniq_region and complex_region and high_MQ_reads:
score += 10
return (round(score, 2))
结果解读
从 cal_level 函数获得的 A-E 的评级，分别表示：
- A: 可信的变异
- B: 可能是真的变异
- C: 不确定真实性的变异
- D: 可能是假的变异
- E: 假的变异
模拟数据验证
采用造数据的形式，模拟出低深度（DP<30)时，该评分机制评级效果。（代码基本在 2s 左右运行完成，共有 20560 种可能的情况）
代码
点击查看代码
with open('res.tsv', 'w') as fo:
fo.write('Uniq\tcomplex\thigh_MQ\tDP\tAD\tAF\tIs_reads\tScore\tLevel\n')
res = ''
for DP in range(1, 30):
for Uniq_region in [0,1]:
for complex_region in [0, 1]:
for high_MQ_reads in [0, 1]:
for AD in range(1, DP+1):
if AD > DP:
continue
AF = round(AD / DP,2)
for Is_reads in range(0, (AD // 2) + 1):
score = call_score(AD, AF, Is_reads, Uniq_region, complex_region, high_MQ_reads)
level = cal_level(score)
res = f"{res}{Uniq_region}\t{complex_region}\t{high_MQ_reads}\t{DP}\t{AD}\t{AF}\t{Is_reads}\t{score}\t{level}\n"
fo.write(res)
结果示例
由于结果行数太多，这里仅列举一些比较有代表性的情况供参考：
点击查看代码
Uniq complex high_MQ DP AD AF Is_reads Score Level
---- ------- ------- -- -- ---- -------- ------- -----
0 0 0 1 1 1.0 0 -96.0 E
0 0 1 1 1 1.0 0 -56.0 E
0 1 0 1 1 1.0 0 -56.0 E
0 1 1 1 1 1.0 0 -16.0 E
1 0 0 1 1 1.0 0 -56.0 E
1 0 1 1 1 1.0 0 -16.0 E
1 1 0 1 1 1.0 0 -16.0 E
1 1 1 1 1 1.0 0 24.0 D
0 0 0 2 1 0.5 0 -122.5 E
0 0 0 2 2 1.0 0 -83.0 E
1 1 1 6 4 0.67 2 63.43 C
1 1 1 6 5 0.83 0 56.45 C
1 1 1 6 5 0.83 1 73.45 C
1 1 1 6 5 0.83 2 80.49 B
1 1 1 6 6 1.0 0 77.0 C
1 1 1 6 6 1.0 1 94.0 B
1 1 1 6 6 1.0 2 101.04 B
1 1 1 6 6 1.0 3 106.44 B
1 1 1 8 7 0.88 1 89.72 B
1 1 1 8 7 0.88 2 96.76 B
1 1 1 8 7 0.88 3 102.16 B
1 1 1 8 8 1.0 0 92.0 B
1 1 1 8 8 1.0 1 109.0 B
1 1 1 8 8 1.0 2 116.04 B
1 1 1 8 8 1.0 3 121.44 A
1 1 1 8 8 1.0 4 126.0 A
拓展验证
对于这类评定机制没有标准，需要反复摸索条件、调试各证据项的评分规律以及后续分级阈值的程序，必然需要一些比较规范的验证机制。（与代码的 test 模块类似）
比如，有一些证据项情况，我们希望在调整完任何一部分的赋分或者分级阈值时，依然如预期的工作着。在代码 opts 中添加所需参数和预期结果，从而可以批量核查指定情况下，程序是否按预期得到了结果。
点击查看代码
def all_check_group():
""" 所有需要的数据核查 """
opts = [
{
'func' : call_score, # 需要运行的函数
'args' : { # 函数的参数
'AD': 10,
'AF': 1,
'Is_reads': 5,
'Uniq_region': 1,
'complex_region': 1,
'high_MQ_reads': 1,
},
'res' : 'A', # 预期结果
'call_back' : cal_level # func 函数运行后的结果调用的函数
},
{
'func' : call_score,
'args' : {
'AD': 6,
'AF': 1,
'Is_reads': 3,
'Uniq_region': 1,
'complex_region': 1,
'high_MQ_reads': 1,
},
'res' : 'B',
'call_back' : cal_level
},
]
pass_test_num = 0
for k in opts:
if 'call_back' in k and k['call_back']:
res = k['call_back'](k['func'](**k['args']))
if res != k['res']:
logging.warning(f"函数 \n\t{k['call_back'].__name__}({k['func'].__name__}(**{k['args']})) \
\n\t结果与预期不符:\n\t预期: {k['res']}\n\t实际: {res}\n")
else:
pass_test_num += 1
continue
res = k['func'](k['args'])
if res != k['res']:
logging.warning(f"函数 {k['func'].__name__}({k['args']}) \
结果与预期不符:\n\t预期: {k['res']}\n\t实际: {res}\n")
else:
pass_test_num += 1
logging.info(f'共 {len(opts)} 个测试, 成功 {pass_test_num} 个, 失败 {len(opts) - pass_test_num} 个.')
总结
该随笔及其部分代码仅仅是一个思路的描述，可拓展的地方很多，比如更全面的证据项，变异的分级、评级机制等都可商榷评估，以使其更为合理。另外该部分不包含这些证据项如何从 bam 中提取出来，后续的工作还有很多，继续加油干吧~~
