HGP

前言

在前面的一系列WGS文章中，我講述了很多基因數(shù)據(jù)分析的來(lái)龍去脈养晋，雖然許多同學(xué)覺(jué)得很有幫助囊拜，但是卻缺了一個(gè)重要的環(huán)節(jié)——沒(méi)有提供實(shí)際可用的數(shù)據(jù)來(lái)實(shí)戰(zhàn)完成具體的流程俱恶，不能得到直觀的體會(huì)。許多讀者也紛紛在后臺(tái)留言反饋這個(gè)問(wèn)題庸毫，特別是在我寫(xiě)了如何入門(mén)生物信息學(xué)的文章之后情況尤甚仔拟。所以我決定要再寫(xiě)一篇文章來(lái)解決這個(gè)問(wèn)題，但碰巧今年年末事情稍多飒赃，在寫(xiě)作的過(guò)程中也曾多次中斷利花，直至今天才完成，各位久等了载佳。本次文章分為上下兩篇炒事，這是第一篇，都是WGS第四節(jié)的延伸蔫慧，本意是結(jié)合具體的數(shù)據(jù)挠乳，讓更多的人能夠 更好地理解整個(gè)WGS數(shù)據(jù)的分析和處理過(guò)程，我也結(jié)合自身的工作經(jīng)驗(yàn)給出一些做項(xiàng)目過(guò)程中的建議藕漱，以作參考欲侮，希望能夠?qū)δ阌懈嗟膸椭Ａ硗饫吡酉聛?lái)我將系統(tǒng)寫(xiě)一個(gè)關(guān)于全基因組關(guān)聯(lián)分析（GWAS）的文章，同時(shí)還會(huì)有更多全面而且緊扣前沿的技術(shù)文章分享出來(lái)刁俭。

那么橄仍，事不宜遲我們馬上開(kāi)始‰蛊荩考慮到實(shí)際的數(shù)據(jù)分析環(huán)境侮繁，我們只在Linux命令行終端（Terminal）中進(jìn)行，執(zhí)行步驟都寫(xiě)在shell如孝，因此不會(huì)有窗口式的操作（使用Mac OS的同學(xué)可以使用Mac自帶的Terminal宪哩，與Linux操作一致），在這篇文章中我們會(huì)用到以下幾個(gè)工具第晰。

• sratoolkit
• bgzip
• tabix
• bwa
• samtools
• GATK 4.0

這些軟件都可以在github上找到（包括GATK）锁孟，需要各位自行安裝。這里補(bǔ)充一句茁瘦，目前GATK4.0的正式版本已經(jīng)發(fā)布品抽，它的使用方式與之前相比有著一些差異（變得更加簡(jiǎn)單，功能也更加豐富了）甜熔，增加了結(jié)構(gòu)性變異檢測(cè)和很多Spark圆恤、Cloud-Only的功能，并集成了Mutect2和picard的所有功能（以及其他很多有用的工具）腔稀，這為我們減少了許多額外的工具盆昙，更加有利于流程的構(gòu)建和維護(hù)羽历，4.0之后的GATK是一個(gè)新的篇章，大家最好是掌握這一個(gè)版本淡喜！另外秕磷，3.x的版本貌似也已經(jīng)不提供下載通道了，如果你還想使用3.x的話可以在公眾號(hào)后臺(tái)回復(fù)“GATK3”拆火，我為你準(zhǔn)備了一個(gè)GATK官方3.7的版本跳夭。我們這里則使用最新的4.0版本。

項(xiàng)目目錄結(jié)構(gòu)

清晰的目錄結(jié)構(gòu)是管理眾多項(xiàng)目的有效途徑们镜，經(jīng)久不忘币叹，隨時(shí)可查。雖然看起來(lái)有些原始模狭，但在Linux終端下面颈抚，我目前還沒(méi)有發(fā)現(xiàn)更好的文件管理辦法。這個(gè)項(xiàng)目的目錄結(jié)構(gòu)嚼鹉，我的建議是按照時(shí)間+項(xiàng)目的規(guī)則來(lái)命名贩汉，下面是我的目錄結(jié)構(gòu)：

