前面我们介绍了GIAB中的标准数据集，包括了7个标准品的数据，分别是HG001（即NA12878），一个德系犹太人家系（HG002/HG003/HG004）和一个中国汉族人家系（HG005/HG006/HG007），内容概要如下：

1. v4和v3版本标准数据集构建方法不同，v4添加了三代数据、基于组装和单体型的变异检测方法，变异位点数目更多；
2. 提供了数据库下载的网页地址和对应的ftp地址，以及通过wget -r下载整个目录的方法；
3. 由于ftp下面包含的数据很多，文中也明确了需要用那个文件进行准确性评估；

这篇文章介绍下，当我们使用自己开发的流程对某一标准品（例如NA12878）完成变异检测后，该如何评估结果的Precision（准确率）、Recall（召回率）和F1-score等指标。

最简单的思路

常用变异检测软件，如GATK、samtools和deepvariant的变异检测结果都是以vcf格式保。vcf文件中前五列分别是变异所在的CHROM染色体名称，POS突变位置，ID突变唯一标识，REF参考基因型，ALT突变基因型，后面FORMAT还记录了是纯和突变还是杂合突变信息。关于vcf格式说明，可以查看：https://samtools.github.io/hts-specs/VCFv4.2.pdf

vcf文件格式示例

这个思路存在的问题

为了说明这个问题，我们先给出一个例子，下图中ALT1和ALT2相对参考REF都存在一个插入，从图中上半部分可以看出属于相同的突变，并且插入序列相同。但是，当对齐的方式不同时，检测到的变异结果就会不同，如图中下半部分所示。

相同变异，由于比对对齐方式不同导致的变异检测结果不同

除此之外，两个相邻的SNV突变，有的软件在vcf中只有一行记录，有的软件会存放两行。所以，只通过位置和基因型来比较两个变异检测结果是存在问题的。
比较正确的方法时，将突变序列替换到突变位置，并且获取突变上下游一段序列，直接比较两段替换后的序列是否相同才是正确的，这种方法称之为“基于单体型的比较”。当然，如果上下游序列过长会降低计算效率，如果取得较短进行比较没有意义，所以下文介绍的软件默认取1k的序列进行比较。

基于单体型比较的软件

在Illuminia官方github上基于单体型比较的工具为：https://github.com/illumina/hap.py，用户可以根据安装说明进行安装。由于软件依赖较多，直接安装会比较麻烦，可以使用docker进行安装，前提是当前操作系统已经安装好docker。当软件下载解压完成后，直接进入解压后的目录，里面必须包含“Dockerfile”文件，然后直接运行docker build .即可，如下所示。

#第一步，下载软件包
wget https://github.com/Illumina/hap.py/archive/refs/tags/v0.3.15.zip
#第二步，解压软件包
unzip hap.py-0.3.15.zip 
#第三步，切换工作目录
cd hap.py-0.3.15
#第四步，构建镜像
docker build .

第四步运行完成后，会出现一串十六进制编码的字符串，这表示构建好镜像的ID（IMAGE ID ），也可以通过docker images命令查看本地镜像仓库中的镜像列表。镜像构建完成后，还可以通过docker查看软件的参数说明：

查看本地镜像仓库中的镜像

