NA12878 全外显子数据分析结果比较

一直都比较有兴趣分析 human DNA 变异，参照 gatk 流程分析了 NA12878 的全外数据，正好有参考变异集，比较下分析准确性，也好参数调整优化之类，吧吧吧，下面就是一些详细步骤粘贴整理

NA12878.vcf 处理

# wget -c ftp://ftp-trace.ncbi.nih.gov/giab/ftp/data/NA12878/Garvan_NA12878_HG001_HiSeq_Exome/project.NIST.hc.snps.indels.vcf
hc_snp_indel=project.NIST.hc.snps.indels.vcf
# wget -c ftp://ftp-trace.ncbi.nih.gov/giab/ftp/data/NA12878/Garvan_NA12878_HG001_HiSeq_Exome/nexterarapidcapture_expandedexome_targetedregions.bed
intervals=nexterarapidcapture_expandedexome_targetedregions.bed

java -jar GenomeAnalysisTK.jar \
    -T SelectVariants \
    -R $genome \
    -V $hc_snp_indel \
    -L $intervals \
    -sn NIST7035 \
    -select 'DP > 5' \
    --excludeNonVariants \
    -o ref.NIST7035.hc.snps.indels.vcf.gz
# 选择特定区间内 DP 大于5，样品 NIST7035 的变异位点

vcf 处理

ref: https://software.broadinstitute.org/gatk/documentation/article.php?id=2806

# SNP
java -jar GenomeAnalysisTK.jar \
    -T SelectVariants \
    -R $genome \
    -L $intervals \
    -selectType SNP \
    -V  NIST7035.vcf.gz \
    -o NIST7035.snp.vcf.gz
java -jar GenomeAnalysisTK.jar \
    -T VariantFiltration \
    -R $genome \
    --filterExpression "DP < 5" --filterName "LowDepth"  \
    --filterExpression "FS > 60.0 || SOR > 3.0" --filterName "StrandBias"  \
    --filterExpression "QD < 2.0  || MQ < 40.0 || MQRankSum < -12.5" \
    --filterName "snp_Filter"  \
    --filterExpression "ReadPosRankSum < -8.0" --filterName "endOfRead"  \
    -V  NIST7035.snp.vcf.gz \
    -o  NIST7035.snp.filter.vcf.gz

# Indel
java -jar GenomeAnalysisTK.jar \
    -T SelectVariants \
    -R $genome \
    -L $intervals \
    -selectType INDEL \
    -V  NIST7035.vcf.gz \
    -o NIST7035.indel.vcf.gz
java -jar GenomeAnalysisTK.jar \
    -T VariantFiltration \
    -R $genome \
    --filterExpression "DP < 5" --filterName "LowDepth"  \
    --filterExpression "QD < 2.0 || FS > 200.0" \
    --filterName "InDels_Filter" \
    --filterExpression "ReadPosRankSum < -20.0" --filterName "endOfRead" \
    -V NIST7035.indel.vcf.gz \
    -o NIST7035.indel.filter.vcf.gz

# MergeVcf
java -jar picard.jar MergeVcfs \
    I=NIST7035.indel.filter.vcf.gz \
    I=NIST7035.snp.filter.vcf.gz \
    O=NIST7035.filter.vcf.gz\
    D=$genome.dict

compare use bcftools

bcftools isec NIST7035.filter.vcf.gz ref.NIST7035.hc.snps.indels.vcf.gz -p bcftools
# ref.vcf.gz, alt.vcf.gz 使用 bgzip 压缩，同时 tabix 进行 index