./201802_wgs_practice/
├── bin
├── input
└── output

頂層的項(xiàng)目名就是20180203_wgs_practice，下面有三個(gè)主目錄：

• input：存儲(chǔ)所有輸入數(shù)據(jù)
• output：存儲(chǔ)所有輸出數(shù)據(jù)
• bin：存放所有執(zhí)行程序和代碼

output只存放結(jié)果數(shù)據(jù)锚赤，它是由input和bin中的數(shù)據(jù)和程序流程生成的匹舞。這樣做的好處是層次分明，流程邏輯清楚线脚，數(shù)據(jù)互不干擾赐稽。

使用E.coli K12完成比對(duì)和變異檢測(cè)

人類(lèi)基因組數(shù)據(jù)很大，參考序列長(zhǎng)度是3Gb浑侥。而一個(gè)人的高深度測(cè)序數(shù)據(jù)往往是這個(gè)數(shù)字的30倍——100Gb姊舵。如果直接用這樣的數(shù)據(jù)來(lái)完成本文的分析，那么許多同學(xué)需要下載大量的原始數(shù)據(jù)寓落。除了下載時(shí)間很長(zhǎng)之外括丁，如果沒(méi)有合適的集群，只是在自己的桌面電腦上干這樣的事情伶选，那么硬盤(pán)空間也將很快不夠用史飞。而且，要在單機(jī)電腦上完成這樣一個(gè)高深度WGS數(shù)據(jù)的分析考蕾，處理對(duì)機(jī)器性能有要求之后祸憋，跑起來(lái)也需要連續(xù)花上差不多140個(gè)小時(shí)——相信大家都等不起呀。

我已經(jīng)等不及了

因此肖卧，為了解決這個(gè)問(wèn)題蚯窥，我找了E.coli K12（一種實(shí)驗(yàn)用的大腸桿菌）的數(shù)據(jù)作為代替，用來(lái)演示 數(shù)據(jù)比對(duì)和變異檢測(cè)這兩個(gè)最消耗計(jì)算資源和存儲(chǔ)空間的步驟。E.coli K12的特點(diǎn)是數(shù)據(jù)很小拦赠，它的基因組長(zhǎng)度只有4.6Mb巍沙，很適合大家用來(lái)快速學(xué)習(xí)WGS的數(shù)據(jù)分析，遇到人類(lèi)的數(shù)據(jù)時(shí)荷鼠，再做替換就行了句携。

下載E.coli K12的參考基因組序列

熟悉的同學(xué)應(yīng)該第一時(shí)間能夠知道，這些物種的基因組參考序列都可以在NCBI上獲取允乐，我們這里也是一樣矮嫉，可以在NCBI網(wǎng)站上直接搜索這個(gè)序列，為了簡(jiǎn)化步驟牍疏，我直接給出E.coli K12參考序列的ftp地址給大家下載之用：

ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz  

你可以在Linux（或者M(jìn)ac OSX）命令行上直接使用wget蠢笋，將這個(gè)fasta下載下來(lái)，由于它很小鳞陨，所以幾秒之后我們就可以得到這個(gè)fasta序列昨寞。

$ wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz 

為了接下來(lái)表達(dá)上的清晰和操作上的方便，我們使用bgzip將這個(gè)序列文件進(jìn)行解壓并把名字重命名為E.coli_K12_MG1655.fa厦滤，這樣就一目了然了援岩。

$ gzip -dc GCF_000005845.2_ASM584v2_genomic.fna.gz > E.coli_K12_MG1655.fa

E.coli K12只有一條完整的染色體，你打開(kāi)文件后將會(huì)看到和我一樣的內(nèi)容：