$ docker run -v `pwd`:/data 4c39d45b49ea /opt/hap.py/bin/hap.py -h         
2023-06-25 05:54:43,150 WARNING  No reference file found at default locations. You can set the environment variable 'HGREF' or 'HG19' to point to a suitable Fasta file.
usage: Haplotype Comparison [-h] [-v] [-r REF] [-o REPORTS_PREFIX]
                            [--scratch-prefix SCRATCH_PREFIX] [--keep-scratch]
                            [-t {xcmp,ga4gh}] [-f FP_BEDFILE]
                            [--stratification STRAT_TSV]
                            [--stratification-region STRAT_REGIONS]
                            [--stratification-fixchr] [-V] [-X]
                            [--no-write-counts] [--output-vtc]
                            [--preserve-info] [--roc ROC] [--no-roc]
                            [--roc-regions ROC_REGIONS]
                            [--roc-filter ROC_FILTER] [--roc-delta ROC_DELTA]
                            [--ci-alpha CI_ALPHA] [--no-json]
                            [--location LOCATIONS] [--pass-only]
                            [--filters-only FILTERS_ONLY] [-R REGIONS_BEDFILE]
                            [-T TARGETS_BEDFILE] [-L] [--no-leftshift]
                            [--decompose] [-D] [--bcftools-norm] [--fixchr]
                            [--no-fixchr] [--bcf] [--somatic]
                            [--set-gt {half,hemi,het,hom,first}]
                            [--filter-nonref] [--convert-gvcf-to-vcf]
                            [--gender {male,female,auto,none}]
                            [--convert-gvcf-truth] [--convert-gvcf-query]
                            [--preprocess-truth] [--usefiltered-truth]
                            [--preprocessing-window-size PREPROCESS_WINDOW]
                            [--adjust-conf-regions] [--no-adjust-conf-regions]
                            [--unhappy] [-w WINDOW]
                            [--xcmp-enumeration-threshold MAX_ENUM]
                            [--xcmp-expand-hapblocks HB_EXPAND]
                            [--threads THREADS]
                            [--engine {xcmp,vcfeval,scmp-somatic,scmp-distance}]
                            [--engine-vcfeval-path ENGINE_VCFEVAL]
                            [--engine-vcfeval-template ENGINE_VCFEVAL_TEMPLATE]
                            [--scmp-distance ENGINE_SCMP_DISTANCE]
                            [--lose-match-distance ENGINE_SCMP_DISTANCE]
                            [--logfile LOGFILE] [--verbose | --quiet]
                            [_vcfs [_vcfs ...]]
..............

使用demo数据进行评估

在下载软件中的example目录下，保存着可供测试的demo数据。在开始测试前，还需要下载参考基因组hg19.chr.fa，并构建其索引文件（samtools faidx hg19.chr.fa），然后才能开始测试，测试完成后会生成以-o参数指定的前缀文件：

#查看当前工作目录下的数据
hap.py-0.3.15 $ ls
CMakeLists.txt  Dockerfile                    example                 hg19.chr.fa      Jenkinsfile  RELEASES.md
configure.sh    Dockerfile.centos6            external                hg19.chr.fa.fai  LICENSE.txt  setup.cfg
doc             Dockerfile.ubuntu-with-tests  happy.requirements.txt  install.py       README.md    src
#测试命令行
hap.py-0.3.15 $ docker run -it -v `pwd`:/data 4c39d45b49ea /opt/hap.py/bin/hap.py /data/example/PG_performance.vcf.gz /data/example/performance.vcf.gz -o /data/test -r /data/hg19.chr.fa

demo运行结果日志

和基准变异集合进行比较

如果需要基于GIAB基准变异集合，对自己流程的分析结果进行准确性评估，还需要添加-f参数，该参数用于指定benchmark对应的bed文件，如果不指定METRIC.Precision值会偏低。
为什么不加会偏低呢？在回答这个问题前，先来说明一下benchmark的bed文件是什么。我们知道在基因组中，有些区域由于重复序列或者GC含量的关系，造成有些区域测序数据量比较低或比对错误比较高，那么这些区域是比较难进行SNV或者InDel检测的。为此，将一些复杂区域剔除，只保留“可检测性”较高的区域就生成了benchmark的bed文件和对应的vcf文件。
以样本HG001在GIAB数据库中NISTv4.2.1版本的“基准变异集合”为例，可以看到两个基因组版本都有基准变异集合的vcf文件和bed文件。

image.png

假设我们用自己的全基因组分析流程对样本HG001进行变异检测，最终得到的变异集合文件为HG001.GATK.vcf.gz，要与基准变异集合进行比较，来评估流程检测结果的Precision（精准率）、Recall（召回率）、F1-score，分析命令如下：

docker run -it -v `pwd`:/data 4c39d45b49ea   #挂载目录，指定镜像
     /opt/hap.py/bin/hap.py    #软件在容器中的路径
     /data/NISTv4.2.1/GRCh37/HG001_GRCh37_1_22_v4.2.1_benchmark.vcf.gz  #基准变异集合
     /data/HG001.GATK.vcf.gz  #自己流程分析到的变异集合
     -o /data/HG001  #输出目录和前缀
     -r /data/hg19.chr.fa   #参考基因组路径
     -f /data/NISTv4.2.1/GRCh37/HG001_GRCh37_1_22_v4.2.1_benchmark.bed  #基准变异集合对应的bed文件

未添加-f参数

添加-f参数后

可以看出添加-f参数Precision、Recall和F1 score都有明显提高，另外UNK标识落在bed区间外的变异数量。