转载

GATK 4.0 WGS germline call variant

GATK升级4.0版了，作为人类call variant的金标准软件，加上其强大的团队，每次重大更新都会给使用者带来一点新的东西（或者说是改变），我也正好整理下，将GATK基本分析流程过渡到4.0版本

GATK4.0最明显的变化是其命令调用发生了改变，可以看看这个就明白了 https://software.broadinstitute.org/gatk/documentation/quickstart

GATK4.0还集成了picard工具以及增加了SortSam功能，这样Germline short variant discovery流程只要GATK一个软件即可完成sam/bam文件到vcf文件全流程

整体上Germline call variant分析思路没有变，所以我按照 Germline short variant discovery (SNPs + Indels) 教程上所示流程图，将GATK4.0基本用法再过一遍，主要是想再对GATK的VQSR质控理解下（之前一直没好好看其说明）

GATK call variant（snp/indel）

测试数据KPGP，下载地址： ftp://ftp.kobic.re.kr/pub/KPGP/2017_release_candidate_update_only/WGS/KPGP-00216/

参考基因组：GATK的hg38文件：

nohup wget wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.fai &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.dict &

使用fastp对测序数据质控，非常好用的一个质控软件，海普洛斯开发的

~/biosoft/fastp/fastp -i KPGP-00216_L1_R1.fq -I KPGP-00216_L1_R2.fq -o KPGP-00216_L1_R1.clean.fq -O KPGP-00216_L1_R2.clean.fq

可以看看质控报告fastp.html，蛮好的用的

压缩下原始数据，节约下空间，暂时不删（删了也可以）

tar cf - KPGP-00216_L1_R1.fq |pigz >KPGP-00216_L1_R1.fq.tar.gz

bwa比对

  1. 先建bwa的索引

    bwa index -a bwtsw -p gatk_hg38 ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta &

  2. 比对

    nohup bwa mem -t 4 -R '@RG/tID:KPGP00216/tPL:illumina/tLB:WGS/tSM:KPGP00216' ~/reference/index/bwa/gatk_hg38/gatk_hg38 KPGP-00216_L1_R1.clean.fq KPGP-00216_L1_R2.clean.fq 1>KPGP-00216_L1.sam 2>/dev/null &

将sam格式转化为bam格式，并且排序

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" SortSam -I KPGP-00216_L1.sam -O KPGP-00216_L1.sorted.bam -SO coordinate --CREATE_INDEX true

标记PCR重复序列

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" MarkDuplicates -I KPGP-00216_L1.sorted.bam -O KPGP-00216_L1.sorted.marked.bam -M KPGP-00216_L1.metrics

这个自己理解下哈，GATK里并没有建议一定要做

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" FixMateInformation -I KPGP-00216_L1.sorted.marked.bam -O KPGP-00216_L1.sorted.marked.fixed.bam -SO coordinate --CREATE_INDEX true

对参考基因组fa生成dict文件（其实在下载基因组时就该将这个dict文件一起下载下来）

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk CreateSequenceDictionary -R Homo_sapiens_assembly38.fasta -O Homo_sapiens_assembly38.dict

BQSR，也就是Recalibration Base Quality Score

  1. 先用FTP链接下GATK的resource database

    主机：ftp.broadinstitute.org  
用户：gsapubftp-anonymous

  2. 下载BQSR需要的注释位点VCF文件

    wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz.tbi
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi

  3. 运行BQSR，先BaseRecalibrator生成一个recal.data.table中间文件，再运行ApplyBQSR

    ~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" BaseRecalibrator -R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta -I KPGP-00216_L1.sorted.marked.fixed.bam -O recal_data.table --known-sites ~/annotation/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz --known-sites ~/annotation/GATK/resources/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz --known-sites ~/annotation/GATK/resources/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz

接下来则是使用HaplotypeCaller对上述的bam文件进行call variant，就是寻找变异的过程，生成一个VCF文件用于后续位点的质控和注释，最简单的单样本call variant如命令所示，从一个bam到一个vcf文件

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx12G HaplotypeCaller /
    -R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta /
    -I KPGP-00216_L1.sorted.marked.fixed.BQSR.bam /
    --dbsnp ~/annotation/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz /
    -O KPGP-00216_L1.raw.vcf