>NC_000913.3 Escherichia coli str. K-12 substr. MG1655, complete genome
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGCTTCTGAACTG
GTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGGTCACTAAATACTTTAACCAATATAGGCATAGCGCACAGAC

接著掏导，我們用samtools為它創(chuàng)建一個(gè)索引享怀，這是為方便其他數(shù)據(jù)分析工具（比如GATK）能夠快速地獲取fasta上的任何序列做準(zhǔn)備。

$ /Tools/common/bin/samtools faidx E.coli_K12_MG1655.fa

這時(shí)會(huì)生成一份E.coli_K12_MG1655.fa.fai文件趟咆。除了方便其他工具之外凹蜈，我們可以通過(guò)這樣的索引來(lái)獲取fasta文件中任意位置的序列或者任意完整的染色體序列∪绦ィ可以很方便地完成對(duì)參考序列（或者任意fasta文件）特定區(qū)域序列的提取。舉個(gè)例子：

$ samtools faidx E.coli_K12_MG1655.fa NC_000913.3:1000000-1000200

我們就獲得了E.coli K12參考序列上的這一段序列：

>NC_000913.3:1000000-1000200
GTGTCAGCTTTCGTGGTGTGCAGCTGGCGTCAGATGACAACATGCTGCCAGACAGCCTGA
AAGGGTTTGCGCCTGTGGTGCGTGGTATCGCCAAAAGCAATGCCCAGATAACGATTAAGC
AAAATGGTTACACCATTTACCAAACTTATGTATCGCCTGGTGCTTTTGAAATTAGTGATC
TCTATTCCACGTCGTCGAGCG

這個(gè)小技巧在特定的時(shí)候非常實(shí)用履植。

下載E.coli K12的測(cè)序數(shù)據(jù)

基因組參考序列準(zhǔn)備好之后计雌，接下來(lái)我們需要下載它的測(cè)序數(shù)據(jù)。E.coli K12作為一種供研究使用的模式生物玫霎，自然已經(jīng)有許多的測(cè)序數(shù)據(jù)在NCBI上了凿滤，在這里我們選擇了其中的1個(gè)數(shù)據(jù)——SRR1770413。這個(gè)數(shù)據(jù)來(lái)自Illumina MiSeq測(cè)序平臺(tái)（不用擔(dān)心平臺(tái)的事情）庶近，read長(zhǎng)度是300bp翁脆，測(cè)序類(lèi)型Pair-End（沒(méi)了解過(guò)PE read同學(xué)可以參考我前面WGS系列的第四節(jié)文章）。你可以在NCBI上直接搜到：這里

在NCBI給出的信息頁(yè)面中鼻种，我們可以清楚地看到這個(gè)數(shù)據(jù)的大蟹捶（如下圖）——差不多200MB，一般家庭網(wǎng)速也能夠較快下載完成。

下載測(cè)序數(shù)據(jù)

從NCBI上下載下來(lái)的測(cè)序數(shù)據(jù)罢缸，不是我們熟悉的fastq格式篙贸，而是SRA（一種NCBI自己設(shè)計(jì)的測(cè)序數(shù)據(jù)存儲(chǔ)格式，具較高的壓縮率）枫疆，我們需要對(duì)其進(jìn)行轉(zhuǎn)換爵川，下文詳述。現(xiàn)在我們先下載息楔，有兩個(gè)下載方式（我在這里告訴大家的方法同樣適用于其他類(lèi)型的數(shù)據(jù)）寝贡，第一個(gè)是如上面所說(shuō)搜索到SRR1770413這個(gè)數(shù)據(jù)的ftq地址，然后直接在命令行中執(zhí)行wget進(jìn)行下載值依，如下：

$ wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR177/SRR1770413/SRR1770413.sra