# output: README.txt
# This file was produced by vcfisec.
# The command line was:   bcftools isec  -p bcftools NIST7035.filter.vcf.gz ref.NIST7035.hc.snps.indels.vcf.gz
#
# Using the following file names:
# bcftools/0000.vcf  for records private to  NIST7035.filter.vcf.gz
# bcftools/0001.vcf  for records private to  ref.NIST7035.hc.snps.indels.vcf.gz
# bcftools/0002.vcf  for records from NIST7035.filter.vcf.gz shared by both  NIST7035.filter.vcf.gz ref.NIST7035.hc.snps.indels.vcf.gz
# bcftools/0003.vcf  for records from ref.NIST7035.hc.snps.indels.vcf.gz shared by both      NIST7035.filter.vcf.gz ref.NIST7035.hc.snps.indels.vcf.gz

cat ref.NIST7035.hc.snps.indels.vcf.gz | grep -v '^#' | wc -l    # 54077

cat 0000.vcf | grep -cv '^#'    # 3863, 假阳性
cat 0001.vcf | grep -cv '^#'    # 2131, 假阴性
cat 0002.vcf | grep -cv '^#'    # 51946
cat 0003.vcf | grep -cv '^#'    # 51946

cat 0000.vcf| grep -v '^#' | grep PASS | wc -l    # 1334
cat 0002.vcf| grep -v '^#' | grep PASS | wc -l    # 47143

不考虑snp、indels过滤情况下，计算假阳性、假阴性率。

假阳性: 3863/54077 = 7.14 %

假阴性: 2131/54077 = 3.94 %

哈哈哈哈，瞎了，没法看的样子。还是考虑下变异过滤结果吧

假阳性: 1334/54077 = 2.46 %

假阴性: (2131+51946-47143)/54077 = 12.82 %

更加让人没话说，这个假阴性太感人，过滤标准不太合适啊

# 先看被过滤掉的这部分阴性
cat 0002.vcf | grep -v '^#' | grep -v PASS | awk '{n[$7]+=1}END{for (i in n){print i,n[i]}}'
# LowDepth 11
# StrandBias;snp_Filter 277
# LowDepth;snp_Filter 2
# InDels_Filter 265
# StrandBias 3098
# LowDepth;StrandBias 2
# snp_Filter 1148

过滤掉位点最严重的两项分别是StrandBias, snp_Filter，对应过滤指标分别为FS > 60.0 || SOR > 3.0, QD < 2.0 || MQ < 40.0 || MQRankSum < -12.5，转换vcf文件到tab格式（脚本：github），对其中一些指标进行分布检查。

先看 StrandBias，对 FS, SOR 两个参数分析，分布结果如下，可看出被过滤的位点基本上都是由于 SOR 参数（作图时FS>60无显示），与 NA12878 变异集相比，需调整 SOR 临界值到更大水平。

再看 snp_Filter 对应的QD < 2.0 || MQ < 40.0 || MQRankSum < -12.5，其中 MQRankSum 对应67个位点，QD 对应557个位点，MQ 对应 1168 个位点，看出 MQ, QD 标准需要进行调整。

以现在质控指标对文件0001.vcf 2131 变异质控，QD 过滤掉949，MQRankSum 过滤掉14，MQ 过滤掉195，FS 过滤掉135，根据这部分参考变异集质量指标，同样 QD, MQ, FS 参数需要进行调整。

另，vcftools 也可以很方便进行两个 vcf 文件的比较，还没 bcftools 对 vcf 文件必须排序、index的要求，输出结果也相应不同，因为我需要看阳性、阴性位点是否被过滤以及过滤情况，就不再用 vcftools 整理了，具体命令可参考如下：

# http://vcftools.github.io/
vcftools --vcf ref.vcf --diff alt.vcf --out prefix
# ref.vcf.gz    # 则需要使用--gzvcf
# alt.vcf.gz    # 则需要使用--gzdiff

质控参数调整思路

根据质控结果比较来看，参控参数都还需要进行再调整，一个初步的参数调整思路是，对高质量变异集各项质控指标分布进行统计，然后进行阈值选择；或者对 NA12878 全基因组水平的变异进行 VQSR 分析，统计高质量变异位点质控指标分布