家和建材廣場：聯(lián)系電話-0373-7621288

位置：河南省新鄉(xiāng)市延津縣人民路與西環(huán)路交叉口向東100米路北私痹。

有家和幸福生活歼冰！

GATK4.0全基因組和全外顯子組分析實戰(zhàn)

文章來源：企鵝號 - 生信知識

前言

GATK是目前業(yè)內(nèi)最權(quán)威获茬、使用最廣的基因數(shù)據(jù)變異檢測工具。相比samtools + bcftools call SNP/Indel谴忧，GATK更加精確哎垦，當然代價是流程復(fù)雜且耗時長。目前已經(jīng)更新到GATK4嫂侍，GATK3在官網(wǎng)已沒有下載通道但其使用仍然很廣泛儿捧。GATK4在核心算法層面并沒太多的修改，但參數(shù)設(shè)置還是有些改變的挑宠，并且取消了RealignerTargetCreator菲盾、IndelRealigner，應(yīng)該是HaplotypeCaller繼承了這部分功能各淀。GATK4使用了新的設(shè)計模式懒鉴，做了很多功能的整合，已經(jīng)把picard完全整合碎浇。本文使用的是GATK4.0.2.0临谱，參考了其他人編寫的GATK3x和GATK4x教程，對全基因組和全外顯子組的每個步驟都進行了講解并放上腳本奴璃，但讀者還是要略作修改才能本地運行成功悉默，所以具有基本的生信知識和編程、Linux技能才能更好的學(xué)習(xí)這篇教程苟穆，下面開始正文抄课。

軟件

需要fastqc、fastp鞭缭、BWA剖膳、samtools、GATK岭辣、annovar（需要學(xué)術(shù)郵箱才能注冊下載）

數(shù)據(jù)質(zhì)控

這部分不贅述吱晒，公司交給你的數(shù)據(jù)肯定質(zhì)量很高，否則不會交付給你沦童，但是為了放心還是要檢查下的仑濒。先用fastqc

wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.7.zip

unzip fastqc_v0.11.7.zip

cd FastQC

chmod755 fastqc

使用：

fastqc wes_1.fq.gz -o fastqc_out_dir/? &&? fastqcwes_2.fq.gz -o fastqc_out_dir/

數(shù)據(jù)過濾使用fastp

安裝：wget http://opengene.org/fastp/fastp

chmod 755 fastp

使用：./fastp -i in.R1.fq.gz -o out.R1.fq.gz -I in.R2.fq.gz -O out.R2.fq.gz

比對

第一步下載參考基因組

for i in $(seq 1 22) X Y M;

do

wget? http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr$.fa.gz &

Done

下載完成后解壓合并

gunzip *.gz

for i in $(seq 1 22) X Y M;

do cat chr$.fa >> hg19.fa;

Done

構(gòu)建索引

bwa index -a bwtsw -p hg19 hg19.fa? &

程序運行時間較長，建議使用nohup 或者 screen

最終生成文件? hg19.amb? hg19.ann? hg19.bwt? hg19.pac hg19.sa

bwa比對

-t偷遗，線程數(shù)墩瞳；

-M , -M 將 shorter split hits 標記為次優(yōu)，以兼容 Picard’s markDuplicates 軟件;

-R 接的是 Read Group的字符串信息氏豌，它是用來將比對的read進行分組的喉酌，這個信息對于我們后續(xù)對比對數(shù)據(jù)進行錯誤率分析和Mark duplicate時非常重要。不設(shè)置-R參數(shù)不會報錯，但使用GATK時是會報錯的泪电。

(1)ID般妙，這是Read Group的分組ID，一般設(shè)置為測序的lane ID

(2) PL相速，指的是所用的測序平臺

(3) SM碟渺，樣本ID

(4) LB，測序文庫的名字

這些信息設(shè)置好之后突诬，在RG字符串中要用制表符（\t）將它們分開苫拍。

上一步生成的SAM文件是文本文件，一般整個文件都非常巨大旺隙，因此绒极，為了有效節(jié)省磁盤空間，一般都會用samtools將它轉(zhuǎn)化為BAM文件（SAM的特殊二進制格式）催束，而且BAM會更加方便于后續(xù)的分析集峦。

samtools? view -b -S abc.sam > abc.bam ##由sam生成bam

samtools view -h abc.bam>abc.sam##由bam生成sam，-h代表是否帶上header

想看bam文件使用命令samtools view -h abc.bam | less

