重測(cè)序分析-使用GATK進(jìn)行SNP和INDEL檢測(cè)
GATK 是 Broad 開(kāi)發(fā)的用于二代測(cè)序數(shù)據(jù)變異檢測(cè)的軟件，后續(xù)推廣到
動(dòng)植物研究中牙丽，是目前最廣泛使用的變異檢測(cè)軟件。本部分介紹如何使
用 GATK 進(jìn)行 SNP 和 INDEL 檢測(cè)及過(guò)濾。
? 軟件：
GATK
? 文件準(zhǔn)備：
基因組文件：genome.fasta
比對(duì)結(jié)果文件：S1.sort.markdup.bam
? 參考腳本：
第一步：創(chuàng)建 tmp 目錄
創(chuàng)建臨時(shí)文件目錄
$ mkdir tmp/
第二步：每個(gè)樣品進(jìn)行 HaplotypeCaller 變異檢測(cè), 基因組較小時(shí)劲厌，可以
整個(gè)基因組進(jìn)行變異檢測(cè)懊渡；基因組較大時(shí)刽射，建議分染色體進(jìn)行變異檢測(cè)军拟，
下面為分染色體運(yùn)行的示例：
創(chuàng)建S1樣品分析結(jié)果目錄
$ mkdir S1
## S1樣品 chr1染色體 $ gatk --java-options "-Xmx10g -Djava.io.tmpdir=./tmp" \ #設(shè)置java參數(shù) HaplotypeCaller \ -R ../01.ref/genome.fasta \ # 參考基因組路徑 -I ../02.mapping//S1.sort.markdup.bam \ # bam文件路徑 -L chr1 \ # 指定染色體名稱(chēng) -ERC GVCF \ # 輸出GVCF文件 -O S1/S1.chr1.g.vcf.gz \ # 指定輸出結(jié)果 1>S1/S1.chr1.HC.log 2>&1 # 保留日志文件
## S1樣品 chr2染色體 $ gatk --java-options "-Xmx10g -Djava.io.tmpdir=./tmp" \ HaplotypeCaller \ -R ../01.ref/genome.fasta \ -I ../02.mapping//S1.sort.markdup.bam \ -L chr2 \ -ERC GVCF \ -O S1/S1.chr2.g.vcf.gz \ 1>S1/S1.chr2.HC.log 2>&1
第三步：多樣本 gvcf 變異檢測(cè)（分染色體）
## 生成chr1染色體map文件 $ ls ./*/*.chr1.g.vcf.gz | awk -F"/" '{print $2"\t"$0}' > gvcf.chr1.map $ cat gvcf.chr1.map S1 ./S1/S1.chr1.g.vcf.gz S2 ./S2/S2.chr1.g.vcf.gz
## 生成chr1染色體的genomeDB $ gatk --java-options "-Xmx10g -Djava.io.tmpdir=./tmp \ -DGATK_STACKTRACE_ON_USER_EXCEPTION=true" \ GenomicsDBImport \ --sample-name-map gvcf.chr1.map \ # 指定map文件名稱(chēng) --genomicsdb-workspace-path genomeDB.chr1 \ # genomeDB目錄名稱(chēng) -L chr1 \ # 指定染色體 --reader-threads 1 \ # 線程數(shù) --batch-size 50 \ # 每次讀取文件個(gè)數(shù) --tmp-dir ./tmp \ 1>chr1.GenomicsDBImport.log 2>&1
## 群體變異檢測(cè) $ gatk --java-options "-Xmx10g -Djava.io.tmpdir=./tmp" \ GenotypeGVCFs \ -R ../01.ref/genome.fasta \ -V gendb://genomeDB.chr1 \ -O chr1.raw.vcf.gz \ 1>chr1.GenotypeGVCFs.log 2>&1
第四步：分染色體 vcf 文件合并成全基因組 vcf 文件
## 生成vcf列表文件 $ awk '{print $1".raw.vcf.gz"}' chr.list > raw_vcf.list $ cat raw_vcf.list chr1.raw.vcf.gz chr2.raw.vcf.gz $ gatk --java-options "-Xmx10g -Djava.io.tmpdir=./tmp" \ MergeVcfs \ -I raw_vcf.list \ -O all.merge_raw.vcf
第五步：分別提取 SNP 和 INDEL 并過(guò)濾
### 提取SNP $ gatk --java-options "-Xmx4g -Djava.io.tmpdir=./tmp" \ SelectVariants \ -R ../01.ref/genome.fasta \ -V all.merge_raw.vcf \ --select-type SNP \ -O all.raw.snp.vcf
### 過(guò)濾SNP（Filter列加標(biāo)記） $ gatk --java-options "-Xmx4g -Djava.io.tmpdir=./tmp" \ VariantFiltration \ -R ../01.ref/genome.fasta \ -V all.raw.snp.vcf \ --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 \ || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \ --filter-name 'SNP_filter' \ -O all.filter.snp.vcf
### 提取過(guò)濾好的SNP $ gatk --java-options "-Xmx4g -Djava.io.tmpdir=./tmp" \ SelectVariants \ -R ../01.ref/genome.fasta \ -V all.filter.snp.vcf \ --exclude-filtered \ -O all.filtered.snp.vcf
### 提取INDEL $ gatk --java-options "-Xmx4g -Djava.io.tmpdir=./tmp" \ SelectVariants \ -R ../01.ref/genome.fasta \ -V all.merge_raw.vcf \ --select-type INDEL \ -O all.raw.indel.vcf
### 過(guò)濾INDEL（Filter列加標(biāo)記） $ gatk --java-options "-Xmx4g -Djava.io.tmpdir=./tmp" \ VariantFiltration \ -R ../01.ref/genome.fasta \ -V all.raw.indel.vcf \ --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 \ || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \ --filter-name 'INDEL_filter' \ -O all.filter.indel.vcf
### 提取過(guò)濾好的INDEL $ gatk --java-options "-Xmx4g -Djava.io.tmpdir=./tmp" \ SelectVariants \ -R ../01.ref/genome.fasta \ -V all.filter.indel.vcf \ --exclude-filtered \ -O all.filtered.indel.vcf