注意圃泡，下載下來(lái)的這個(gè)SRA文件雖然只有一份，但是里面其實(shí)存了read1和read2的測(cè)序數(shù)據(jù)鳞滨，我們要把它解出來(lái)洞焙，轉(zhuǎn)換成為我們所需的fastq格式。這個(gè)時(shí)候拯啦，我們就需要用到NCBI的官方工具包sratoolkit澡匪，大家下載對(duì)應(yīng)系統(tǒng)的版本，直接解壓之后就可以使用了褒链。

sratoolkit工具包

sratoolkit是一個(gè)工具包唁情，所有的執(zhí)行程序都在它解壓后的bin文件夾下，我們要把SRA轉(zhuǎn)換為fastq甫匹，只需直接通過(guò)工具包中的fastq-dump即可完成甸鸟。

$ /Tools/sratoolkit/2.8.2/bin/fastq-dump --split-files SRR1770413.sra

然后我們就會(huì)得到這個(gè)E.coli K12數(shù)據(jù)的read1和read2了：

SRR1770413_1.fastq
SRR1770413_2.fastq

另一個(gè)數(shù)據(jù)下載的方法就是用上面說(shuō)到的sratoolkit，我也比較推薦這個(gè)方法兵迅，操作很簡(jiǎn)單抢韭。同樣是使用fastq-dump（沒(méi)錯(cuò)，與上面一樣恍箭，區(qū)別在于下載的時(shí)候刻恭，輸入的是數(shù)據(jù)的SRA編號(hào)），它可以在我們下載的過(guò)程中就直接將SRA轉(zhuǎn)換為兩個(gè)fastq扯夭。

$ /Tools/sratoolkit/2.8.2/bin/fastq-dump --split-files SRR1770413

下載完成后鳍贾，我們最好用bgzip（不推薦gzip）將其壓縮為.gz文件，這樣可以節(jié)省空間交洗，而且也不會(huì)對(duì)接下來(lái)的數(shù)據(jù)分析產(chǎn)生影響骑科。

$ /Tools/common/bin/bgzip -f SRR1770413_1.fastq
$ /Tools/common/bin/bgzip -f SRR1770413_1.fastq 

至此，E.coli K12相關(guān)的數(shù)據(jù)我們就都準(zhǔn)備好了构拳，先看一眼我現(xiàn)在的目錄結(jié)構(gòu)咆爽。

./201802_wgs_practice/
├── bin
├── input
│   ├── E.coli
│       ├── fasta
│       │   ├── E.coli_K12_MG1655.fa
│       │   ├── E.coli_K12_MG1655.fa.fai
│       │   └── work.log.sh
│       └── fastq
│           ├── SRR1770413_1.fastq.gz
│           ├── SRR1770413_2.fastq.gz
│           └── work.log.sh   
├── output
└── work.log.sh

其中各個(gè)目錄下的 work.log.sh記錄了我在該目錄下的所有重要操作——這是我的個(gè)人習(xí)慣梁棠，目的是方便以后反查數(shù)據(jù)的需要。

數(shù)據(jù)準(zhǔn)備完畢之后伍掀，接下來(lái)就可以進(jìn)行具體的分析了掰茶。

質(zhì)控

質(zhì)控是必須做的，我們需要完整認(rèn)識(shí)原始的測(cè)序數(shù)據(jù)質(zhì)量到底如何蜜笤，該步驟不能省略濒蒋。我專(zhuān)門(mén)為此單獨(dú)寫(xiě)了一篇文章（WGS系列第三節(jié)），在正式的數(shù)據(jù)分析過(guò)程中把兔，大家可以參考它來(lái)完成數(shù)據(jù)的質(zhì)控嫉柴，然后再進(jìn)行接下來(lái)的分析僵闯。本篇文章為了控制篇幅和盡可能扣住核心內(nèi)容碳柱，就不再對(duì)此深入展開(kāi)禀倔，大家如果碰到問(wèn)題，可以到在后臺(tái)留言缕贡，或者加入我的交流圈（知識(shí)星球：解螺旋技術(shù)交流圈）和更多有經(jīng)驗(yàn)的人一起交流翁授。