bam文件是一個重要的文件抠刺，每一列的意義都需要理解塔淤。

前五列分別為：reads名flag染色體位置比對質(zhì)量。第十列時read序列速妖，11列是堿基質(zhì)量值高蜂。第1，10罕容，11列可以提取出來還原成我們的測序數(shù)據(jù)fastq格式备恤。有人說第五列比對質(zhì)量值為0時代表該read可以比對到基因組多個位置，說法是否正確我不確定锦秒。

第二列的flag包含很多信息露泊，flag解釋網(wǎng)站 http://broadinstitute.github.io/picard/explain-flags.html ，打開此網(wǎng)站

輸入flag值便可以看到包含的信息旅择，其中read unmapped是4惭笑，mate unmapped是8，如果雙端reads都不能map到參考基因組上生真，flag值就包括4和8沉噩，和為12≈埃可以用命令samtools view –f 12 wes.bam > unmapped.sam 提取雙端reads都沒有比對成功的reads川蒙，samtools view –F 12 wes.bam > mapped.sam 代表過濾雙端都比對不上的reads。

sort 排序

BWA比對后輸出的bam文件是沒順序的长已，比對后得到的結(jié)果文件中畜眨，每一條記錄之間位置的先后順序是亂的昼牛，我們后續(xù)去重復(fù)等步驟都需要在比對記錄按照順序從小到大排序下來才能進行。

做類似分析的時候在文件名字將所做的關(guān)鍵操作包含進去康聂，因為這樣即使過了很長時間匾嘱，當你再去看這個文件的時候也能夠立刻知道當時對它做了什么。

bam信息統(tǒng)計

做到這一步需要對序列比對情況進行統(tǒng)計早抠，如果比對情況很差需要查找原因。

覆蓋度撬讽，深度等信息統(tǒng)計

覆蓋度和深度是我們關(guān)心的重要參數(shù)蕊连，如果是全外顯子組可以用picard（已經(jīng)整合到GATK4中）進行統(tǒng)計。

生成參考基因組的dict文件

$GATK CreateSequenceDictionary -R hg19.fa -O hg19.dict

生成interval

外顯子組是用試劑盒捕獲再進行測序游昼，不同試劑盒捕獲的區(qū)域不同甘苍，要下載相應(yīng)的包含捕獲區(qū)域bed文件，本文用的是安捷倫的捕獲區(qū)域的bed文件烘豌。

覆蓋度载庭，深度等信息統(tǒng)計

On targeted bases相對總bases的比例（PCT_USABLE_BASES_ON_BAIT）

On and near targeted bases相對總bases的比例（PCT_SELECTED_BASES）

MEAN_TARGET_COVERAGE平均覆蓋度

如果是全基因組wgs，運行以下命令廊佩，

生成的wgs.metrics包含多種信息囚聚，可自行查閱研究。

去除由于PCR擴增引起的重復(fù)reads

在NGS測序之前都需要先構(gòu)建測序文庫：通過物理（超聲）打斷或者化學(xué)試劑（酶切）切斷原始的DNA序列标锄，然后選擇特定長度范圍的序列去進行PCR擴增并上機測序顽铸。這個過程中產(chǎn)生的重復(fù)reads，增大了變異檢測結(jié)果的假陰率和假陽率Ａ匣省谓松！原因如下：

1.PCR反應(yīng)過程中也會帶來新的堿基錯誤。發(fā)生在前幾輪的PCR擴增發(fā)生的錯誤會在后續(xù)的PCR過程中擴大践剂，同樣帶來假的變異鬼譬；

2. PCR反應(yīng)可能會對包含某一個堿基的DNA模版擴增更加劇烈（這個現(xiàn)象稱為PCR Bias）

3.如果某個變異位點的變異堿基都是來自于PCR重復(fù)，而我們卻認為它深度足夠判斷是真的變異位點逊脯，這個結(jié)論其實有很大可能是假陽性优质。

利用picard標記重復(fù)序列

查看被標記重復(fù)的reads

samtools view –f 1024wes.sorted.MarkDuplicates.bam| less

為何是1024去查閱flag解釋網(wǎng)站。

下一步為wes.sorted.MarkDuplicates.bam創(chuàng)建索引文件男窟，它的作用能夠讓我們可以隨機訪問這個文件中的任意位置盆赤，而且后面的步驟也要求這個BAM文件一定要有索引.

samtools index wes.sorted.MarkDuplicates.bam

