医学、机器学习等等，在统计结果时时长会用到这两个指标来说明数据的特性。

定义

敏感性：在金标准判断有病（阳性）人群中，检测出阳性的几率。真阳性。（检测出确实有病的能力）

特异性：在金标准判断无病（阴性）人群中，检测出阴性的几率。真阴性。（检测出确实没病的能力）

：得到了阳性结果，但这个阳性结果是假的。即在金标准判断无病（阴性）人群中，检测出为阳性的几率。（没病，但却检测结果说有病），为误诊率。

假阴性率：得到了阴性结果，但这个阴性结果是假的。即在金标准判断有病（阳性）人群中，检测出为阴性的几率。（有病，但却检测结果说没病），为漏诊率。

计算方法

Sensitivity and specificity：

 

True Positive （真正, TP）被模型预测为正的正样本；可以称作判断为真的正确率True Negative（真负 , TN）被模型预测为负的负样本 ；可以称作判断为假的正确率False Positive （假正, FP）被模型预测为正的负样本；可以称作误报率False Negative（假负 , FN）被模型预测为负的正样本；可以称作漏报率True Positive Rate（真正率 , TPR）或灵敏度（sensitivity） TPR = TP /（TP + FN） 正样本预测结果数 / 正样本实际数 True Negative Rate（真负率 , TNR）或特指度（specificity） TNR = TN /（TN + FP） 负样本预测结果数 / 负样本实际数 False Positive Rate （假正率, FPR） FPR = FP /（FP + TN） 被预测为正的负样本结果数 /负样本实际数False Negative Rate（假负率 , FNR） FNR = FN /（TP + FN） 被预测为负的正样本结果数 / 正样本实际数

　　

假阳性率=假阳性人数÷金标准阴性人数

即： 假阳性率=b÷（b+d）

		金标准	金标准
		阳性（+）	阴性（-）	合计
某筛检方法	阳性（+）	a	b	a+b
某筛检方法	阴性（-）	c	d	c+d
合计		a+c	b+d	N

公式为:假阳性率=b/(b+d)×100%

(b:筛选为阳性,而标准分类为阴性的例数;d:阴性一致例数)

假阴性率=假阴性人数÷金标准阳性人数

即： β=c÷（a+c）

---

终于要用到这个玩意了，很激动，主要统计假阴假阳性率。

我的任务：

1. 评估Pacbio MHC variation calling 结果（CCS/non-CCS）与Hiseq数据结果的一致性。

2. 分别在不同深度梯度的区域完成以上评估，推断PB MHC做variation calling的最低深度。

这里要将一个位点分为SNP、REF 和 LowQual，然后只去 SNP 和 REF 进行统计。

这是我一下午写出来的统计代码：