比對(duì)

首先是比對(duì)。所謂比對(duì)就是把測(cè)序數(shù)據(jù)定位到參考基因組上晾咪，確定每一個(gè)read在基因組中的位置收擦。這里，我們依然用目前使用最廣的BWA來(lái)完成這個(gè)工作谍倦。在正式比對(duì)之前塞赂，需要先為參考序列構(gòu)建BWA比對(duì)所需的FM-index（比對(duì)索引）。

$ /Tools/common/bin/bwa index E.coli_K12_MG1655.fa

由于這個(gè)序列很短昼蛀，只需幾秒就可以完成這個(gè)索引文件的構(gòu)建（對(duì)于人類(lèi)基因組則需要3個(gè)小時(shí)的時(shí)間）宴猾。創(chuàng)建完畢之后，將多出5份以E.coli_K12_MG1655.fa為前綴的序列索引文件叼旋。

E.coli_K12_MG1655.fa.amb
E.coli_K12_MG1655.fa.ann
E.coli_K12_MG1655.fa.bwt
E.coli_K12_MG1655.fa.pac
E.coli_K12_MG1655.fa.sa

現(xiàn)在我們使用bwa完成比對(duì)仇哆，用samtools完成BAM格式轉(zhuǎn)換、排序并標(biāo)記PCR重復(fù)序列夫植。步驟分解如下：

#1 比對(duì)
time /Tools/common/bin/bwa mem -t 4 -R '@RG\tID:foo\tPL:illumina\tSM:E.coli_K12' /Project/201802_wgs_practice/input/E.coli/fasta/E.coli_K12_MG1655.fa /Project/201802_wgs_practice/input/E.coli/fastq/SRR1770413_1.fastq.gz /Project/201802_wgs_practice/input/E.coli/fastq/SRR1770413_2.fastq.gz | /Tools/common/bin/samtools view -Sb - > /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam && echo "** bwa mapping done **"

#2 排序
time /Tools/common/bin/samtools sort -@ 4 -m 4G -O bam -o /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.bam /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam && echo "** BAM sort done"

rm -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam

#3 標(biāo)記PCR重復(fù)
time /Tools/common/bin/gatk/4.0.1.2/gatk MarkDuplicates -I /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.bam -O /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup.bam -M /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup_metrics.txt && echo "** markdup done **"

#4 刪除不必要文件(可選)
rm -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam
rm -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.bam

#5 創(chuàng)建比對(duì)索引文件
time /Tools/common/bin/samtools index /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup.bam && echo "** index done **"

從上面的命令大家也可以看到税产，我嚴(yán)格按照上文提到的項(xiàng)目目錄規(guī)范來(lái)執(zhí)行（步驟中涉及到的數(shù)據(jù)路徑也都盡可能使用全路徑），這個(gè)比對(duì)的shell存放在bin目錄下偷崩，名稱(chēng)是bwa_and_markdup.sh，從名稱(chēng)也能夠一眼可以看出是要做什么的撞羽。

簡(jiǎn)單解釋一下這個(gè)shell所做的事情：首先利用bwa mem比對(duì)模塊將E.coli K12質(zhì)控后的測(cè)序數(shù)據(jù)定位到其參考基因組上（我們這里設(shè)置了4個(gè)線程來(lái)完成比對(duì)阐斜，根據(jù)電腦性能可以適當(dāng)調(diào)大），同時(shí)通過(guò)管道（'|' 操作符）將比對(duì)數(shù)據(jù)流引到samtools轉(zhuǎn)換為BAM格式（SAM的二進(jìn)制壓縮格式）诀紊，然后重定向('>'操作符)輸出到文件中保存下來(lái)谒出。