變異檢測

開始檢測變異前還是要做一些準備工作，首先是重新校正堿基質(zhì)量值（BQSR）歉眷。

變異檢測是一個極度依賴測序堿基質(zhì)量值牺六，因為這個質(zhì)量值是衡量我們測序出來的這個堿基到底有多正確的重要指標。它來自于測序圖像數(shù)據(jù)的base calling汗捡，因此淑际，基本上是由測序儀和測序系統(tǒng)來決定的畏纲，計算出來的堿基質(zhì)量值要么高于真實結(jié)果，要么低于真實結(jié)果春缕。

BQSR（Base Quality Score Recalibration）這個步驟就是為此而存在的盗胀，這一步非常重要。它主要是通過機器學(xué)習(xí)的方法構(gòu)建測序堿基的錯誤率模型锄贼，然后對這些堿基的質(zhì)量值進行相應(yīng)的調(diào)整票灰。

這里包含了兩個步驟：

第一步，BaseRecalibrator宅荤，這里計算出了所有需要進行重校正的read和特征值屑迂，然后把這些信息輸出為一份校準表文件（wes.recal_data.table）

第二步，ApplyBQSR冯键，這一步利用第一步得到的校準表文件（wes.recal_data.table）重新調(diào)整原來BAM文件中的堿基質(zhì)量值惹盼，并使用這個新的質(zhì)量值重新輸出一份新的BAM文件。

-R $GENOME -I wes.sorted.MarkDuplicates.bam \

--known-sites $hg19_VCF/dbsnp_138.hg19.vcf \

-O wes.recal_data.table

（注意：此步驟以及后面幾個步驟中外顯子數(shù)據(jù)要加上外顯子捕獲區(qū)域的bed文件惫确，并把-ip設(shè)為reads長手报，全基因組數(shù)據(jù)則不需要加-L 和 -ip）

$GATK --java-options "-Xmx8G -Djava.io.tmpdir=./" ApplyBQSR \

-R $GENOME -I wes.sorted.MarkDuplicates.bam \

-bqsr wes.recal_data.table \

-O wes.sorted.MarkDuplicates.BQSR.bam

變異檢測前確定bam文件是否符合GATK要求，運行

$GATK? ValidateSamFile -I wes.sorted.MarkDuplicates.BQSR.bam改化，如果顯示 no error掩蛤，則可以用HaplotypeCaller call SNP/Indel。

利用samtools為hg19.fa創(chuàng)建一個索引 samtools faidx hg19.fa

利用這個索引可以查看參考基因組任何位置所袁，如運行

利用HaplotypeCaller檢測突變還需利用許多數(shù)據(jù)庫盏档，這些數(shù)據(jù)庫包含以前研究過的突變位點，利用這些位點提高變異檢測的準確率燥爷。

運行上面命令下載所需數(shù)據(jù)庫蜈亩。注意下載文件均為壓縮文件，需要解壓才能使用前翎。

下面開始HaplotypeCaller突變檢測稚配，

HaplotypeCaller的應(yīng)用有兩種做法，區(qū)別在于是否生成中間文件gVCF：

（1）直接進行HaplotypeCaller港华，這適合于單樣本道川，只執(zhí)行一次HaplotypeCaller。如果多樣本立宜，你每增加一個樣本數(shù)據(jù)都需要重新運行這個HaplotypeCaller冒萄，而這個時候算法需要重新去讀取所有人的BAM文件，浪費大量時間橙数；

（2）每個樣本先各自生成gVCF尊流，然后再進行群體joint-genotype。gVCF全稱是genome VCF灯帮，是每個樣本用于變異檢測的中間文件崖技，格式類似于VCF逻住，它把joint-genotype過程中所需的所有信息都記錄在這里面，文件無論是大小還是數(shù)據(jù)量都遠遠小于原來的BAM文件迎献。這樣一旦新增加樣本也不需要再重新去讀取所有人的BAM文件了瞎访，只需為新樣本生成一份gVCF，然后重新執(zhí)行這個joint-genotype就行了吁恍。

第一種方法

-R $GENOME? -I wes.sorted.MarkDuplicates.BQSR.bam \

第二種方法

#1 生成中間文件gvcf

-R $GENOME --emit-ref-confidence GVCF \

-I wes.sorted.MarkDuplicates.BQSR.bam -D $hg19_VCF/dbsnp_138.hg19.vcf -L S31285117_Regions.bed -ip 90 -O wes.gvcf

