敏感性、特异性(sensitivity and specificity)| 假阳性、假阴性 | FDR | 第一类错误、第二类错误 | ROC | AUC
这些概念确实很难记忆,长时间不用很容易忘记。此文可以帮你快速回忆起这些概念,同时不涉及任何绕口的专业术语。
1. 记住基本框架,金标准和预测结果,没有这两个概念就没有敏感性和特异性了。以上指标都是用于衡量我们(预测)方法的效果的。考虑两个极端的方面,一、我的诊断/预测方法里的阳性包含了所有的阳性(金标准)个体,想实现这个很简单,就是任何一个人来看病我都说你得了艾滋病,这会带来什么结果呢?正面的就是所有的艾滋病个体我都可以给你判断出来,负面的呢?就是同时一大堆正常人给误判为艾滋病了;二、我的诊断/预测方法的阴性包含了所有的阴性(金标准)个体,一个极端的诊断方法就是任何一个人来看病我都说你没有艾滋病,这回漏掉了很多的阳性个体。当然大部分的诊断和预测方法都不会这么极端,核心意思就是我们必须同时考虑我们方法在阳性和阴性(犯第一类错误和第二类错误)两方面的综合表现。
2. 画表是你看到这几个指标的第一反应,横轴是金标准的阳性和阴性,纵轴是方法的阳性和阴性,我们预测的任何一个个体都会落到以下四个象限中的一个:
看图说话:
真阳性:左上格,金标准和预测都为真;
真阴性:右下格,金标准和预测都为假;
假阳性:右上格,就是我们的预测结果是假的阳性(预测出来的阳性其实是金标准的阴性);
假阴性:左下格,就是我们的预测结果是假的阴性(预测出来的阴性其实是金标准的阳性);
敏感性sensitivity:顾名思义,就是我们的方法对阳性的敏感度,换句话说就是我给你一堆金标准的阳性结果,你能判断出多少来。就像缉毒犬一样,给你一堆含毒品的行李箱,你能闻出多少出来。这个值命名为敏感性还是很形象的。
特异性specificity:顾名思义,就是你的方法能指定地判断出阴性的能力,我给你一堆金标准的阴性结果,你会不会误判出一些阳性的结果。这个命名还是欠缺形象性,叫准阴性更贴切。
FDR:常用于多重检验下的p-value矫正;伪发现率,其意义为是 错误拒绝(拒绝真的(原)假设)的个数占所有被拒绝的原假设个数的比例的期望值。就是我的所有预测阳性中有多少是误判的假阳性。以缉毒犬为例,我的狗闻出了一堆箱子,里面有多少是不含毒品的。在差异基因检测中,使用FDR就是用于控制假阳性的比率(通常为0.05),是不是又把所有知识串起来了。
第一类错误、第二类错误:显然我们的方法里有且仅有两类错误,通常我们先考虑我们预测的阳性结果,里面有多少错误(其实是金标准的阴性),这些错误就是第一类错误,也就是上面的假阳性。在考虑我们预测的阴性结果,里面有多少是假阴性,也就是有多少第二类错误。
False positive rate (FPR):假阳性率,给你一堆金标准的阴性,有多少是假的阴性,等于1-specificity,1-准阴性不就是假阴性吗。
Confusion matrix混淆矩阵 / Contingency table列联表 (这个表的两种名称)
基于sensitivity和specificity衍生出来的两个概念:只能用于二分类模型的评价。
怎么全面地评价一个二分类模型的好坏,模型的其他指标都依赖一个threshold,单一的threshold是有偏的。
ROC:receiver operating characteristic,每个分类器作出的预测呢,都是基于一个probability score的。一般默认的threshold呢都是0.5,如果probability>0.5,那么这个sample被模型分成正例了哈,反之则是反例。ROC曲线是一系列threshold(比如逻辑斯蒂分类问题,取多少阈值来划分类别)下的(FPR,TPR)数值点的连线。不同阈值下模型的sensitivity和specificity不同。
AUC:areas-under-the-curve,曲线下的面积(AUC)越大,或者说曲线更接近左上角(true positive rate=1, false positive rate=0)。AUC的意义:分别随机从正负样本集中抽取一个正样本,一个负样本,正样本的预测值大于负样本的概率。Wilcoxon-Mann-Witney Test。
- (0,0):TP=0,FP=0,可以发现该分类器预测所有的样本都为负样本(Negative)
- (1,1):TN=0,FN=0,可以发现该分类器预测所有的样本都为正样本(Positive)
- (0,1):FP=0,FN=0,它将所有的样本都正确分类
- (1,0):TP=0,TN=0,它将所有的样本都错误分类
继续深入:
精确度Precision:很容易与敏感性搞混,分子是一样的,但分母不同,精确度的分母是我们方法预测出来的阳性,敏感性是金标准的阳性。精确度表示我们预测了一堆阳性结果,里面有多少是靠谱的。
准确度Accuracy:和Precision中文意思是一样的,Precision指的是测量结果的一致性。Accuracy指的是测量结果和真实结果的接近程度。
Recall:TPR,真阳性率,和sensitivity一个东西。
F1 score:F1-score(均衡平均数)是综合考虑了模型查准率和查全率的计算结果,取值更偏向于取值较小的那个指标。F1-score越大自然说明模型质量更高。但是还要考虑模型的泛化能力,F1-score过高但不能造成过拟合,影响模型的泛化能力。
以下是我曾经的一个项目,不感兴趣的可以不看。
终于要用到这个玩意了,很激动,主要统计假阴假阳性率。
我的任务:
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的最低深度?
统计不同深度下的假阴假阳性率,看在什么深度下其达到饱和。
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 一个奇形怪状的面试题:Bean中的CHM要不要加volatile?
· [.NET]调用本地 Deepseek 模型
· 一个费力不讨好的项目,让我损失了近一半的绩效!
· .NET Core 托管堆内存泄露/CPU异常的常见思路
· PostgreSQL 和 SQL Server 在统计信息维护中的关键差异
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 如何使用 Uni-app 实现视频聊天(源码,支持安卓、iOS)
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)