自己在研究生的第一個項目就是處理重測序數(shù)據(jù)拭卿，在做的過程中參考了很多資料骡湖，現(xiàn)在把整個重測序的上游分析流程分享給大家，一些關(guān)鍵環(huán)節(jié)已經(jīng)幫大家踩過坑了峻厚，相信跟著學習實踐一定能夠收獲很多响蕴。制作不易，希望大家多多關(guān)注點贊支持；萏摇浦夷！
0.準備階段
在開始之前，我們需要做一些準備工作辜王，主要是部署好相關(guān)的軟件和工具劈狐。我們在這個重測序數(shù)據(jù)分析過程中用到的所有軟件都是開源的，它們的代碼全部都能夠在github上找到呐馆，具體如下：
BWA: 這是最權(quán)威肥缔，使用最廣的NGS數(shù)據(jù)比對軟件；conda install bowtie2
Samtools: 是一個專門用于處理比對數(shù)據(jù)的工具摹恰；conda install -c bioconda bwa #-c 即channel辫继，下載生信軟件都要加上-c bioconda software
Picard: 它是目前最著名的組學研究中心-Broad研究所開發(fā)的一款強大的NGS數(shù)據(jù)處理工具，功能方面和Samtools有些重疊俗慈，但更多的是互補姑宽，它是由java編寫的，我們直接下載最新的.jar包就行了闺阱。conda install -c bioconda picard
GATK :目前GATK有3.x和4.x兩個不同的版本炮车，代碼在github上也是分開的。4.x是今年新推出的酣溃，在核心算法層面并沒太多的修改瘦穆，但使用了新的設(shè)計模式，做了很多功能的整合赊豌，是更適合于大規(guī)模集群和云運算的版本扛或，后續(xù)GATK團隊也將主要維護4.x的版本，而且它的代碼是100%開源的碘饼，這和3.x只有部分開源的情況不同熙兔”妫看得出GATK今年的這次升級是為了應(yīng)對接下來越來越多的大規(guī)模人群測序數(shù)據(jù)而做出的改變，但現(xiàn)階段4.x版本還不穩(wěn)定住涉，真正在使用的人和機構(gòu)其實也還不多麸锉。短期來看，3.x版本還將在業(yè)內(nèi)繼續(xù)使用一段時間舆声；其次花沉，3.x對于絕大分部的分析需求來說是完全足夠的。我們在這里以GATK4作為流程的重要工具進行分析流程的構(gòu)建媳握。conda install -c bioconda gatk4
事實上碱屁，對于構(gòu)造重測序分析流程來說，以上這個四個工具就完全足夠了毙芜。它們的安裝都非常簡單忽媒，除了BWA和Samtools由C編寫的，安裝時需要進行編譯之外腋粥，另外兩個只要保證系統(tǒng)中的java是1.8.x版本及以上的，那么直接下載jar包就可以使用了架曹。操作系統(tǒng)方面推薦linux（集群）或者Mac OS隘冲。
1.原始數(shù)據(jù)質(zhì)控
這一部分考慮到現(xiàn)在公司給出的測序數(shù)據(jù)已經(jīng)進行了質(zhì)控，到我們手里基本都是clean data,網(wǎng)上關(guān)于質(zhì)控的資源也比較多绑雄，基本上通過fastqc和fastp就可以進行質(zhì)控和數(shù)據(jù)清洗展辞。所以這里不再贅述。
2.數(shù)據(jù)預處理
序列比對
先問一個問題：為什么需要比對万牺？
我們已經(jīng)知道NGS測序下來的短序列（read）存儲于FASTQ文件里面罗珍。雖然它們原本都來自于有序的基因組，但在經(jīng)過DNA建庫和測序之后脚粟，文件中不同read之間的前后順序關(guān)系就已經(jīng)全部丟失了覆旱。因此，F(xiàn)ASTQ文件中緊挨著的兩條read之間沒有任何位置關(guān)系核无，它們都是隨機來自于原本基因組中某個位置的短序列而已扣唱。
因此，我們需要先把這一大堆的短序列捋順团南，一個個去跟該物種的 參考基因組【注】比較噪沙，找到每一條read在參考基因組上的位置，然后按順序排列好吐根，這個過程就稱為測序數(shù)據(jù)的比對正歼。這 也是核心流程真正意義上的第一步，只有完成了這個序列比對我們才有下一步的數(shù)據(jù)分析拷橘。
【注】參考基因組：指該物種的基因組序列局义，是已經(jīng)組裝成的完整基因組序列齐疙，常作為該物種的標準參照物，比如人類基因組參考序列旭咽，fasta格式贞奋。
序列比對本質(zhì)上是一個尋找最大公共子字符串的過程。大家如果有學過生物信息學的話穷绵，應(yīng)該或多或少知道BLAST轿塔，它使用的是動態(tài)規(guī)劃的算法來尋找這樣的子串，但在面對巨量的短序列數(shù)據(jù)時仲墨，類似BLAST這樣的軟件實在太慢了勾缭！因此，需要更加有效的數(shù)據(jù)結(jié)構(gòu)和相應(yīng)的算法來完成這個搜索定位的任務(wù)目养。
我們這里將用于流程構(gòu)建的BWA就是其中最優(yōu)秀的一個俩由，它將BW(Burrows-Wheeler)壓縮算法和后綴樹相結(jié)合，能夠讓我們以較小的時間和空間代價癌蚁，獲得準確的序列比對結(jié)果幻梯。
以下我們就開始流程的搭建。
首先努释，我們需要為參考基因組的構(gòu)建索引——這其實是在為參考序列進行Burrows Wheeler變換（wiki: 塊排序壓縮）碘梢，以便能夠在序列比對的時候進行快速的搜索和定位。
$ bwa index genome.fasta #對基因組文件進行bwa索引
完成之后伐蒂，你會看到類似如下幾個以genome.fasta為前綴的文件：
這些就是在比對時真正需要被用到的文件煞躬。這一步完成之后，我們就可以將read比對至參考基因組了：
.├── genome.fasta.amb
├── genome.fasta.ann
├── genome.fasta.bwt
├── genome.fasta.pac
└── genome.fasta.sa
這些就是在比對時真正需要被用到的文件逸邦。這一步完成之后恩沛，我們就可以將read比對至參考基因組了：
$ 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
大伙如果以前沒使用過這個比對工具的話，那么可能不明白上面參數(shù)的含義缕减。我們這里調(diào)用的是bwa的mem比對模塊雷客，在解釋這樣做之前，我們不妨先看一下bwa mem的官方用法說明烛卧，它就一句話：
Usage: bwa mem [options] <idxbase> <in1.fq> [in2.fq]
其中佛纫，[options]是一系列可選的參數(shù)，暫時不多說总放。這里的 < idxbase>要輸入的是參考基因組的BW索引文件呈宇，我們上面通過bwa index構(gòu)建好的那幾個以human.fasta為前綴的文件便是；< in1.fq>和 [in2.fq]輸入的是質(zhì)控后的fastq文件局雄。但這里輸入的時候為什么會需要兩個fq（in1.fq和in2.fq）呢甥啄？我們上面的例子也是有兩個：read_1.fq.gz和read_2.fq.gz。這是因為這是雙末端測序（也稱Pair-End）的情況炬搭，那什么是“雙末端測序”呢蜈漓？這兩個fq之間的關(guān)系又是什么穆桂？這個我需要簡單解析一下。
我們已經(jīng)知道NGS是短讀長的測序技術(shù)融虽，一次測出來的read的長度都不會太長享完，那為了盡可能把一個DNA序列片段盡可能多地測出來，既然測一邊不夠有额，那就測兩邊般又，于是就有了一種從被測DNA序列兩端各測序一次的模式，這就被稱為雙末端測序（Pair-End Sequencing巍佑，簡稱PE測序）茴迁。如下圖是Pair-End測序的示意圖，中間灰色的是被測序的DNA序列片段，左邊黃色帶箭頭和右邊藍色帶箭頭的分別是測序出來的read1和read2序列，這里假定它們的長度都是100bp扬跋。雖然很多時候Pair-End測序還是無法將整個被測的DNA片段完全測通，但是它依然提供了極其有用的信息倦卖，比如，我們知道每一對的read1和read2都來自于同一個DNA片段筹吐，read1和read2之間的距離是這個DNA片段的長度糖耸，而且read1和read2的方向剛好是相反的（這里排除mate-pair的情形）等，這些信息對于后面的變異檢測等分析來說都是非常有用的丘薛。
Pair-End 測序
另外，在read1在fq1文件中位置和read2在fq2文件中的文件中的位置是相同的邦危，而且read ID之間只在末尾有一個’/1’或者’/2’的差別洋侨。
read1 ID和read2 ID的差別
既然有雙末端測序，那么與之對應(yīng)的就有單末端測序（Single End Sequecing倦蚪，簡稱SE測序）希坚，即只測序其中一端。因此陵且，我們在使用bwa比對的時候裁僧，實際上，in2.fq是非強制性的（所以用方括號括起來）慕购，只有是雙末端測序的數(shù)據(jù)時才需要添加聊疲。
回到上面我們的例子，大伙可以看到我這里除了用法中提到的參數(shù)之外沪悲，還多了2個額外的參數(shù)获洲，分別是：-t，線程數(shù)殿如，我們在這里使用4個線程贡珊；-R 接的是 Read Group的字符串信息最爬，這是一個非常重要的信息，以@RG開頭门岔，它是用來將比對的read進行分組的爱致。不同的組之間測序過程被認為是相互獨立的，這個信息對于我們后續(xù)對比對數(shù)據(jù)進行錯誤率分析和Mark duplicate時非常重要寒随。在Read Group中糠悯，有如下幾個信息非常重要：
(1) ID，這是Read Group的分組ID牢裳，一般設(shè)置為測序的lane ID（不同lane之間的測序過程認為是獨立的）逢防，下機數(shù)據(jù)中我們都能看到這個信息的，一般都是包含在fastq的文件名中蒲讯；
(2) PL忘朝，指的是所用的測序平臺，這個信息不要隨便寫判帮！特別是當我們需要使用GATK進行后續(xù)分析的時候局嘁，更是如此！這是一個很多新手都容易忽視的一個地方晦墙，在GATK中悦昵，PL只允許被設(shè)置為：ILLUMINA，SLX晌畅，SOLEXA但指，SOLID，454抗楔，LS454棋凳，COMPLETE，PACBIO连躏，IONTORRENT剩岳，CAPILLARY，HELICOS或UNKNOWN這幾個信息入热∨淖兀基本上就是目前市場上存在著的測序平臺，當然勺良，如果實在不知道绰播，那么必須設(shè)置為UNKNOWN，名字方面不區(qū)分大小寫郑气。如果你在分析的時候這里沒設(shè)置正確幅垮，那么在后續(xù)使用GATK過程中可能會碰到類似如下的錯誤：
ERROR MESSAGE: The platform (xx) associated with read group GATKSAMReadGroupRecord @RG:xx is not a recognized platform.
這個時候你需要對比對文件的header信息進行重寫，就會稍微比較麻煩。
我們上面的例子用的是PL:illumina忙芒。如果你的數(shù)據(jù)是CG測序的那么記得不要寫成CG示弓！而要寫COMPLETE。
(3) SM呵萨，樣本ID奏属，同樣非常重要，有時候我們測序的數(shù)據(jù)比較多的時候潮峦，那么可能會分成多個不同的lane分布測出來囱皿，這個時候SM名字就是可以用于區(qū)分這些樣本。
(4) LB忱嘹，測序文庫的名字嘱腥，這個重要性稍微低一些，主要也是為了協(xié)助區(qū)分不同的group而存在拘悦。文庫名字一般可以在下機的fq文件名中找到齿兔，如果上面的lane ID足夠用于區(qū)分的話，也可以不用設(shè)置LB础米；
除了以上這四個之外分苇，還可以自定義添加其他的信息，不過如無特殊的需要屁桑，對于序列比對而言医寿，這4個就足夠了。這些信息設(shè)置好之后蘑斧，在RG字符串中要用制表符（\t）將它們分開靖秩。
最后在我們的例子中，我們將比對的輸出結(jié)果直接重定向到一份sample_name.sam文件中竖瘾，這類文件是BWA比對的標準輸出文件盆偿，它的具體格式我會在下一篇文章中進行詳細說明。但SAM文件是文本文件准浴，一般整個文件都非常巨大，因此捎稚，為了有效節(jié)省磁盤空間乐横，一般都會用samtools將它轉(zhuǎn)化為BAM文件（SAM的特殊二進制格式），而且BAM會更加方便于后續(xù)的分析今野。所以我們上面比對的命令可以和samtools結(jié)合并改進為：
$ 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
我們通過管道(“|”)把比對的輸出如同引導水流一樣導流給samtools去處理葡公，上面samtools view的-b參數(shù)指的就是輸出為BAM文件，這里需要注意的地方是-b后面的’-‘条霜，它代表就是上面管道引流過來的數(shù)據(jù)催什，經(jīng)過samtools轉(zhuǎn)換之后我們再重定向為sample_name.bam。
關(guān)于BWA的其他參數(shù)宰睡，我這里不打算對其進行一一解釋蒲凶，在絕大多數(shù)情況下气筋，采用默認是合適的做法。
[Tips] BWA MEM比對模塊是有一定適用范圍的：它是專門為長read比對設(shè)計的旋圆，目的是為了解決宠默，第三代測序技術(shù)這種能夠產(chǎn)生長達幾十kb甚至幾Mbp的read情況。一般只有當read長度≥70bp的時候灵巧，才推薦使用搀矫，如果比這個要小，建議使用BWA ALN模塊刻肄。
對bam文件進行進行重新排序（reorder）
由BWA生成的sam文件時按字典式排序法進行的排序（lexicographically）進行排序的（chr10瓤球，chr11…chr19，chr1敏弃，chr20…chr22卦羡，chr2，chr3…chrM权她，chrX虹茶，chrY），但是GATK在進行callsnp的時候是按照染色體組型（karyotypic）進行的（chrM隅要，chr1蝴罪，chr2…chr22，chrX步清，chrY）要门，因此要對原始sam文件進行reorder±。可以使用picard-tools中的ReorderSam完成欢搜。
Picard的SortSam需指定一個tmp目錄，用于存放中間文件谴轮，中間文件會很大炒瘟，above 10G.注意指定目錄的空間大小。
排序
以上第步，我們就完成了read比對的步驟疮装。接下來是排序：
排序這一步我們也是通過使用samtools來完成的，命令很簡單：
Usage: samtools sort [options...] [in.bam]
但在執(zhí)行之前粘都，我們有必要先搞明白為什么需要排序廓推，為什么BWA比對后輸出的BAM文件是沒順序的！原因就是FASTQ文件里面這些被測序下來的read是隨機分布于基因組上面的翩隧，第一步的比對是按照FASTQ文件的順序把read逐一定位到參考基因組上之后樊展，隨即就輸出了，它不會也不可能在這一步里面能夠自動識別比對位置的先后位置重排比對結(jié)果。因此专缠，比對后得到的結(jié)果文件中雷酪，每一條記錄之間位置的先后順序是亂的，我們后續(xù)去重復等步驟都需要在比對記錄按照順序從小到大排序下來才能進行藤肢，所以這才是需要進行排序的原因太闺。對于我們的例子來說，這個排序的命令如下：
$ time samtools sort -@ 4 -m 4G -O bam -o sample_name.sorted.bam sample_name.bam
其中嘁圈，-@省骂，用于設(shè)定排序時的線程數(shù)，我們設(shè)為4最住；-m钞澳，限制排序時最大的內(nèi)存消耗，這里設(shè)為4GB涨缚；-O 指定輸出為bam格式轧粟；-o 是輸出文件的名字，這里叫sample_name.sorted.bam脓魏。我會比較建議大伙在做類似分析的時候在文件名字將所做的關(guān)鍵操作包含進去兰吟，因為這樣即使過了很長時間，當你再去看這個文件的時候也能夠立刻知道當時對它做了什么茂翔；最后就是輸入文件——sample_name.bam混蔼。
【注意】排序后如果發(fā)現(xiàn)新的BAM文件比原來的BAM文件稍微小一些，不用覺得驚訝珊燎，這是壓縮算法導致的結(jié)果惭嚣，文件內(nèi)容是沒有損失的。
去除重復序列（或者標記重復序列）
在排序完成之后我們就可以開始執(zhí)行去除重復（準確來說是 去除PCR重復序列）的步驟了悔政。
首先晚吞，我們需要先理解什么是重復序列，它是如何產(chǎn)生的谋国，以及為什么需要去除掉槽地？要回答這幾個問題，我們需要再次理解在建庫和測序時到底發(fā)生了什么芦瘾。
我們知道在NGS測序之前都需要先構(gòu)建測序文庫：通過物理（超聲）打斷或者化學試劑（酶切）切斷原始的DNA序列闷盔，然后選擇特定長度范圍的序列去進行PCR擴增并上機測序。
因此旅急，這里重復序列的來源實際上就是由PCR過程中所引入的。因為所謂的PCR擴增就是把原來的一段DNA序列復制多次牡整。可是為什么需要PCR擴增呢藐吮？如果沒有擴增不就可以省略這一步了嗎？
情況確實如此，但是很多時候我們構(gòu)建測序文庫時能用的細胞量并不會非常充足谣辞，而且在打斷的步驟中也會引起部分DNA的降解迫摔，這兩點會使整體或者局部的DNA濃度過低，這時如果直接從這個溶液中取樣去測序就很可能漏掉原本基因組上的一些DNA片段泥从，導致測序不全句占。而PCR擴增的作用就是為了把這些微弱的DNA多復制幾倍乃至幾十倍，以便增大它們在溶液中分布的密度躯嫉，使得能夠在取樣時被獲取到纱烘。所以這里大家需要記住一個重點，PCR擴增原本的目的是為了增大微弱DNA序列片段的密度祈餐，但由于整個反應(yīng)都在一個試管中進行擂啥，因此其他一些密度并不低的DNA片段也會被同步放大，那么這時在取樣去上機測序的時候帆阳，這些DNA片段就很可能會被重復取到相同的幾條去進行測序（下圖為PCR擴增示意圖）哺壶。
PCR擴增示意圖：PCR擴增是一個指數(shù)擴增的過程，圖中原本只有一段雙鏈DNA序列蜒谤，在經(jīng)過3輪PCR后就被擴增成了8段
看到這里山宾，你或許會覺得，那沒必要去除不也應(yīng)該可以嗎鳍徽？因為即便擴增了多幾次资锰，不也同樣還是原來的那一段DNA嗎？直接用來分析對結(jié)果也不會有影響把ⅰ台妆！難道不是嗎？
會有影響胖翰，而且有時影響會很大接剩！最直接的后果就是同時增大了變異檢測結(jié)果的假陰和假陽率。主要有幾個原因：
- DNA在打斷的那一步會發(fā)生一些損失萨咳，主要表現(xiàn)是會引發(fā)一些堿基發(fā)生顛換變換（嘌呤-變嘧啶或者嘧啶變嘌呤）懊缺，帶來假的變異。PCR過程會擴大這個信號培他，導致最后的檢測結(jié)果中混入了假的結(jié)果鹃两；
- PCR反應(yīng)過程中也會帶來新的堿基錯誤。發(fā)生在前幾輪的PCR擴增發(fā)生的錯誤會在后續(xù)的PCR過程中擴大舀凛，同樣帶來假的變異俊扳；
- 對于真實的變異，PCR反應(yīng)可能會對包含某一個堿基的DNA模版擴增更加劇烈（這個現(xiàn)象稱為PCR Bias）猛遍。如果反應(yīng)體系是對含有reference allele的模板擴增偏向強烈馋记，那么變異堿基的信息會變小号坡，從而會導致假陰。
PCR對真實的變異檢測和個體的基因型判斷都有不好的影響梯醒。GATK宽堆、Samtools、Platpus等這種利用貝葉斯原理的變異檢測算法都是認為所用的序列數(shù)據(jù)都不是重復序列（即將它們和其他序列一視同仁地進行變異的判斷茸习，所以帶來誤導）畜隶，因此必須要進行標記（去除）或者使用PCR-Free的測序方案（這個方案目前正變得越來越流行，特別是對于RNA-Seq來說尤為重要号胚，現(xiàn)在著名的基因組學研究所——Broad Institute籽慢，基本都是使用PCR-Free的測序方案）。
那么具體是如何做到去除這些PCR重復序列的呢涕刚？我們可以拋開任何工具嗡综，仔細想想，既然PCR擴增是把同一段DNA序列復制出很多份杜漠，那么這些序列在經(jīng)過比對之后它們一定會定位到基因組上相同的位置极景，比對的信息看起來也將是一樣的！于是驾茴，我們就可以根據(jù)這個特點找到這些重復序列了盼樟！
事實上，現(xiàn)有的工具包括Samtools和Picard中去除重復序列的算法也的確是這么做的锈至。不同的地方在于晨缴，samtools的rmdup是直接將這些重復序列從比對BAM文件中刪除掉，而Picard的MarkDuplicates默認情況則只是在BAM的FLAG信息中標記出來峡捡，而不是刪除击碗，因此這些重復序列依然會被留在文件中，只是我們可以在變異檢測的時候識別到它們们拙，并進行忽略稍途。
考慮到盡可能和現(xiàn)在主流的做法一致（但我并不是說主流的做法就一定是對的，要分情況看待砚婆，只是主流的做法容易被做成生產(chǎn)流程而已）械拍，我們這里也用Picard來完成這個事情：
picard MarkDuplicates I=sample_name.sorted.bam O=sample_name.sorted.markdup.bam M=sample_name.markdup_metrics.txt
這里只把重復序列在輸出的新結(jié)果中標記出來，但不刪除装盯。如果我們非要把這些序列完全刪除的話可以這樣做：java -jar 把參數(shù)REMOVE_DUPLICATES設(shè)置為ture坷虑，那么重復序列就被刪除掉，不會在結(jié)果文件中留存埂奈。我比較建議使用第一種做法迄损，只是標記出來，并留存這些序列账磺，以便在你需要的時候還可以對其做分析海蔽。
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)建索引文件，它的作用能夠讓我們可以隨機訪問這個文件中的任意位置党窜，而且后面的“局部重比對”步驟也要求這個BAM文件一定要有索引，命令如下：
samtools index sample_name.sorted.markdup.bam
完成之后借宵，會生成一份sample_name.sorted.markdup.bam.bai文件幌衣，這就是上面這份BAM的index。
同時壤玫，還要對參考基因組進行GATK的索引,也就是準備參考基因組.fai和.dict文件
gatk CreateSequenceDictionary -R genome.fa -O genome.dict && samtools faidx genome.fa
生成gvcf文件并合并
有多個個體的情況下豁护，首先我們用HaplotypeCaller命令給每一個個體生成gvcf文件，然后用CombineGVCFs按染色體合并gvcf文件欲间。
我是用shell生成的批量腳本然后放到服務(wù)器上一起跑的楚里，就會比較快。批量寫腳本什么語言都可以猎贴，就按照你的需求循環(huán)改變名稱班缎、個體號、染色體號等各種字符就行
#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(染色體號) -O chrx.merged.g.vcf.gz
生成基因型文件(這一步變成vcf了）
gatk GenotypeGVCFs -R ref.fa -V chrx.merged.g.vcf.gz -O chrx.genotype.vcf.gz
過濾
上一步得到的是rawdata她渴，我們還需要對原始數(shù)據(jù)做變異質(zhì)控达址。質(zhì)控是指通過一定的標準，最大可能地剔除假陽性的結(jié)果趁耗，并盡可能地保留最多的正確數(shù)據(jù)沉唠。
SNP過濾有兩種情況，一種是僅根據(jù)位點質(zhì)量信息（測序深度苛败，回帖質(zhì)量等）對SNP進行粗過濾满葛。另一種是考慮除了測序質(zhì)量以外的信息進行的過濾。 接下來我會分別介紹到這兩種過濾罢屈。
從測序位點質(zhì)量上看嘀韧，在GATK HaplotypeCaller之后，首選的質(zhì)控方案是GATK VQSR儡遮， 原理是利用自身數(shù)據(jù)和已知變異位點集的overlap乳蛾，通過GMM模型構(gòu)建一個分類器來對變異數(shù)據(jù)進行打分，從而評估每個位點的可行度鄙币。
雖然經(jīng)乘嘁叮看到VQSR的教程，但是很可惜目前應(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文件合并
請查看http://www.reibang.com/p/3d86879f6f5c