按照我的电脑配置，用了21h(组装的台式机，非服务器~)，为了加快速度，可以采用下述分染色体来call variant(也就是一个染色体生成一个VCF文件，22条常染色体 + XY性染色体，总共24个VCF文件)，当然不是一个个运行，使用shell命令行for循环加上 & 提交后台运行，如下所示： 

for i in {1..22} X Y;do (~/biosoft/GATK4.0/gatk-4.0.5.1/gatk HaplotypeCaller -R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta -I KPGP-00216_L1.sorted.marked.fixed.BQSR.bam --dbsnp ~/annotation/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz -L chr$i -O KPGP-00216_L1.chr$i.raw.vcf &);done

然后再用GatherVcfs将多个染色体的vcf文件合并

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk GatherVcfs $(for i in {1..22} X Y;do echo "-I KPGP-00216_L1.chr${i}.raw.vcf";done) -O KPGP-00216_L1_gather.raw.vcf

如果你有多个样本（或者大样本的一个队列），那么首选用g.vcf，其是GTAK HaplotypeCaller生成第一个中间文件，可以用于后续的GenotypeGVCFs进行joint genotyping；g.vcf简单的说是储存了所有位点的变异信息（不做任何可信度等的过滤）的vcf文件；其使用方式就上述HaplotypeCaller工具上加上 -ERC GVCF ，然后输出是XXX.g.vcf

有一点需要注意的是：

按照以前的思路，将多个样本的g.vcf文件一起通过的GenotypeGVCFs进行joint genotyping；但是！GATK 4.0的GenotypeGVCFs只支持a single single-sample GVCF，a single multi-sample GVCF created by CombineGVCFs 以及a GenomicsDB workspace created by GenomicsDBImport；所以之前的方法已经失效了，你在用GenotypeGVCFs前需要将多个样本的g.vcf文件用CombineGVCFs方式或者GenomicsDBImport方式合并成一个文件，前者（比较传统）是一个总的g.vcf文件，后者是一个GenomicsDB（XX.db）

这里对这个先不做详细说明啦，后续笔记会做尝试

另外对于多样本是分开cal还是joint call，可以参考下： 生信笔记：call snp是应该一起call还是分开call？

GTAK的VQSR质控

接下来则的对获得的VCF文件做质控，对于人物种的WGS测序，首选VQSR（Variant Quality Score Recalibration），其原理可以看 Variant Quality Score Recalibration (VQSR)-High-level overview ，个人理解简单的说：

  1. GATK认为VQSR比根据各种annotations进行hard-filtering过滤要好，减少了人为阈值的局限性，避免了一刀切的弊端，从而在sensitivity和specificity之间达到一定的平衡
  2. VQSR根据机器学习算法从highly validated变异位点数据集（每个位点的annotation profile，一般使用5-8个annotation）中获取到good variants/bad variants
  3. 根据上述的位点从我们自己数据集中挑选出一个变异子集（probably true positives）来建模训练，获得一个可识别good variants的模型；bad variants的模型也是如此获得
  4. 然后根据上述获得的模型，对自己数据集的每个变异位点进行一个总的打分
  5. 最后根据设定的sensitivity阈值对变异位点进行过滤

更简单的说：

利用自身数据和已知变异位点集的overlap，通过GMM模型构建一个分类器来对变异数据进行打分，从而评估每个位点的可行度

如果官网文档原理看的不是很清楚，可以参考这个（个人觉得讲的很好） GATK4.0和全基因组数据分析实践（下） 中的VQSR原理部分

所以可看出，VQSR对你的数据集的变异位点的数目有一定的要求；如果你的变异位点数目过少（如样本数较少的WES或者小panel测序），那么是无法做VQSR的；GATK对于样本不达标的WES，建议用1000 Genomes Project中的数据代替（将1000G的bam文件和WES样本的bam生成g.vcf文件，再一起做joint call）

GTAK的VQSR主要分两步骤进行：

  • VariantRecalibrator 对应我上述理解原理的1-4步骤，对每个变异位点打分注释VQSLOD value，从而生成一个recalibration文件以及一个xxx.plots.R.pdf(Gaussian mixture model plots)
  • ApplyRecalibration 对应5步骤，根据recalibration文件生成recalibrated VCF文件，并且根据过滤参数进行过滤（标记PASS）