#!/usr/bin/env python# Author: LI ZHIXINimport sysimport pysamfrom collections import OrderedDictdef classify_DP(depth):    if depth > 101:        return 21    return ((depth-1)//5+1)def parse_rec(rec):    sample = list(rec.samples)[0]    # filter the Invalid line    if not ('GQ' or 'GT' or 'DP') in rec.samples[sample].keys() or len(rec.alleles) <= 1:        # continue        return 1, "LowQual", rec.pos    # filter the LowQual    if rec.samples[sample]['GQ'] < 30:        return rec.samples[sample]['DP'], "LowQual", rec.pos    # filter the indel    flag = 0    for one in rec.alleles:        if len(one) != len(rec.ref):            flag = 1    if flag == 1:        return rec.samples[sample]['DP'], "LowQual", rec.pos    if rec.samples[sample]['GT'] != (0, 0): # rec.qual > 30        # variation_dict[rec.pos] = ["snp", rec.alleles]        return rec.samples[sample]['DP'], "snp", rec.pos      elif rec.samples[sample]['GT'] == (0, 0):        # variation_dict[rec.pos] = ["ref", rec.alleles]        return rec.samples[sample]['DP'], "ref", rec.posdef read_gvcf(gvcf_file_path):    variation_dict = OrderedDict()    for i in range(1,22):        variation_dict[i] = {}        for j in ('LowQual', 'snp', 'ref'):            variation_dict[i][j] = []    # pos_list = []    gvcf_file = pysam.VariantFile(gvcf_file_path)    for rec in gvcf_file.fetch('chr6',28477796,33448354):        DP, pos_type, pos = parse_rec(rec)        if DP < 1 or DP > 20:            continue        # DP = classify_DP(DP)        variation_dict[DP][pos_type].append(pos)        # print(pos, DP, pos_type)    gvcf_file.close()    # return variation_dict, pos_list    return variation_dictdef read_hiseq_gvcf(gvcf_file_path):    variation_dict = OrderedDict()    # for i in range(1,22):    # variation_dict[i] = {}    for j in ('LowQual', 'snp', 'ref'):        variation_dict[j] = []    # pos_list = []    gvcf_file = pysam.VariantFile(gvcf_file_path)    for rec in gvcf_file.fetch('chr6',28477796,33448354):        DP, pos_type, pos = parse_rec(rec)        DP = classify_DP(DP)        variation_dict[pos_type].append(pos)        # print(pos, DP, pos_type)    gvcf_file.close()    # return variation_dict, pos_list    return variation_dictdef show_dict_diff_DP(Hiseq_unified_variation_dict, PB_non_CCS_variation_dict, outf, outf2):    for DP in range(1,21):        Hiseq_snp = set(Hiseq_unified_variation_dict['snp'])        Hiseq_ref = set(Hiseq_unified_variation_dict['ref'])        Hiseq_lowqual = set(Hiseq_unified_variation_dict['LowQual'])        PB_snp = PB_non_CCS_variation_dict[DP]['snp']        PB_ref = PB_non_CCS_variation_dict[DP]['ref']        PB_lowqual = PB_non_CCS_variation_dict[DP]['LowQual']        total = set(PB_snp + PB_ref + PB_lowqual)        Hiseq_snp = total & Hiseq_snp        Hiseq_ref = total & Hiseq_ref        Hiseq_lowqual = total & Hiseq_lowqual        PB_snp = set(PB_snp)        PB_ref = set(PB_ref)        PB_lowqual = set(PB_lowqual)        a = len(Hiseq_snp & PB_snp)        b = len(Hiseq_ref & PB_snp)        c = len(Hiseq_lowqual & PB_snp)        d = len(Hiseq_snp & PB_ref)        e = len(Hiseq_ref & PB_ref)        f = len(Hiseq_lowqual & PB_ref)        g = len(Hiseq_snp & PB_lowqual)        h = len(Hiseq_ref & PB_lowqual)        i = len(Hiseq_lowqual & PB_lowqual)        Low_total = (g+h+i)/(a+b+c+d+e+f+g+h+i)        if (a+b) == 0:            PPV = "NA"        else:            PPV = a/(a+b)            PPV = "%.4f"%(PPV)        if (a+d) == 0:            TPR = "NA"        else:            TPR = a/(a+d)            TPR = "%.4f"%(TPR)        print(str(DP)+" :\n", a,b,c,"\n",d,e,f,"\n",g,h,i,"\n", file=outf2, sep='\t', end='\n')        print(DP, TPR, PPV, "%.4f"%Low_total, file=outf, sep='\t', end='\n')with open("./depth_stat.txt", "w") as outf:    print("Depth", "TPR", "PPV", "Low_total", file=outf, sep='\t', end='\n')    outf2 = open("raw.txt", "w")    Hiseq_unified_variation_dict = read_hiseq_gvcf("./hiseq_call_gvcf/MHC_Hiseq.unified.gvcf.gz")    PB_non_CCS_variation_dict = read_gvcf("./non_CCS_PB_call_gvcf/MHC_non_CCS.unified.gvcf.gz")    show_dict_diff_DP(Hiseq_unified_variation_dict, PB_non_CCS_variation_dict, outf, outf2)    outf2

　　

又碰到一个高级python语法：在双层循环中如何退出外层循环？ 我用了一个手动的flag，有其他好方法吗？

如何统计下机数据的覆盖度和深度？当然要比对之后才能统计，而且还要对比对做一些处理。

在计算一个位点是否是SNP、indel、Ref时，不仅要考虑ref、alts、qual、GQ，而且必须要把GT、DP考虑在内，所以说还是比较复杂的。

 

最后如何分析第二个问题，call variation的最低深度？

统计不同深度下的假阴假阳性率，看在什么深度下其达到饱和。