#2 通過gvcf檢測變異

若有多個樣本的gvcf文件扒秸，運行$GATK CombineGVCFs? -V 1.gvcf –V 2.gvcf ……? -O final.gvcf ，再用final.gvcf運行下一步

變異質(zhì)控和過濾

質(zhì)控的含義和目的是指通過一定的標準冀瓦，最大可能地剔除假陽性的結(jié)果鸦采，并盡可能地保留最多的正確數(shù)據(jù)。

第一種方法 GATK VQSR咕幻，它通過機器學(xué)習(xí)的方法利用多個不同的數(shù)據(jù)特征訓(xùn)練一個模型（高斯混合模型）對變異數(shù)據(jù)進行質(zhì)控，使用VQSR需要具備以下兩個條件：

第一顶霞，需要一個精心準備的已知變異集肄程，它將作為訓(xùn)練質(zhì)控模型的真集。比如选浑，Hapmap蓝厌、OMNI，1000G和dbsnp等這些國際性項目的數(shù)據(jù)古徒，這些可以作為高質(zhì)量的已知變異集拓提。

第二，要求新檢測的結(jié)果中有足夠多的變異隧膘，不然VQSR在進行模型訓(xùn)練的時候會因為可用的變異位點數(shù)目不足而無法進行代态。適合全基因組分析。

第二種方法通過過濾指標過濾疹吃。

QualByDepth（QD）

FisherStrand (FS)

StrandOddsRatio (SOR)

RMSMappingQuality (MQ)

MappingQualityRankSumTest (MQRankSum)

ReadPosRankSumTest (ReadPosRankSum)

GATK VQSR

-resource hapmap,known=false,training=true,truth=true,prior=15.0:$hg19_VCF/hapmap_3.3.hg19.sites.vcf \

-resource omini,known=false,training=true,truth=false,prior=12.0:$hg19_VCF/1000G_omni2.5.hg19.sites.vcf \

-resource dbsnp,known=true,training=false,truth=false,prior=6.0:$hg19_VCF/dbsnp_138.hg19.vcf \

-an QD -an MQ -an MQRankSum -an ReadPosRankSum -an SOR -mode SNP \

$GATK ApplyVQSR -V wes.raw.vcf -O wes.snps.VQSR.vcf \

此方法要求新檢測的結(jié)果中有足夠多的變異蹦疑，不然VQSR在進行模型訓(xùn)練的時候會因為可用的變異位點數(shù)目不足而無法進行∪唬可能很多非人的物種在完成變異檢測之后沒法使用GATK VQSR的方法進行質(zhì)控歉摧，一些小panel、外顯子測序腔呜，由于最后的變異位點不夠叁温，也無法使用VQSR。全基因組分析或多個樣本的全外顯子組分析適合用此方法核畴。

通過過濾指標過濾

# 使用SelectVariants膝但，選出SNP

$GATK SelectVariants -select-type SNP -V wes.GATK.vcf -O wes.GATK.snp.vcf

# 為SNP作過濾

# 使用SelectVariants，選出Indel

$GATK SelectVariants -select-type INDEL -V wes.GATK.vcf -O wes.GATK.indel.vcf

# 為Indel作過濾

VCF文件

VCF文件是記錄突變信息的重要文件

前五列信息為

1. 染色體(Chromosome)

2. 起始位置(Start)

3. 結(jié)束位置(End)

4. 參考等位基因(Reference Allele)

5. 替代等位基因(Alternative Allele)

ANNOVAR注釋時主要也是利用前五列信息對數(shù)據(jù)庫進行比對膛檀，注釋變異锰镀。Info和Format信息同樣重要娘侍，比如DP代表測序深度，這些內(nèi)容的含義再vcf文件的開頭都有介紹泳炉，請仔細閱讀并理解相應(yīng)內(nèi)容的意義憾筏。

突變注釋

vcf文件中保存的突變位點不能直接使用，只有根據(jù)已有數(shù)據(jù)庫進行注釋花鹅，才能知道該位點的有何功能氧腰，是否與疾病相關(guān)以及其它信息。annovar是突變注釋的常用軟件刨肃，ANNOVAR是一個perl編寫的命令行工具古拴，能在安裝了perl解釋器的多種操作系統(tǒng)上執(zhí)行。允許多種輸入文件格式真友，包括最常被使用的VCF格式黄痪。輸出文件也有多種格式，包括注釋過的VCF文件盔然、用tab或者逗號分隔的txt文件桅打。ANNOVAR能快速注釋遺傳變異并預(yù)測其功能，這個軟件需要edu郵箱注冊才能下載愈案。