所以当你tranche sensitivity（–truth-sensitivity-filter-level）设置为99.9%时，则表示将VQSLOD值高于整体99.9%的变异位点标记为PASS，表示通过过滤阈值，剩下的位点则被认为是假阳性的了；这个参数可以根据你的具体需求来设定，看你是需要more specific or more sensitive（tranche从90到100，specific随之降低，而sensitive随之升高）；或者生成多个tranche的结果，从中挑选你满意的阈值（可以看结果文件SNP中的xxx.tranches，而indel是没这个文件的；人全基因组Ti/Tv（转换和颠换的比例）的值一般在2.0-2.1…）

GATK的VQSR用的annotations指标有以下几种（VQSR的文档中还提到了InbreedingCoeff，这个是在大于10个样本中才使用的一个指标；如果样本过少或者是closely related samples(such as a family)的话，建议剔除；而且QD不适合用于外显子测序的数据的VQSR）：

  • Coverage（DP）
  • QualByDepth（QD）
  • FisherStrand (FS)
  • StrandOddsRatio (SOR)
  • RMSMappingQuality (MQ)
  • MappingQualityRankSumTest (MQRankSum)
  • ReadPosRankSumTest (ReadPosRankSum)

按照官方教程，SNP的VQSR过滤，选用的resource datasets为：

  • HapMap，hapmap_3.3.hg38.vcf.gz， truth=true 表示VQSR将这个数据集中的变异位点作为真位点true sites， training=true 表示VQSR将true sites用于训练recalibration model，并赋予这些变异位点prior likelihood值为 Q15 (96.84%)
  • Omni，1000G_omni2.5.hg38.vcf.gz， truth=true ， training=false （文档中写着是true，参数建议中写着的是false。。。我就按照参数上的来了）， Q12 (93.69%)
  • 1000G，1000G_phase1.snps.high_confidence.hg38.vcf.gz， truth=false 表示VQSR考虑到在1000G数据集中的不仅包含了true variants还有false positives， training=true ， Q10 (90%)
  • dbSNP，dbsnp_146.hg38.vcf.gz， truth=false 表示VQSR未将dbSNP数据集中的位点作为可信数据集， training=false 表示不用于训练数据集， known=true 表示stratify output metrics such as Ti/Tv ratio by whether variants are present in dbsnp or not， Q2 (36.90%)

INDEL的VQSR过滤，选用的resource datasets为：

  • Mills，Mills_and_1000G_gold_standard.indels.hg38.vcf.gz， truth=true ， training=true ， Q12 (93.69%)
  • dbSNP，dbsnp_146.hg38.vcf.gz， truth=false ， training=false ， known=true ， Q2 (36.90%)

其中几个resource datasets已经在BQSR中下载了，剩下的也是按照其方法下载即可

对于GATK的VQSR，snp和indel过滤是分开的，两者参数整体上差不多，但有点略微的区别

SNP的VQSR的过滤，一般的参数可以如下：

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk VariantRecalibrator /
-R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta /
-V KPGP-00216_L1_gather.raw.vcf /
--resource hapmap,known=false,training=true,truth=true,prior=15.0:/home/anlan/annotation/GATK/resources/bundle/hg38/hapmap_3.3.hg38.vcf.gz /
--resource omni,known=false,training=true,truth=false,prior=12.0:/home/anlan/annotation/GATK/resources/bundle/hg38/1000G_omni2.5.hg38.vcf.gz /
--resource 1000G,known=false,training=true,truth=false,prior=10.0:/home/anlan/annotation/GATK/resources/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz /
--resource dbsnp,known=true,training=false,truth=false,prior=2.0:/home/anlan/annotation/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz /
-an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an DP /
-mode SNP /
-O KPGP-00216_L1.snp.recal /
--tranches-file KPGP-00216_L1.snp.tranches /
--rscript-file KPGP-00216_L1.snp.plots.R

这里还有个参数要提一下， -tranche 默认是输出[100,99.9,99.0,90.0]4个tranche阈值的统计结果，如果想看其他阈值的结果，需要自行加上；结果就是看KPGP-00216_L1.snp.tranches（还有图形展示的KPGP-00216_L1.snp.tranches.pdf），而KPGP-00216_L1.snp.recal文件则是用于ApplyVQSR的

Applying recalibration/filtering to SNPs， --truth-sensitivity-filter-level 可以根据自身需求来设定 

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk ApplyVQSR /
-R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta /
-V KPGP-00216_L1_gather.raw.vcf /
-O KPGP-00216_L1.snp.VQSR.vcf /
--truth-sensitivity-filter-level 99.5 /
--tranches-file KPGP-00216_L1.snp.tranches /
--recal-file KPGP-00216_L1.snp.recal /
-mode SNP

