重測序分析-使用GATK進行SNP和INDEL檢測

GATK 是 Broad 開發(fā)的用于二代測序數(shù)據(jù)變異檢測的軟件，后續(xù)推廣到
動植物研究中，是目前最廣泛使用的變異檢測軟件。本部分介紹如何使
用 GATK 進行 SNP 和 INDEL 檢測及過濾。

? 軟件：
GATK
? 文件準備：
基因組文件：genome.fasta
比對結(jié)果文件：S1.sort.markdup.bam
? 參考腳本：
第一步：創(chuàng)建 tmp 目錄

創(chuàng)建臨時文件目錄

$ mkdir tmp/
第二步：每個樣品進行 HaplotypeCaller 變異檢測, 基因組較小時，可以
整個基因組進行變異檢測；基因組較大時，建議分染色體進行變異檢測，
下面為分染色體運行的示例：

創(chuàng)建S1樣品分析結(jié)果目錄

$ mkdir S1

## S1樣品 chr1染色體 $ gatk --java-options "-Xmx10g -Djava.io.tmpdir=./tmp" \ #設置java參數(shù) HaplotypeCaller \ -R ../01.ref/genome.fasta \ # 參考基因組路徑 -I ../02.mapping//S1.sort.markdup.bam \ # bam文件路徑 -L chr1 \ # 指定染色體名稱 -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 變異檢測（分染色體）

## 生成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文件名稱 --genomicsdb-workspace-path genomeDB.chr1 \ # genomeDB目錄名稱 -L chr1 \ # 指定染色體 --reader-threads 1 \ # 線程數(shù) --batch-size 50 \ # 每次讀取文件個數(shù) --tmp-dir ./tmp \ 1>chr1.GenomicsDBImport.log 2>&1

## 群體變異檢測 $ 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 并過濾
### 提取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

### 過濾SNP（Filter列加標記） $ 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

### 提取過濾好的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

### 過濾INDEL（Filter列加標記） $ 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

### 提取過濾好的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