-R 設(shè)置Read Group信息，雖然我在以前的文章中已經(jīng)反復(fù)強(qiáng)調(diào)過(guò)它的重要性，但這里還是再說(shuō)一次笤喳，它是read數(shù)據(jù)的組別標(biāo)識(shí)为居，并且其中的ID，PL和SM信息在正式的項(xiàng)目中是不能缺少的(如果樣本包含多個(gè)測(cè)序文庫(kù)的話杀狡，LB信息也不要省略)蒙畴，另外由于考慮到與GATK的兼容關(guān)系，PL（測(cè)序平臺(tái)）信息不能隨意指定呜象，必須是：ILLUMINA膳凝，SLX，SOLEXA恭陡，SOLID蹬音，454，LS454休玩，COMPLETE著淆，PACBIO，IONTORRENT拴疤，CAPILLARY永部，HELICOS或UNKNOWN這12個(gè)中的一個(gè)。

接著用samtools對(duì)原始的比對(duì)結(jié)果按照參考序列位置從小到大進(jìn)行排序（同樣是4個(gè)線程）遥赚，只有這個(gè)步驟完成之后才可以繼續(xù)往下扬舒。

然后，我們使用GATK標(biāo)記出排完序的數(shù)據(jù)中的PCR重復(fù)序列凫佛。這個(gè)步驟完成后讲坎，如無(wú)特殊需要，我們就可以直接刪除前面那兩個(gè)BAM文件了（原始比對(duì)結(jié)果和排序后的結(jié)果）——后續(xù)幾乎不會(huì)再用到那兩份文件了愧薛。關(guān)于標(biāo)記PCR重復(fù)序列的操作比較簡(jiǎn)單晨炕，不再細(xì)說(shuō)（如果希望了解更多有關(guān)重復(fù)序列特征的信息可以回看WGS系列第四節(jié)中的內(nèi)容）。

最后毫炉，我們?cè)儆胹amtools為E_coli_K12.sorted.markdup.bam創(chuàng)建索引瓮栗。我認(rèn)為不論是否有后續(xù)分析，為BAM文件創(chuàng)建索引應(yīng)該作為一個(gè)常規(guī)步驟瞄勾，它可以讓我們快速地訪問(wèn)基因組上任意位置的比對(duì)情況费奸，這一點(diǎn)非常有助于我們隨時(shí)了解數(shù)據(jù)。

至于每個(gè)步驟最前面的time进陡，則是用于記錄執(zhí)行時(shí)間的愿阐，有助于我們清楚地知道每一個(gè)分析過(guò)程都花了多少時(shí)間，當(dāng)需要優(yōu)化流程的時(shí)候這個(gè)信息會(huì)很有用趾疚。

變異檢測(cè)

接下來(lái)是用GATK完成變異檢測(cè)缨历。但在開(kāi)始之前之前我們還需要先為E.coli K12的參考序列生成一個(gè).dict文件以蕴，這可以通過(guò)調(diào)用CreateSequenceDictonary模塊來(lái)完成(這是原來(lái)picard的功能)。

$ /Tools/common/bin/gatk/4.0.1.2/gatk CreateSequenceDictionary -R E.coli_K12_MG1655.fa -O E.coli_K12_MG1655.dict && echo "** dict done **"

唯一需要注意的是.dict文件的名字前綴需要和fasta的一樣辛孵，并跟它在同一個(gè)路徑下丛肮，這樣GATK才能夠找到。

OK魄缚，現(xiàn)在我們就可以進(jìn)行變異檢測(cè)了宝与，同樣使用GATK 4.0的HaplotypeCaller模塊來(lái)完成。由于我們只有一個(gè)樣本鲜滩，要完成這個(gè)工作其實(shí)很簡(jiǎn)單伴鳖，直接輸入比對(duì)文件和參考序列就行了，但是考慮到實(shí)際的情況徙硅，我想告訴大家一個(gè)更好的方式（雖然這會(huì)多花些時(shí)間）榜聂，就是：先為每個(gè)樣本生成一個(gè)GVCF，然后再用GenotypeGVCFs對(duì)這些GVCF進(jìn)行joint calling嗓蘑，如下 须肆，我把命令都寫(xiě)在gatk.sh中，并執(zhí)行桩皿。

