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

定义

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

计算方法

Sensitivity and specificity：完整定义

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
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 进行统计。

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

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
#!/usr/bin/env python
# Author: LI ZHIXIN
import sys
import pysam
from collections import OrderedDict
def 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.pos
def 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_dict
def 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_dict
def 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的最低深度？

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

敏感性、特异性、假阳性、假阴性（sensitivity and specificity）相关推荐

  1. 真阳性假阳性假阴性分割可视化

    1.分割掩码二值化 分割掩码转化为图像格式时会在分割边界处有很多灰度像素点,这将导致后续利用分割掩码和预测分割掩码进行处理时会在边界处出现很多噪声点,因此利用阈值将分割掩码转换为二值图,消除边界上的灰 ...

  2. 机器学习基础（一）混淆矩阵，真阳性(TP)，真阴性(TN)，假阳性(FP)，假阴性(FN)以及敏感性(Sensitivity)和特异性(Specificity)

    机器学习基础(一) 混淆矩阵 真阳性,真阴性,假阳性,假阴性 敏感性,特异性 混淆矩阵 混淆矩阵如下图:这里以是否有心脏病举例(二分类举例),列代表机器学习算法所做的预测,有心脏病还是没有心脏病,行代 ...

  3. 敏感性、特异性（sensitivity and specificity）| 假阳性、假阴性 | FDR | 第一类错误、第二类错误 | ROC | AUC...

    这些概念确实很难记忆,长时间不用很容易忘记.此文可以帮你快速回忆起这些概念,同时不涉及任何绕口的专业术语. 1. 记住基本框架,金标准和预测结果,没有这两个概念就没有敏感性和特异性了.以上指标都是用于 ...

  4. 准确率，召回率，mAP，ROC，AUC，特异性，敏感性，假阴性，假阳性

    P/R和ROC是两个不同的评价指标和计算方式,一般情况下,检索用准确率.召回率.F1.AP和mAP,分类.识别等用ROC和AUC(特异性和敏感性是ROC当中的一个部分). 准确率.召回率.F1.AP和 ...

  5. 敏感性、特异性、假阳性、假阴性

    敏感性.特异性.假阳性.假阴性是医学领域常用的评估指标. 敏感性:在金标准判断有病(阳性)的人群中,检测出阳性的几率.真阳性(检测出确实有病的能力) TPR = TP / ( TP+FN ) = TP ...

  6. 阳性，阴性，假阳性，假阴性，敏感度，特异性

    一般从医学角度说,阳性(positive),代表有病或者有病毒,阴性(negative),代表正常. 假阳性(false positive)是指因为种种原因把不具备阳性症状的人检测出阳性的结果.其实就 ...

  7. Sensitivity and specificity 敏感性与特异性

    Sensitivity and specificity 敏感性与特异性 敏感性 =真正例率=灵敏度(此说法不准确,翻译问题) 特异性=假正例率 参考 https://blog.csdn.net/roc ...

  8. True/False Positive True/False Negative (真/假阳性, 真/假阴性)

    真阳性 (True Positive) - 本身是对的,你当成对的了 假阳性 (False Positive) - 本身是错的,你当成对的了 真阴性 (True Negative) - 本身是错的,你 ...

  9. fNIRS中的假阳性和假阴性：问题、挑战和方法

    导读 本文强调了在进行功能性近红外光谱(fNIRS)研究时需要考虑和解决的一个重要问题,即无意中测量非神经血管耦合引起的fNIRS血流动力学反应的可能性.这些可能被误解为大脑活动,即"假阳性 ...

最新文章

  1. SQL Server 2008 安装过程中遇到“性能计数器注册表”..
  2. 金融数据分析之pdfplumber提取年报PDF关键数据（其他PDF数据通用）
  3. 2020-11-18 Ubuntu 安装 Chrome
  4. 关于在Eclipse里面启动了服务,但是localhost:8080无法访问的问题:
  5. QQ邮箱鸡肋存储型XSS漏洞利用
  6. acm 算法竞赛 时间
  7. 如何用Python批量打印PDF文档、Word文档、Excel表格、图片呢？
  8. ecu故障现象_发动机各传感器故障现象总结
  9. url采集工具_2022年1月6日更新：关键词URL采集工具最新版
  10. 搭建3款远程开发环境：Pycharm、Jupyter notebook以及code-server
  11. 世界卫生组织关于糖尿病、眼部疾病的相关数据整理
  12. 前端基础之《NodeJS（2）—模块化》
  13. ARKit之路－Depth API
  14. 我，程序员，马上35岁...
  15. Y05 - 017、猜小埋年龄游戏
  16. sh_update_eop 更新eop
  17. “天问一号”火星车命名由来
  18. HDU - 2153 仙人球的残影
  19. 超融合市场的战争远未结束，谁将最终胜出？
  20. 超简单的用PS(PhotoShop)转换png为ico，不用下什么插件什么玩意儿的

热门文章

  1. linux系统与window区别,Linux和windows操作系统有哪些区别
  2. python实现简单的图书管理系统
  3. 电脑视频转换成mp4格式，视频格式转换器转换
  4. 2022.04.11【读书笔记】|单细胞转录组概述
  5. appstore上架助手
  6. java为word、excel、pdf、ppt、图片添加图片水印（文字水印同理）
  7. 程序员是什么又代表这多少角色？你想过吗？
  8. java 半小时_java获取当前时间加半小时之后的时间
  9. pdf文件预览 浏览器窗口名修改
  10. 迭代法求一元三次方程