INDEL的VQSR过滤，这里 -an 参数跟SNP有点不同，没有MQ，并且加了 --max-gaussians 4 用于设定Gaussians（clusters of variants that have similar properties）的数目，即减少聚类的组数，从而使得每个组的变异位点数目达到要求 

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk VariantRecalibrator /
-R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta /
-V KPGP-00216_L1_gather.raw.vcf /
--max-gaussians 4 /
--resource mills,known=false,training=true,truth=true,prior=12.0:/home/anlan/annotation/GATK/resources/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz /
--resource dbsnp,known=true,training=false,truth=false,prior=2.0:/home/anlan/annotation/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz /
-an QD -an DP -an FS -an SOR -an ReadPosRankSum -an MQRankSum /
-mode INDEL /
-O KPGP-00216_L1.indel.recal /
--tranches-file KPGP-00216_L1.indel.tranches /
--rscript-file KPGP-00216_L1.indel.plots.R

Applying recalibration/filtering to INDEL

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk ApplyVQSR /
-R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta /
-V KPGP-00216_L1_gather.raw.vcf /
-O KPGP-00216_L1.indel.VQSR.vcf /
--truth-sensitivity-filter-level 99.0 /
--tranches-file KPGP-00216_L1.indel.tranches /
--recal-file KPGP-00216_L1.indel.recal /
-mode INDEL

如果VQSR不能使用的话，那么只能用硬过滤（hard-filtering）了，这个注释指标每个人都可能用的不一样，可以先从GATK官方建议的指标入手，然后再看看文献中应用的指标，最后再结合自身需要来设定

以上是个人对GATK的一些理解，有些理解的不够好的地方有请联系我纠正下（尤其的VQSR那块，有些地方还云里雾里的）

参考资料：

Variant Quality Score Recalibration (VQSR)

(howto) Recalibrate variant quality scores = run VQSR

GATK流程

(How to) Consolidate GVCFs for joint calling with GenotypeGVCFs

Germline short variant discovery (SNPs + Indels)

本文出自于 http://www.bioinfo-scrounger.com 转载请注明出处

原文  http://www.bioinfo-scrounger.com/archives/622
正文到此结束
所属分类： Java 编程技术

热门推荐

相关文章

阿里云首购8折
说给你听

本文目录
随机标签
书籍教程
springboot-demo Java入门教程 bootstrap3 CSS Apache基础教程 php ionic 教程 Python mysql教程 eclipse Ubuntu VPS系统配置 AngularJS 教程 MongoDB教程 Struts2教程 Redis教程 Spring教程 Git教程 Jenkins进阶系列 openfire参考指南 Java设计模式 HBase教程 java-demo Maven教程 springcloud-demo hibernate教程 Docker 教程 memcached教程 Quartz指南 Hive教程 ANTLR教程 SpringCloud Ant教程 java实例教程 Hazelcast教程 XStream教程 Elastic-Job-Lite 深入浅出MyBatis ibaties教程 SVN教程 rabittmq教程 Hadoop教程 solr教程 WebService CXF学习 JPA教程 ActiveMQ中文指南 Java内存模型 dubbo教程 python3-demo Linux入门视频教程
近期评论
  1. 1 使用 Docker 打包和运行 Vue 应用
  2. 2 Dependency Check命令行方式扫描jar包的安全漏洞
  3. 3 maven使用Dependency-Check来扫描安全漏洞
  4. 4 使用sonarqube分析代码质量
  5. 5 使用screw来对比数据库表和字段差异
  6. 6 Spring Boot 实现防盗链功能
  7. 7 在Spring Boot中的实现国际化（i18n）
  8. 8 一个用于 SMB/CIFS 网络的 Java 库-JCIFS
  9. 9 Spring Boot 3.x嵌入MongoDB 进行测试
  10. 10 java如何解析和生成yaml文件？
网站信息
  • 文章总数：82,737 篇
  • 文件总数：284,312 个
  • 标签总数：2,416 个
  • 分类总数：85 个
  • 留言数量：2,559 条
  • 在线人数：663 人
  • 运行天数：4,439天
  • 最后更新：2024年12月22日15点
Loading...

其他链接

关于本站

本站定位：个人技术类博客

本站作用：写博客、记日志、闲聊扯淡鼓捣技术。

问题交流

[HBLOG]公众号 HBLOG HBLOG