#1 生成中間文件gvcf
time /Tools/common/bin/gatk/4.0.1.2/gatk HaplotypeCaller \
  -R /Project/201802_wgs_practice/input/E.coli/fasta/E.coli_K12_MG1655.fa \
  --emit-ref-confidence GVCF \
  -I /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup.bam \
  -O /Project/201802_wgs_practice/output/E.coli/E_coli_K12.g.vcf && echo "** gvcf done **"

#2 通過(guò)gvcf檢測(cè)變異
time /Tools/common/bin/gatk/4.0.1.2/gatk GenotypeGVCFs \
  -R /Project/201802_wgs_practice/input/E.coli/fasta/E.coli_K12_MG1655.fa \
  -V /Project/201802_wgs_practice/output/E.coli/E_coli_K12.g.vcf \
  -O /Project/201802_wgs_practice/output/E.coli/E_coli_K12.vcf && echo "** vcf done **"

很快我們就獲得了E.coli K12這個(gè)樣本初步的變異結(jié)果——E_coli_K12.vcf豌汇。之所以非要分成兩個(gè)步驟，是因?yàn)槲蚁虢璐烁嬖V大家泄隔，變異檢測(cè)不是一個(gè)樣本的事情拒贱，有越多的同類(lèi)樣本放在一起joint calling結(jié)果將會(huì)越準(zhǔn)確，而如果樣本足夠多的話佛嬉，在低測(cè)序深度的情況下也同樣可以獲得完整并且準(zhǔn)確的結(jié)果逻澳，而這樣的分步方式是應(yīng)對(duì)多樣本的好方法。

最后暖呕，我們用bgzip對(duì)這個(gè)VCF進(jìn)行壓縮斜做，并用tabix為它構(gòu)建索引，方便以后的分析湾揽。

#1 壓縮 
time /Tools/common/bin/bgzip -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.vcf

#2 構(gòu)建tabix索引
time /Tools/common/bin/tabix -p vcf /Project/201802_wgs_practice/output/E.coli/E_coli_K12.vcf.gz

bgzip壓縮完成之后瓤逼，原來(lái)的VCF文件會(huì)被自動(dòng)刪除。

為了保持一致库物，現(xiàn)在再看一下完成到這里之后我們的目錄長(zhǎng)什么樣了霸旗，供大家對(duì)照。

./201802_wgs_practice/
├── bin
│   ├── bwa_and_markdup.sh
│   └── gatk.sh
├── input
│   └── E.coli
│       ├── fasta
│       │   ├── E.coli_K12_MG1655.dict
│       │   ├── E.coli_K12_MG1655.fa
│       │   ├── E.coli_K12_MG1655.fa.amb
│       │   ├── E.coli_K12_MG1655.fa.ann
│       │   ├── E.coli_K12_MG1655.fa.bwt
│       │   ├── E.coli_K12_MG1655.fa.fai
│       │   ├── E.coli_K12_MG1655.fa.pac
│       │   ├── E.coli_K12_MG1655.fa.sa
│       │   └── work.log.sh
│       └── fastq
│           ├── SRR1770413_1.fastq.gz
│           ├── SRR1770413_2.fastq.gz
│           └── work.log.sh
├── output
│   └── E.coli
│       ├── E_coli_K12.g.vcf
│       ├── E_coli_K12.g.vcf.idx
│       ├── E_coli_K12.sorted.markdup.bam
│       ├── E_coli_K12.sorted.markdup.bam.bai
│       ├── E_coli_K12.sorted.markdup_metrics.txt
│       ├── E_coli_K12.vcf
│       └── E_coli_K12.vcf.idx
└── work.log.sh

