敏感性、特異性（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。

參考：如何理解機器學習和統計中的AUC？

• (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的最低深度？

統計不同深度下的假陰假陽性率，看在什么深度下其達到飽和。