敏感性、特异性（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 进行统计。
#!/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的最低深度？
统计不同深度下的假阴假阳性率，看在什么深度下其达到饱和。