http://www.openbioinformatics.org/annovar/annovar_download_form.php

軟件下載后無需安裝直接使用挺尾，但要下載數(shù)據(jù)庫才能對突變位點進行注釋，annotate_variation.pl可以方便快捷下載數(shù)據(jù)庫站绪。

perl annotate_variation.pl -buildver hg19 -downdb -webfrom annovar refGene humandb/

# -buildver 表示version

# -downdb 下載數(shù)據(jù)庫的指令

# -webfrom annovar 從annovar提供的鏡像下載遭铺，不加此參數(shù)將尋找數(shù)據(jù)庫本身的源

# humandb/ 存放于humandb/目錄下

其它數(shù)據(jù)庫下載

perl annotate_variation.pl--buildver hg19 --downdb gwascatalog humandb/ &

perl annotate_variation.pl--buildver hg19 --downdb ljb26_all --webfrom annovarhumandb/ &

perl annotate_variation.pl--buildver hg19 --downdb esp6500siv2_ea --webfromannovar humandb/ &

perl annotate_variation.pl--buildver hg19 --downdb esp6500siv2_all --webfromannovar humandb/ &

perl annotate_variation.pl--buildver hg19 --downdb 1000g2014oct humandb/ &

perl annotate_variation.pl--buildver hg19 --downdb cytoBand humandb/ &

perl annotate_variation.pl--buildver hg19 --downdb avsift -webfrom annovarhumandb/ &

perl annotate_variation.pl--buildver hg19 --downdb snp138 humandb/ &

perl annotate_variation.pl--buildver hg19 --downdb genomicSuperDups humandb/ &

perl annotate_variation.pl--buildver hg19 --downdbphastConsElements46wayhumandb/ &

perl annotate_variation.pl--buildver hg19 --downdb tfbs humandb/

annovar有三種注釋方式Gene-based Annotation(基于基因的注釋)，Region-based Annotation（基于區(qū)域的注釋）恢准，F(xiàn)ilter-based Annotation（基于過濾的注釋）魂挂。看起來很復(fù)雜實際做起來很簡單馁筐，用table_annovar.pl進行注釋锰蓬，可一次性完成三種類型的注釋。

注釋的第一步是把要注釋的文件轉(zhuǎn)化成annovar需要的文件眯漩，

convert2annovar.pl #將多種格式轉(zhuǎn)為.avinput的程序

avinput文件由tab分割芹扭，最重要的地方為前5列，分別是:

1. 染色體(Chromosome)

2. 起始位置(Start)

3. 結(jié)束位置(End)

4. 參考等位基因(Reference Allele)

5. 替代等位基因(Alternative Allele)

annovar主要也是利用前五列信息對數(shù)據(jù)庫進行比對赦抖，注釋變異舱卡。

SNP注釋：

用table_annovar.pl進行注釋，可一次性完成三種類型的注釋队萤。

table_annovar.pl snp.avinput $humandb -buildver hg19 -out snpanno -remove -protocol refGene,cytoBand,genomicSuperDups,esp6500siv2_all -operation g,r,r,f -nastring . –csvout

# -buildver hg19 表示使用hg19版本

# -out snpanno 表示輸出文件的前綴為snpanno

# -remove 表示刪除注釋過程中的臨時文件

# -protocol 表示注釋使用的數(shù)據(jù)庫轮锥，用逗號隔開，且要注意順序

# -operation 表示對應(yīng)順序的數(shù)據(jù)庫的類型（g代表gene-based要尔、r代表region-based舍杜、f代表filter-based）新娜，用逗號隔開，注意順序

# -nastring . 表示用點號替代缺省的值

# -csvout 表示最后輸出.csv文件

Indel注釋同上

最后生成的注釋文件如下

總結(jié)

本教程基本按照GATK4 Best PracticesGermline short variant discovery (SNPs + Indels)既绩，GATK3與GATK4的分析思路也幾乎沒有變化概龄，除了GATK4取消了RealignerTargetCreator 和 IndelRealigner 這兩步，另外一些命令發(fā)生了變化饲握。在外顯子分析步驟中記得加上外顯子捕獲區(qū)域的bed文件（-L參數(shù)）私杜，這樣可以防止檢測出捕獲區(qū)域以外的突變，還可以顯著提高運行速度