如果大家仔細(xì)看過(guò)WGS系列第四節(jié)的話戚揭，會(huì)發(fā)現(xiàn)我這里缺少了兩個(gè)步驟：重比對(duì)和BQSR定硝。沒(méi)有執(zhí)行BQSR是因?yàn)镋.coli K12沒(méi)有那些必須的known變異集（或者有但我沒(méi)找到），所以無(wú)法進(jìn)行毫目；但沒(méi)有重比對(duì)蔬啡，則是因?yàn)槲以贕ATK 4.0中沒(méi)發(fā)現(xiàn)IndelRealigner這個(gè)功能，雖然我們使用GATK HaplotypeCaller或者M(jìn)utect2的話確實(shí)可以省略這個(gè)步驟镀虐，但如果是其他軟件來(lái)進(jìn)行變異檢測(cè)那么該步驟依然十分重要箱蟆，我目前不太清楚為何GATK 4.0沒(méi)有將這個(gè)功能單獨(dú)分離出來(lái)。

后面要談到的就是變異的質(zhì)控了刮便。很遺憾我們這個(gè)E.coli K12的變異結(jié)果并不適合通過(guò)VQSR來(lái)進(jìn)行過(guò)濾空猜，原因上面也提到了一些，它不像人類(lèi)的基因組數(shù)據(jù)恨旱，有著一套適合用來(lái)訓(xùn)練過(guò)濾模型的已知變異集（dbSNP辈毯，1000G，Hapmap和omini等）搜贤。其實(shí)這種情況有時(shí)候我們?cè)诠ぷ髦幸矔?huì)碰到谆沃，比如有些捕獲測(cè)序（Panel測(cè)序數(shù)據(jù)，甚至外顯子測(cè)序）的數(shù)據(jù)仪芒，由于它的區(qū)域較小唁影，獲得的變異也不多，導(dǎo)致最終沒(méi)法滿(mǎn)足VQSR進(jìn)行模型訓(xùn)練時(shí)所需的最低變異數(shù)要求掂名，那時(shí)你也不能通過(guò)這個(gè)方式協(xié)助變異質(zhì)控据沈。那么碰到這種情況的時(shí)候該怎么辦？我將這部分的內(nèi)容放在了下一篇文章中饺蔑，在那里我們?cè)賮?lái)討論這個(gè)問(wèn)題锌介。我也會(huì)告訴大家變異質(zhì)控的基本邏輯，而不是簡(jiǎn)單羅列一個(gè)命令猾警，同時(shí)也會(huì)再用NA12878這個(gè)人的數(shù)據(jù)來(lái)進(jìn)一步告訴大家如何比較和評(píng)估變異結(jié)果孔祸。

小結(jié)

至此，這個(gè)篇文章的上半部分就到此為止了肿嘲。除了那些重要的內(nèi)容之外融击，在上文中，你會(huì)看到我反復(fù)提到了創(chuàng)建“索引”這個(gè)事情雳窟，比如為fasta尊浪，為BAM，為VCF封救。我為什么非要反復(fù)強(qiáng)調(diào)這個(gè)事情不可呢拇涤？因?yàn)槲野l(fā)現(xiàn)許多初學(xué)者并不知道索引的作用，當(dāng)被問(wèn)到如何從巨大的比對(duì)文件或者變異文件中提取某個(gè)信息時(shí)誉结，總是要走彎路——努力寫(xiě)程序去提取鹅士，既慢又費(fèi)力，結(jié)果還不一定好惩坑，甚至有些有一定經(jīng)驗(yàn)的同學(xué)也不知道使用bgzip和tabix的好處掉盅，因此我才反復(fù)在文章里提及也拜。

歡迎關(guān)注我的公眾號(hào)：堿基礦工（helixminer），更及時(shí)了解更多信息