自己在研究生的第一個(gè)項(xiàng)目就是處理重測(cè)序數(shù)據(jù)，在做的過程中參考了很多資料，現(xiàn)在把整個(gè)重測(cè)序的上游分析流程分享給大家，一些關(guān)鍵環(huán)節(jié)已經(jīng)幫大家踩過坑了，相信跟著學(xué)習(xí)實(shí)踐一定能夠收獲很多。制作不易，希望大家多多關(guān)注點(diǎn)贊支持??！
0.準(zhǔn)備階段
在開始之前，我們需要做一些準(zhǔn)備工作，主要是部署好相關(guān)的軟件和工具。我們?cè)谶@個(gè)重測(cè)序數(shù)據(jù)分析過程中用到的所有軟件都是開源的，它們的代碼全部都能夠在github上找到，具體如下：
BWA: 這是最權(quán)威，使用最廣的NGS數(shù)據(jù)比對(duì)軟件；conda install bowtie2
Samtools: 是一個(gè)專門用于處理比對(duì)數(shù)據(jù)的工具；conda install -c bioconda bwa #-c 即channel，下載生信軟件都要加上-c bioconda software
Picard: 它是目前最著名的組學(xué)研究中心-Broad研究所開發(fā)的一款強(qiáng)大的NGS數(shù)據(jù)處理工具，功能方面和Samtools有些重疊，但更多的是互補(bǔ)，它是由java編寫的，我們直接下載最新的.jar包就行了。conda install -c bioconda picard
GATK :目前GATK有3.x和4.x兩個(gè)不同的版本，代碼在github上也是分開的。4.x是今年新推出的，在核心算法層面并沒太多的修改，但使用了新的設(shè)計(jì)模式，做了很多功能的整合，是更適合于大規(guī)模集群和云運(yùn)算的版本，后續(xù)GATK團(tuán)隊(duì)也將主要維護(hù)4.x的版本，而且它的代碼是100%開源的，這和3.x只有部分開源的情況不同?？吹贸鯣ATK今年的這次升級(jí)是為了應(yīng)對(duì)接下來(lái)越來(lái)越多的大規(guī)模人群測(cè)序數(shù)據(jù)而做出的改變，但現(xiàn)階段4.x版本還不穩(wěn)定，真正在使用的人和機(jī)構(gòu)其實(shí)也還不多。短期來(lái)看，3.x版本還將在業(yè)內(nèi)繼續(xù)使用一段時(shí)間；其次，3.x對(duì)于絕大分部的分析需求來(lái)說是完全足夠的。我們?cè)谶@里以GATK4作為流程的重要工具進(jìn)行分析流程的構(gòu)建。conda install -c bioconda gatk4
事實(shí)上，對(duì)于構(gòu)造重測(cè)序分析流程來(lái)說，以上這個(gè)四個(gè)工具就完全足夠了。它們的安裝都非常簡(jiǎn)單，除了BWA和Samtools由C編寫的，安裝時(shí)需要進(jìn)行編譯之外，另外兩個(gè)只要保證系統(tǒng)中的java是1.8.x版本及以上的，那么直接下載jar包就可以使用了。操作系統(tǒng)方面推薦linux（集群）或者M(jìn)ac OS。
1.原始數(shù)據(jù)質(zhì)控
這一部分考慮到現(xiàn)在公司給出的測(cè)序數(shù)據(jù)已經(jīng)進(jìn)行了質(zhì)控，到我們手里基本都是clean data,網(wǎng)上關(guān)于質(zhì)控的資源也比較多，基本上通過fastqc和fastp就可以進(jìn)行質(zhì)控和數(shù)據(jù)清洗。所以這里不再贅述。
2.數(shù)據(jù)預(yù)處理
序列比對(duì)
先問一個(gè)問題：為什么需要比對(duì)？
我們已經(jīng)知道NGS測(cè)序下來(lái)的短序列（read）存儲(chǔ)于FASTQ文件里面。雖然它們?cè)径紒?lái)自于有序的基因組，但在經(jīng)過DNA建庫(kù)和測(cè)序之后，文件中不同read之間的前后順序關(guān)系就已經(jīng)全部丟失了。因此，F(xiàn)ASTQ文件中緊挨著的兩條read之間沒有任何位置關(guān)系，它們都是隨機(jī)來(lái)自于原本基因組中某個(gè)位置的短序列而已。
因此，我們需要先把這一大堆的短序列捋順，一個(gè)個(gè)去跟該物種的 參考基因組【注】比較，找到每一條read在參考基因組上的位置，然后按順序排列好，這個(gè)過程就稱為測(cè)序數(shù)據(jù)的比對(duì)。這 也是核心流程真正意義上的第一步，只有完成了這個(gè)序列比對(duì)我們才有下一步的數(shù)據(jù)分析。
【注】參考基因組：指該物種的基因組序列，是已經(jīng)組裝成的完整基因組序列，常作為該物種的標(biāo)準(zhǔn)參照物，比如人類基因組參考序列，fasta格式。
序列比對(duì)本質(zhì)上是一個(gè)尋找最大公共子字符串的過程。大家如果有學(xué)過生物信息學(xué)的話，應(yīng)該或多或少知道BLAST，它使用的是動(dòng)態(tài)規(guī)劃的算法來(lái)尋找這樣的子串，但在面對(duì)巨量的短序列數(shù)據(jù)時(shí)，類似BLAST這樣的軟件實(shí)在太慢了！因此，需要更加有效的數(shù)據(jù)結(jié)構(gòu)和相應(yīng)的算法來(lái)完成這個(gè)搜索定位的任務(wù)。
我們這里將用于流程構(gòu)建的BWA就是其中最優(yōu)秀的一個(gè)，它將BW(Burrows-Wheeler)壓縮算法和后綴樹相結(jié)合，能夠讓我們以較小的時(shí)間和空間代價(jià)，獲得準(zhǔn)確的序列比對(duì)結(jié)果。
以下我們就開始流程的搭建。
首先，我們需要為參考基因組的構(gòu)建索引——這其實(shí)是在為參考序列進(jìn)行Burrows Wheeler變換（wiki: 塊排序壓縮），以便能夠在序列比對(duì)的時(shí)候進(jìn)行快速的搜索和定位。
$ bwa index genome.fasta #對(duì)基因組文件進(jìn)行bwa索引
完成之后，你會(huì)看到類似如下幾個(gè)以genome.fasta為前綴的文件：
這些就是在比對(duì)時(shí)真正需要被用到的文件。這一步完成之后，我們就可以將read比對(duì)至參考基因組了：
.├── genome.fasta.amb
├── genome.fasta.ann
├── genome.fasta.bwt
├── genome.fasta.pac
└── genome.fasta.sa
這些就是在比對(duì)時(shí)真正需要被用到的文件。這一步完成之后，我們就可以將read比對(duì)至參考基因組了：
$ bwa mem -t 4 -R '@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name' /path/to/genome.fasta read_1.fq.gz read_2.fq.gz | samtools view -S -b > sample.bam
大伙如果以前沒使用過這個(gè)比對(duì)工具的話，那么可能不明白上面參數(shù)的含義。我們這里調(diào)用的是bwa的mem比對(duì)模塊，在解釋這樣做之前，我們不妨先看一下bwa mem的官方用法說明，它就一句話：
Usage: bwa mem [options] <idxbase> <in1.fq> [in2.fq]
其中，[options]是一系列可選的參數(shù)，暫時(shí)不多說。這里的 < idxbase>要輸入的是參考基因組的BW索引文件，我們上面通過bwa index構(gòu)建好的那幾個(gè)以human.fasta為前綴的文件便是；< in1.fq>和 [in2.fq]輸入的是質(zhì)控后的fastq文件。但這里輸入的時(shí)候?yàn)槭裁磿?huì)需要兩個(gè)fq（in1.fq和in2.fq）呢？我們上面的例子也是有兩個(gè)：read_1.fq.gz和read_2.fq.gz。這是因?yàn)檫@是雙末端測(cè)序（也稱Pair-End）的情況，那什么是“雙末端測(cè)序”呢？這兩個(gè)fq之間的關(guān)系又是什么？這個(gè)我需要簡(jiǎn)單解析一下。
我們已經(jīng)知道NGS是短讀長(zhǎng)的測(cè)序技術(shù)，一次測(cè)出來(lái)的read的長(zhǎng)度都不會(huì)太長(zhǎng)，那為了盡可能把一個(gè)DNA序列片段盡可能多地測(cè)出來(lái)，既然測(cè)一邊不夠，那就測(cè)兩邊，于是就有了一種從被測(cè)DNA序列兩端各測(cè)序一次的模式，這就被稱為雙末端測(cè)序（Pair-End Sequencing，簡(jiǎn)稱PE測(cè)序）。如下圖是Pair-End測(cè)序的示意圖，中間灰色的是被測(cè)序的DNA序列片段，左邊黃色帶箭頭和右邊藍(lán)色帶箭頭的分別是測(cè)序出來(lái)的read1和read2序列，這里假定它們的長(zhǎng)度都是100bp。雖然很多時(shí)候Pair-End測(cè)序還是無(wú)法將整個(gè)被測(cè)的DNA片段完全測(cè)通，但是它依然提供了極其有用的信息，比如，我們知道每一對(duì)的read1和read2都來(lái)自于同一個(gè)DNA片段，read1和read2之間的距離是這個(gè)DNA片段的長(zhǎng)度，而且read1和read2的方向剛好是相反的（這里排除mate-pair的情形）等，這些信息對(duì)于后面的變異檢測(cè)等分析來(lái)說都是非常有用的。
Pair-End 測(cè)序
另外，在read1在fq1文件中位置和read2在fq2文件中的文件中的位置是相同的，而且read ID之間只在末尾有一個(gè)’/1’或者’/2’的差別。
read1 ID和read2 ID的差別
既然有雙末端測(cè)序，那么與之對(duì)應(yīng)的就有單末端測(cè)序（Single End Sequecing，簡(jiǎn)稱SE測(cè)序），即只測(cè)序其中一端。因此，我們?cè)谑褂胋wa比對(duì)的時(shí)候，實(shí)際上，in2.fq是非強(qiáng)制性的（所以用方括號(hào)括起來(lái)），只有是雙末端測(cè)序的數(shù)據(jù)時(shí)才需要添加。
回到上面我們的例子，大伙可以看到我這里除了用法中提到的參數(shù)之外，還多了2個(gè)額外的參數(shù)，分別是：-t，線程數(shù)，我們?cè)谶@里使用4個(gè)線程；-R 接的是 Read Group的字符串信息，這是一個(gè)非常重要的信息，以@RG開頭，它是用來(lái)將比對(duì)的read進(jìn)行分組的。不同的組之間測(cè)序過程被認(rèn)為是相互獨(dú)立的，這個(gè)信息對(duì)于我們后續(xù)對(duì)比對(duì)數(shù)據(jù)進(jìn)行錯(cuò)誤率分析和Mark duplicate時(shí)非常重要。在Read Group中，有如下幾個(gè)信息非常重要：
(1) ID，這是Read Group的分組ID，一般設(shè)置為測(cè)序的lane ID（不同lane之間的測(cè)序過程認(rèn)為是獨(dú)立的），下機(jī)數(shù)據(jù)中我們都能看到這個(gè)信息的，一般都是包含在fastq的文件名中；
(2) PL，指的是所用的測(cè)序平臺(tái)，這個(gè)信息不要隨便寫！特別是當(dāng)我們需要使用GATK進(jìn)行后續(xù)分析的時(shí)候，更是如此！這是一個(gè)很多新手都容易忽視的一個(gè)地方，在GATK中，PL只允許被設(shè)置為：ILLUMINA，SLX，SOLEXA，SOLID，454，LS454，COMPLETE，PACBIO，IONTORRENT，CAPILLARY，HELICOS或UNKNOWN這幾個(gè)信息。基本上就是目前市場(chǎng)上存在著的測(cè)序平臺(tái)，當(dāng)然，如果實(shí)在不知道，那么必須設(shè)置為UNKNOWN，名字方面不區(qū)分大小寫。如果你在分析的時(shí)候這里沒設(shè)置正確，那么在后續(xù)使用GATK過程中可能會(huì)碰到類似如下的錯(cuò)誤：
ERROR MESSAGE: The platform (xx) associated with read group GATKSAMReadGroupRecord @RG:xx is not a recognized platform.
這個(gè)時(shí)候你需要對(duì)比對(duì)文件的header信息進(jìn)行重寫，就會(huì)稍微比較麻煩。
我們上面的例子用的是PL:illumina。如果你的數(shù)據(jù)是CG測(cè)序的那么記得不要寫成CG！而要寫COMPLETE。
(3) SM，樣本ID，同樣非常重要，有時(shí)候我們測(cè)序的數(shù)據(jù)比較多的時(shí)候，那么可能會(huì)分成多個(gè)不同的lane分布測(cè)出來(lái)，這個(gè)時(shí)候SM名字就是可以用于區(qū)分這些樣本。
(4) LB，測(cè)序文庫(kù)的名字，這個(gè)重要性稍微低一些，主要也是為了協(xié)助區(qū)分不同的group而存在。文庫(kù)名字一般可以在下機(jī)的fq文件名中找到，如果上面的lane ID足夠用于區(qū)分的話，也可以不用設(shè)置LB；
除了以上這四個(gè)之外，還可以自定義添加其他的信息，不過如無(wú)特殊的需要，對(duì)于序列比對(duì)而言，這4個(gè)就足夠了。這些信息設(shè)置好之后，在RG字符串中要用制表符（\t）將它們分開。
最后在我們的例子中，我們將比對(duì)的輸出結(jié)果直接重定向到一份sample_name.sam文件中，這類文件是BWA比對(duì)的標(biāo)準(zhǔn)輸出文件，它的具體格式我會(huì)在下一篇文章中進(jìn)行詳細(xì)說明。但SAM文件是文本文件，一般整個(gè)文件都非常巨大，因此，為了有效節(jié)省磁盤空間，一般都會(huì)用samtools將它轉(zhuǎn)化為BAM文件（SAM的特殊二進(jìn)制格式），而且BAM會(huì)更加方便于后續(xù)的分析。所以我們上面比對(duì)的命令可以和samtools結(jié)合并改進(jìn)為：
$ bwa mem -t 4 -R '@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name' /path/to/human.fasta read_1.fq.gz read_2.fq.gz | samtools view -S -b - > sample_name.bam
我們通過管道(“|”)把比對(duì)的輸出如同引導(dǎo)水流一樣導(dǎo)流給samtools去處理，上面samtools view的-b參數(shù)指的就是輸出為BAM文件，這里需要注意的地方是-b后面的’-‘，它代表就是上面管道引流過來(lái)的數(shù)據(jù)，經(jīng)過samtools轉(zhuǎn)換之后我們?cè)僦囟ㄏ驗(yàn)閟ample_name.bam。
關(guān)于BWA的其他參數(shù)，我這里不打算對(duì)其進(jìn)行一一解釋，在絕大多數(shù)情況下，采用默認(rèn)是合適的做法。
[Tips] BWA MEM比對(duì)模塊是有一定適用范圍的：它是專門為長(zhǎng)read比對(duì)設(shè)計(jì)的，目的是為了解決，第三代測(cè)序技術(shù)這種能夠產(chǎn)生長(zhǎng)達(dá)幾十kb甚至幾Mbp的read情況。一般只有當(dāng)read長(zhǎng)度≥70bp的時(shí)候，才推薦使用，如果比這個(gè)要小，建議使用BWA ALN模塊。
對(duì)bam文件進(jìn)行進(jìn)行重新排序（reorder）
由BWA生成的sam文件時(shí)按字典式排序法進(jìn)行的排序（lexicographically）進(jìn)行排序的（chr10，chr11…chr19，chr1，chr20…chr22，chr2，chr3…chrM，chrX，chrY），但是GATK在進(jìn)行callsnp的時(shí)候是按照染色體組型（karyotypic）進(jìn)行的（chrM，chr1，chr2…chr22，chrX，chrY），因此要對(duì)原始sam文件進(jìn)行reorder。可以使用picard-tools中的ReorderSam完成。
Picard的SortSam需指定一個(gè)tmp目錄，用于存放中間文件，中間文件會(huì)很大，above 10G.注意指定目錄的空間大小。
排序
以上，我們就完成了read比對(duì)的步驟。接下來(lái)是排序：
排序這一步我們也是通過使用samtools來(lái)完成的，命令很簡(jiǎn)單：
Usage: samtools sort [options...] [in.bam]
但在執(zhí)行之前，我們有必要先搞明白為什么需要排序，為什么BWA比對(duì)后輸出的BAM文件是沒順序的！原因就是FASTQ文件里面這些被測(cè)序下來(lái)的read是隨機(jī)分布于基因組上面的，第一步的比對(duì)是按照FASTQ文件的順序把read逐一定位到參考基因組上之后，隨即就輸出了，它不會(huì)也不可能在這一步里面能夠自動(dòng)識(shí)別比對(duì)位置的先后位置重排比對(duì)結(jié)果。因此，比對(duì)后得到的結(jié)果文件中，每一條記錄之間位置的先后順序是亂的，我們后續(xù)去重復(fù)等步驟都需要在比對(duì)記錄按照順序從小到大排序下來(lái)才能進(jìn)行，所以這才是需要進(jìn)行排序的原因。對(duì)于我們的例子來(lái)說，這個(gè)排序的命令如下：
$ time samtools sort -@ 4 -m 4G -O bam -o sample_name.sorted.bam sample_name.bam
其中，-@，用于設(shè)定排序時(shí)的線程數(shù)，我們?cè)O(shè)為4；-m，限制排序時(shí)最大的內(nèi)存消耗，這里設(shè)為4GB；-O 指定輸出為bam格式；-o 是輸出文件的名字，這里叫sample_name.sorted.bam。我會(huì)比較建議大伙在做類似分析的時(shí)候在文件名字將所做的關(guān)鍵操作包含進(jìn)去，因?yàn)檫@樣即使過了很長(zhǎng)時(shí)間，當(dāng)你再去看這個(gè)文件的時(shí)候也能夠立刻知道當(dāng)時(shí)對(duì)它做了什么；最后就是輸入文件——sample_name.bam。
【注意】排序后如果發(fā)現(xiàn)新的BAM文件比原來(lái)的BAM文件稍微小一些，不用覺得驚訝，這是壓縮算法導(dǎo)致的結(jié)果，文件內(nèi)容是沒有損失的。
去除重復(fù)序列（或者標(biāo)記重復(fù)序列）
在排序完成之后我們就可以開始執(zhí)行去除重復(fù)（準(zhǔn)確來(lái)說是 去除PCR重復(fù)序列）的步驟了。
首先，我們需要先理解什么是重復(fù)序列，它是如何產(chǎn)生的，以及為什么需要去除掉？要回答這幾個(gè)問題，我們需要再次理解在建庫(kù)和測(cè)序時(shí)到底發(fā)生了什么。
我們知道在NGS測(cè)序之前都需要先構(gòu)建測(cè)序文庫(kù)：通過物理（超聲）打斷或者化學(xué)試劑（酶切）切斷原始的DNA序列，然后選擇特定長(zhǎng)度范圍的序列去進(jìn)行PCR擴(kuò)增并上機(jī)測(cè)序。
因此，這里重復(fù)序列的來(lái)源實(shí)際上就是由PCR過程中所引入的。因?yàn)樗^的PCR擴(kuò)增就是把原來(lái)的一段DNA序列復(fù)制多次。可是為什么需要PCR擴(kuò)增呢？如果沒有擴(kuò)增不就可以省略這一步了嗎？
情況確實(shí)如此，但是很多時(shí)候我們構(gòu)建測(cè)序文庫(kù)時(shí)能用的細(xì)胞量并不會(huì)非常充足，而且在打斷的步驟中也會(huì)引起部分DNA的降解，這兩點(diǎn)會(huì)使整體或者局部的DNA濃度過低，這時(shí)如果直接從這個(gè)溶液中取樣去測(cè)序就很可能漏掉原本基因組上的一些DNA片段，導(dǎo)致測(cè)序不全。而PCR擴(kuò)增的作用就是為了把這些微弱的DNA多復(fù)制幾倍乃至幾十倍，以便增大它們?cè)谌芤褐蟹植嫉拿芏?，使得能夠在取樣時(shí)被獲取到。所以這里大家需要記住一個(gè)重點(diǎn)，PCR擴(kuò)增原本的目的是為了增大微弱DNA序列片段的密度，但由于整個(gè)反應(yīng)都在一個(gè)試管中進(jìn)行，因此其他一些密度并不低的DNA片段也會(huì)被同步放大，那么這時(shí)在取樣去上機(jī)測(cè)序的時(shí)候，這些DNA片段就很可能會(huì)被重復(fù)取到相同的幾條去進(jìn)行測(cè)序（下圖為PCR擴(kuò)增示意圖）。
PCR擴(kuò)增示意圖：PCR擴(kuò)增是一個(gè)指數(shù)擴(kuò)增的過程，圖中原本只有一段雙鏈DNA序列，在經(jīng)過3輪PCR后就被擴(kuò)增成了8段
看到這里，你或許會(huì)覺得，那沒必要去除不也應(yīng)該可以嗎？因?yàn)榧幢銛U(kuò)增了多幾次，不也同樣還是原來(lái)的那一段DNA嗎？直接用來(lái)分析對(duì)結(jié)果也不會(huì)有影響??！難道不是嗎？
會(huì)有影響，而且有時(shí)影響會(huì)很大！最直接的后果就是同時(shí)增大了變異檢測(cè)結(jié)果的假陰和假陽(yáng)率。主要有幾個(gè)原因：
- DNA在打斷的那一步會(huì)發(fā)生一些損失，主要表現(xiàn)是會(huì)引發(fā)一些堿基發(fā)生顛換變換（嘌呤-變嘧啶或者嘧啶變嘌呤），帶來(lái)假的變異。PCR過程會(huì)擴(kuò)大這個(gè)信號(hào)，導(dǎo)致最后的檢測(cè)結(jié)果中混入了假的結(jié)果；
- PCR反應(yīng)過程中也會(huì)帶來(lái)新的堿基錯(cuò)誤。發(fā)生在前幾輪的PCR擴(kuò)增發(fā)生的錯(cuò)誤會(huì)在后續(xù)的PCR過程中擴(kuò)大，同樣帶來(lái)假的變異；
- 對(duì)于真實(shí)的變異，PCR反應(yīng)可能會(huì)對(duì)包含某一個(gè)堿基的DNA模版擴(kuò)增更加劇烈（這個(gè)現(xiàn)象稱為PCR Bias）。如果反應(yīng)體系是對(duì)含有reference allele的模板擴(kuò)增偏向強(qiáng)烈，那么變異堿基的信息會(huì)變小，從而會(huì)導(dǎo)致假陰。
PCR對(duì)真實(shí)的變異檢測(cè)和個(gè)體的基因型判斷都有不好的影響。GATK、Samtools、Platpus等這種利用貝葉斯原理的變異檢測(cè)算法都是認(rèn)為所用的序列數(shù)據(jù)都不是重復(fù)序列（即將它們和其他序列一視同仁地進(jìn)行變異的判斷，所以帶來(lái)誤導(dǎo)），因此必須要進(jìn)行標(biāo)記（去除）或者使用PCR-Free的測(cè)序方案（這個(gè)方案目前正變得越來(lái)越流行，特別是對(duì)于RNA-Seq來(lái)說尤為重要，現(xiàn)在著名的基因組學(xué)研究所——Broad Institute，基本都是使用PCR-Free的測(cè)序方案）。
那么具體是如何做到去除這些PCR重復(fù)序列的呢？我們可以拋開任何工具，仔細(xì)想想，既然PCR擴(kuò)增是把同一段DNA序列復(fù)制出很多份，那么這些序列在經(jīng)過比對(duì)之后它們一定會(huì)定位到基因組上相同的位置，比對(duì)的信息看起來(lái)也將是一樣的！于是，我們就可以根據(jù)這個(gè)特點(diǎn)找到這些重復(fù)序列了！
事實(shí)上，現(xiàn)有的工具包括Samtools和Picard中去除重復(fù)序列的算法也的確是這么做的。不同的地方在于，samtools的rmdup是直接將這些重復(fù)序列從比對(duì)BAM文件中刪除掉，而Picard的MarkDuplicates默認(rèn)情況則只是在BAM的FLAG信息中標(biāo)記出來(lái)，而不是刪除，因此這些重復(fù)序列依然會(huì)被留在文件中，只是我們可以在變異檢測(cè)的時(shí)候識(shí)別到它們，并進(jìn)行忽略。
考慮到盡可能和現(xiàn)在主流的做法一致（但我并不是說主流的做法就一定是對(duì)的，要分情況看待，只是主流的做法容易被做成生產(chǎn)流程而已），我們這里也用Picard來(lái)完成這個(gè)事情：
picard MarkDuplicates I=sample_name.sorted.bam O=sample_name.sorted.markdup.bam M=sample_name.markdup_metrics.txt
這里只把重復(fù)序列在輸出的新結(jié)果中標(biāo)記出來(lái)，但不刪除。如果我們非要把這些序列完全刪除的話可以這樣做：java -jar 把參數(shù)REMOVE_DUPLICATES設(shè)置為ture，那么重復(fù)序列就被刪除掉，不會(huì)在結(jié)果文件中留存。我比較建議使用第一種做法，只是標(biāo)記出來(lái)，并留存這些序列，以便在你需要的時(shí)候還可以對(duì)其做分析。
picard MarkDuplicates REMOVE_DUPLICATES=true I=sample_name.sorted.bam O=sample_name.sorted.markdup.bam M=sample_name.markdup_metrics.txt
創(chuàng)建索引
這一步完成之后，我們需要為sample_name.sorted.markdup.bam創(chuàng)建索引文件，它的作用能夠讓我們可以隨機(jī)訪問這個(gè)文件中的任意位置，而且后面的“局部重比對(duì)”步驟也要求這個(gè)BAM文件一定要有索引，命令如下：
samtools index sample_name.sorted.markdup.bam
完成之后，會(huì)生成一份sample_name.sorted.markdup.bam.bai文件，這就是上面這份BAM的index。
同時(shí)，還要對(duì)參考基因組進(jìn)行GATK的索引,也就是準(zhǔn)備參考基因組.fai和.dict文件
gatk CreateSequenceDictionary -R genome.fa -O genome.dict && samtools faidx genome.fa
生成gvcf文件并合并
有多個(gè)個(gè)體的情況下，首先我們用HaplotypeCaller命令給每一個(gè)個(gè)體生成gvcf文件，然后用CombineGVCFs按染色體合并gvcf文件。
我是用shell生成的批量腳本然后放到服務(wù)器上一起跑的，就會(huì)比較快。批量寫腳本什么語(yǔ)言都可以，就按照你的需求循環(huán)改變名稱、個(gè)體號(hào)、染色體號(hào)等各種字符就行
#1 生成gvcf文件
gatk HaplotypeCaller -R ref.fa -I sample_name.sorted.markdup.bam --emit-ref-confidence GVCF --min-base-quality-score 10 -O chrx.g.vcf.gz
#2 gvcf文件按染色體合并
ls chrx.g.vcf.gz > chrx_gvcf.list
gatk CombineGVCFs -R ref.fa -V chrx_gvcf.list -L X(染色體號(hào)) -O chrx.merged.g.vcf.gz
生成基因型文件(這一步變成vcf了）
gatk GenotypeGVCFs -R ref.fa -V chrx.merged.g.vcf.gz -O chrx.genotype.vcf.gz
過濾
上一步得到的是rawdata，我們還需要對(duì)原始數(shù)據(jù)做變異質(zhì)控。質(zhì)控是指通過一定的標(biāo)準(zhǔn)，最大可能地剔除假陽(yáng)性的結(jié)果，并盡可能地保留最多的正確數(shù)據(jù)。
SNP過濾有兩種情況，一種是僅根據(jù)位點(diǎn)質(zhì)量信息（測(cè)序深度，回帖質(zhì)量等）對(duì)SNP進(jìn)行粗過濾。另一種是考慮除了測(cè)序質(zhì)量以外的信息進(jìn)行的過濾。 接下來(lái)我會(huì)分別介紹到這兩種過濾。
從測(cè)序位點(diǎn)質(zhì)量上看，在GATK HaplotypeCaller之后，首選的質(zhì)控方案是GATK VQSR， 原理是利用自身數(shù)據(jù)和已知變異位點(diǎn)集的overlap，通過GMM模型構(gòu)建一個(gè)分類器來(lái)對(duì)變異數(shù)據(jù)進(jìn)行打分，從而評(píng)估每個(gè)位點(diǎn)的可行度。
雖然經(jīng)?？吹絍QSR的教程，但是很可惜目前應(yīng)該只有人類上可以做（還是需要高質(zhì)量的已知變異集），所以我們其他物種只能做hard-filtering硬過濾了。
挑選snp
gatk SelectVariants -R genome.fa -variant sample.vcf -O sample_snp.vcf -select-type SNP
過濾snp
gatk VariantFiltration -variant sample_snp.vcf -O sample_snp_filter.vcf -R genome.fa
vcf文件合并
請(qǐng)查看http://www.itdecent.cn/p/3d86879f6f5c
