參考：堿基礦工
從零開始完整學(xué)習(xí)全基因組測序數(shù)據(jù)分析：第4節(jié) 構(gòu)建WGS主流程

GATK4.0和全基因組數(shù)據(jù)分析實踐（上）

前言

由于人類基因組數(shù)據(jù)過大（我的小小服務(wù)器也只有50G）禀综，而人的參考序列長度就有3Gb，且高深度測序結(jié)果往往動輒100G。因此作者也嘗試用E.coli K12 的數(shù)據(jù)作為替代，一來它的基因不是很復(fù)雜，數(shù)據(jù)小轧房，它的基因組長度只有4.6Mb，二來也能很好的跑一遍流程，熟悉過程七兜。
咱們開始吧～

準(zhǔn)備階段

首先需要在服務(wù)器上部署相關(guān)的軟件，包括以下四個軟件福扬。
1）BWA (Burrow-Wheeler Aligner)：非常權(quán)威的NGS數(shù)據(jù)比對軟件腕铸，https://github.com/lh3/bwa。
2）Samtools：一個專門用于處理比對數(shù)據(jù)的軟件铛碑，https://github.com/samtools/samtools狠裹。
3）Picard：也是一款功能強大的NGS數(shù)據(jù)處理工具，功能和samtools 很類似汽烦，但更多的是與其互補涛菠，http://broadinstitute.github.io/picard/。
4）GATK：是一款非常權(quán)威的基因數(shù)據(jù)變異檢測工具撇吞。目前有3.x和4.x兩個不同的版本俗冻。4.x在核心算法層面并沒太多的修改，采用了新的設(shè)計模式牍颈，后續(xù)GATK團隊也將主要維護4.x的版本迄薄。但目前還不如3.x 版本穩(wěn)定，且短期來看3.x版本還將在業(yè)內(nèi)繼續(xù)使用一段時間煮岁。這里作者也以GATK3.8（最新版本）作為流程的重要工具進行分析流程的構(gòu)建讥蔽。https://software.broadinstitute.org/gatk/download/

對于WGS 分析來說涣易，這四個工具就夠了。

安裝

關(guān)于conda 安裝冶伞，參見：http://www.reibang.com/p/632f35c48137

BWA 的安裝

conda  create -n wgs  python=2 bwa # 創(chuàng)建環(huán)境及conda 安裝
source activate wgs # 激活環(huán)境
# 這樣我們安裝的包就只在這一個環(huán)境下了新症，不會使linux 工作環(huán)境污染。

samtools

conda install samtools

Picard

conda install picard

GATK

conda install gatk

sratoolkit

wget https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.10.5/sratoolkit.2.10.5-centos_linux64.tar.gz # 下載

添加內(nèi)容的環(huán)境變量响禽，vi ~/.bashrc 并在底部加上文件目錄

db-config --interactive # 按照提示完成配置

項目結(jié)構(gòu)設(shè)計

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

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

獲取E.coli K12的參考基因組序列

一般來說基因組信息都可以從ensemble 或者ncbi兩個數(shù)據(jù)庫獲取徒爹，除此之外也還有很多的數(shù)據(jù)庫，可以參考以下：
http://www.reibang.com/p/6c45620b2578

• 這里我是通過ncbi獲取芋类。

  
  
  
  

  可以通過ftp下載得到fasta 格式的序列瀑焦。

  
  

這里我們可以利用wget 獲得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

通過gzip解壓縮并將名字重命名為E.coli_K12_MG1655.fa。

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

我們可以查看一下它的前十行

$head -10 E.coli_K12_MG1655.fa 
>NC_000913.3 Escherichia coli str. K-12 substr. MG1655, complete genome
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGCTTCTGAACTG
GTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGGTCACTAAATACTTTAACCAATATAGGCATAGCGCACAGAC
AGATAAAAATTACAGAGTACACAACATCCATGAAACGCATTAGCACCACCATTACCACCACCATCACCATTACCACAGGT
AACGGTGCGGGCTGACGCGTACAGGAAACACAGAAAAAAGCCCGCACCTGACAGTGCGGGCTTTTTTTTTCGACCAAAGG
TAACGAGGTAACAACCATGCGAGTGTTGAAGTTCGGCGGTACATCAGTGGCAAATGCAGAACGTTTTCTGCGTGTTGCCG
ATATTCTGGAAAGCAATGCCAGGCAGGGGCAGGTGGCCACCGTCCTCTCTGCCCCCGCCAAAATCACCAACCACCTGGTG
GCGATGATTGAAAAAACCATTAGCGGCCAGGATGCTTTACCCAATATCAGCGATGCCGAACGTATTTTTGCCGAACTTTT
GACGGGACTCGCCGCCGCCCAGCCGGGGTTCCCGCTGGCGCAATTGAAAACTTTCGTCGATCAGGAATTTGCCCAAATAA
AACATGTCCTGCATGGCATTAGTTTGTTGGGGCAGTGCCCGGATAGCATCAACGCTGCGCTGATTTGCCGTGGCGAGAAA

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

ncbi 的SRA 梗肝，Sequence Read Archive榛瓮，數(shù)據(jù)庫http://www.ncbi.nlm.nih.gov/sra/ 是目前最常用的存儲測序數(shù)據(jù)的數(shù)據(jù)庫。
E.coli K12 作為一種非常常用的模式生物巫击，已經(jīng)有很多測序信息在ncbi上了禀晓，這里選擇其中一個測序信息SRR1770413。

小插曲顷级，PE測序與SE測序

在之前的測序技術(shù)的介紹中凫乖，我們已經(jīng)知道包括illumina 在內(nèi)的NGS測序都是短讀長的測序，每次測出來的read 長度都不太長弓颈，為了盡可能的把DNA的序列片段多測出來帽芽，一邊測不夠，就測兩邊翔冀。
于是就有了一種從被測DNA 序列兩端各測序一次的模式导街，被稱為雙末端測序（Pair-End Sequencing，簡稱PE測序）纤子。

灰色的是被測序的DNA 片段搬瑰，左邊黃色與右邊藍色分別是測序出來的read1 與read2 序列，上圖中假定它們的長度都是100bp控硼。
雖然很多時候Pair-End測序還是無法完全的將這個被測DNA 的片段完全測完泽论，但它依然提供了非常有用的信息。

比如當(dāng)我們知道每一對的read1與read2 都來自于同一個DNA片段卡乾，read1 與read2 之間的距離是這個DNA片段的長度翼悴，且read1 與read2 的方向是剛好相反的。（這些信息對于后期的變異檢測等分析非常有用）

除此之外说订，read1 在fq1 文件中位置和read2 在fq2 文件中位置是相同的抄瓦，而且read ID 之間只在末尾有一個'/1'或者'/2'的差別潮瓶。

堿基礦工

對應(yīng)雙末端陶冷，自然也有單末端測序（Single End Sequecing钙姊，簡稱SE測序），它只測序其中一端埂伦。因此煞额，在使用如bwa進行比對時，in2.fq是非強制性的（所以用方括號括起來）沾谜，只有在雙末端測序時才需要添加膊毁。

Usage: bwa mem [options] <idxbase> <in1.fq> [in2.fq]

繼續(xù)下載SRA 文件

我們可以在通過sratoolkit 工具直接下載prefetch SRRxxxxx某個SRR 的run 編號。（但實測速度非常慢）
或者也可以在NCBI 的SRR 網(wǎng)站上下載相關(guān)的數(shù)據(jù)基跑。

下載完成后由于SRR 得到的是壓縮的文件婚温，我們需要把里面的read1和read2的測序數(shù)據(jù)解壓出來，變成需要的fastq格式媳否。

fastq-dump --split-files SRR1770413.1

準(zhǔn)備就緒栅螟，現(xiàn)在目錄下就有這些東西了～

$tree ..
..
|-- bin
|-- input
|   |-- E.coli_K12_MG1655.fa
|   |-- E.coli_K12_MG1655.fa.fai
|   |-- GCF_000005845.2_ASM584v2_genomic.fna.gz
|   |-- SRR1770413.1
|   |-- SRR1770413.1_1.fastq
|   `-- SRR1770413.1_2.fastq
`-- output

正式開始

• 將目錄內(nèi)容打理一下，再來看看有哪些東西篱竭。

  
  

質(zhì)控

關(guān)于質(zhì)控力图，是不可以省略的。
可以借助Trimmomatic掺逼。
參見http://www.reibang.com/p/56977fc52ecf

構(gòu)建索引

首先我們需要為該序列構(gòu)建一個索引吃媒，以便能夠在序列比對的時候進行快速的搜索和定位。這里使用samtools吕喘。

samtools faidx E.coli_K12_MG1655.fa

根據(jù)samtools 幫助界面可知赘那，faidx 還可以對序列內(nèi)容進行提取。

faidx index/extract FASTA

因而我們也可以用它來提取特定區(qū)域的序列內(nèi)容氯质。

$samtools faidx E.coli_K12_MG1655.fa NC_000913.3:1000000-1000200
>NC_000913.3:1000000-1000200
GTGTCAGCTTTCGTGGTGTGCAGCTGGCGTCAGATGACAACATGCTGCCAGACAGCCTGA
AAGGGTTTGCGCCTGTGGTGCGTGGTATCGCCAAAAGCAATGCCCAGATAACGATTAAGC
AAAATGGTTACACCATTTACCAAACTTATGTATCGCCTGGTGCTTTTGAAATTAGTGATC
TCTATTCCACGTCGTCGAGCG

序列比對

通過NGS 儀器測序得到的漓概，而NGS 在測定這些序列的過程中，非常粗暴地將序列給沖散成一個個短序列(read)病梢，保存在FASTQ 文件中胃珍。因此文件中不同的read 之間是沒有位置關(guān)系的。

因此首先要做的就是通過序列比對的方式蜓陌，將這些read 和物種的參考基因組（比如人的全基因組）比較觅彰，從而找到不同read 在參考基因組上本來的位置，將它們捋順排好钮热。

在這次實踐中填抬，我們就需要將獲得的E.coli 測序信息，與原始參考基因組序列進行比對隧期，確定每一個read的位置飒责。

即便是相比起雙序列比對已經(jīng)有所優(yōu)化的BLAST 算法赘娄，面對海量的短數(shù)據(jù)時，執(zhí)行的耗費的時間也太長了宏蛉。而這里使用的便是BWA遣臼，其將BW(Burrows-Wheeler)壓縮算法和后綴樹相結(jié)合，可以耗費較小的時間與空間拾并，獲得較為準(zhǔn)確的序列比對結(jié)果揍堰。
關(guān)于blast 的知識：BLAST 序列比對

使用bwa 構(gòu)建索引

首先需要為參考序列構(gòu)建bwa 比對所需要的FM-index 比對索引。

$bwa index E.coli_K12_MG1655.fa
[bwa_index] Pack FASTA... 0.03 sec
[bwa_index] Construct BWT for the packed sequence...
[bwa_index] 2.27 seconds elapse.
[bwa_index] Update BWT... 0.04 sec
[bwa_index] Pack forward-only FASTA... 0.02 sec
[bwa_index] Construct SA from BWT and Occ... 1.26 sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa index E.coli_K12_MG1655.fa
[main] Real time: 3.904 sec; CPU: 3.625 sec

由于該序列比較短嗅义，所以很快（4s 不到屏歹，相比來說人類基因組就得要3h了）索引就創(chuàng)建完畢了。

這時文件中便會多出來幾個**.fa.** 的文件之碗。ps：其中E.coli_K12_MG1655.fa.fai 是我們之前用samtools 工具獲得的蝙眶。

$ls -l
total 12544
-rw-r--r-- 1 root root 4699745 May 24 10:48 E.coli_K12_MG1655.fa
-rw-r--r-- 1 root root      12 May 27 16:17 E.coli_K12_MG1655.fa.amb
-rw-r--r-- 1 root root      98 May 27 16:17 E.coli_K12_MG1655.fa.ann
-rw-r--r-- 1 root root 4641732 May 27 16:17 E.coli_K12_MG1655.fa.bwt
-rw-r--r-- 1 root root      29 May 24 10:54 E.coli_K12_MG1655.fa.fai
-rw-r--r-- 1 root root 1160415 May 27 16:17 E.coli_K12_MG1655.fa.pac
-rw-r--r-- 1 root root 2320880 May 27 16:17 E.coli_K12_MG1655.fa.sa

接著我們使用mem BWA-MEM algorithm 這個算法，將E.coli K12質(zhì)控后的測序數(shù)據(jù)定位到其參考基因組上褪那。ps：一般MEM模塊 適合長read的比對幽纷，目的是解決第三代測序技術(shù)這種能夠產(chǎn)生長達幾十kb甚至幾Mbp的read情況。而通常小于等于70bp 的情況下武通，一般使用ALN模塊霹崎。
為了方便顯示，就不使用代碼塊而轉(zhuǎn)為引用冶忱。

$bwa mem -t 4 -R '@RG\tID:foo\tPL:illumina\tSM:E.coli_K12' ~/20200524_wgs_practice/input/E.col/fasta/E.coli_K12_MG1655.fa ~/20200524_wgs_practice/input/E.col/fastq/SRR1770413.1_1.fastq ~/20200524_wgs_practice/input/E.col/fastq/SRR1770413.1_2.fastq|samtools view -Sb - > ~/20200524_wgs_practice/output/E_coli_K12.bam && echo "** bwa mapping done **"

解讀bwa 語法

由bwa 說明界面可知：

Usage: bwa mem [options] <idxbase> <in1.fq> [in2.fq]

[options] 為一系列可選擇項尾菇，比如這里的-t 4, 表示四線程，-R 則指的是Read Group的字符串信息囚枪，以派诬。而<idxbase>輸入的便是參考基因組的bw 索引文件，便是之前通過bwa 操作建立的索引链沼。

• 特別強調(diào)的-R
@RG\tID:foo\tPL:illumina\[tSM:E.coli_K12'](tsm:E.coli_K12') 在比對語法中默赂，我們以@RG開頭，表示為對比對的read 進行分組括勺。關(guān)于read group有以下幾個參數(shù)：

1）ID缆八。一般為測序的lane ID。
2）PL疾捍，為測序的平臺奈辰，需要準(zhǔn)確填寫。在GATK中乱豆，PL只被允許寫為ILLUMINA,SLX,SOLEXA,SOLID,454,LS454,COMPLETE,PACBIO,IONTORRENT,CAPILLARY,HELICOS或UNKNOWN 等信息之一奖恰。如果這一步?jīng)]有設(shè)定正確的PL名稱，則后續(xù)使用GATK過程則可能會出現(xiàn)報錯。
3）SM：為樣本ID瑟啃，因為我們在測序時论泛，由于樣本量較多，可能會分成許多不同的lane 被分別測出來蛹屿，這時候可以使用SM 區(qū)分這些樣本屁奏。
4）LB：測序文庫的名字。一般如果ID 足夠用于區(qū)分蜡峰，則可以不設(shè)定LB了袁。
對于序列比對而言朗恳，設(shè)定這四個參數(shù)就足夠了湿颅。除此之外，在RG中還需要用制表符\t將各個內(nèi)容區(qū)分開粥诫。

而<in1.fq>與[in2.fq] 則表示質(zhì)控后的fastq 文件（也就是參考序列）油航，因為這是雙末端測序，因此其結(jié)果有兩個文件怀浆。

管道后的內(nèi)容

而其中管道|后面的內(nèi)容谊囚，則是將對比數(shù)據(jù)引至samtools工具，通過samtools 將文件轉(zhuǎn)換為bam 格式（sam格式的二進制形式）执赡，其中samtools view -b便表示輸出為bam镰踏。另外-Sb與 >中間的-表示被引流的數(shù)據(jù)。接著重定向> 輸出到文件并保存沙合。

對比對結(jié)果進行排序

samtools sort -@ 4 -m 4G -O bam -o E_coli_K12.bam E_coli_K12.sorted.bam
其中-@ 4 表示排序時的線程奠伪，設(shè)定為4。
-m 表示限制排序時最大的內(nèi)存消耗首懈，設(shè)定為4G绊率。
-O bam 表示指定輸出格式為bam。
-o表示輸出文件名究履，為‘E_coli_K12.sorted.bam’滤否。

由于這里使用后我的會提示，Error: samtools sort: couldn't allocate memory for bam_mem最仑，嘗試不設(shè)置內(nèi)存藐俺。

samtools sort -O bam -o E_coli_K12.sorted.bam E_coli_K12.bam

$ll
total 481124
-rw-r--r-- 1 root root 284944979 May 27 16:41 E_coli_K12.bam
-rw-r--r-- 1 root root 207714622 May 27 18:26 E_coli_K12.sorted.bam

去除重復(fù)序列

完成了比對之后，就可以開始進行去重復(fù)了泥彤。（去除PCR 重復(fù)序列）

重復(fù)序列的產(chǎn)生

首先欲芹，什么是重復(fù)序列？

我們知道全景，在進行NGS 測序時耀石，其中一個步驟便是橋式PCR擴增。而進行該步驟目的就是為了彌補測序文庫由于構(gòu)建過程的種種問題而導(dǎo)致濃度不夠的缺陷，因而會使測序過程極有可能漏掉某些原本基因組上的一些DNA片段滞伟，導(dǎo)致測序不全揭鳞。

而正因為整個反應(yīng)都是同步進行，無差別式擴增梆奈， 所以某些本身密度就不是很低的DNA序列也會被同步的放大野崇，而這些DNA 片段就非常有可能會被重復(fù)去進行測序。

重復(fù)序列的主要影響

其主要影響在于增大了變異的結(jié)果亩钟。
1）DNA 在打斷過程可能會引發(fā)一些堿基的顛換（嘌呤與嘧啶互換）乓梨，而PCR 過程會放大該結(jié)果，以產(chǎn)生變異清酥。
2）PCR 這一過程本身可能也會帶來新的堿基錯誤扶镀。發(fā)生在前幾輪的PCR 擴增反應(yīng)發(fā)生的錯誤會在后續(xù)的過程里繼續(xù)放大，帶來新的變異焰轻。
3）對于真實的變異臭觉，PCR 反應(yīng)可能會對包含某一個堿基的DNA模版擴增而更加劇烈（PCR bias），因此如果反應(yīng)對參考等位基因的模版偏向的話辱志，可能會導(dǎo)致變異堿基的信號變小蝠筑，因而導(dǎo)致假陰性結(jié)果。

因此揩懒，重復(fù)序列對于真實的變異檢測與真實序列的檢測都有不好的影響什乙。GATK辽俗、Samtools 這些工具都會默認序列數(shù)據(jù)為非重復(fù)序列沪羔，因此重復(fù)序列的去除是必不可少的操作。（或者直接采用PCR-free 的方案）

進行重復(fù)序列的去除工作

現(xiàn)有的工具包括Samtools间涵、Picard和悦、gatk 都可以用于去除重復(fù)序列退疫。
samtools 的rmdup 是直接將重復(fù)序列從比對BAM文件中刪除。picard與gatk 的MarkDuplicates則是在BAM的FLAG信息中標(biāo)記出來鸽素，不作刪除褒繁，但我們可以在變異檢測時識別它們，并忽略馍忽。

• 我們可以用gatk 來創(chuàng)建

gatk MarkDuplicates \
-I E_coli_K12.sorted.bam\
-O E_coli_K12.sorted.markdup.bam \ 
-M E_coli_K12.sorted.markdup_metrics.txt

$tree .
.
|-- E_coli_K12.bam
|-- E_coli_K12.sorted.bam
|-- E_coli_K12.sorted.markdup.bam
`-- E_coli_K12.sorted.markdup_metrics.txt

0 directories, 4 files

• 類似的也可以使用picard

java -jar picard.jar MarkDuplicates \ 
  I=E_coli_K12.sorted.bam \
  O=E_coli_K12.sorted.markdup.bam \
  M=E_coli_K12.sorted.markdup_metrics.txt

picard 也提供了直接刪除重復(fù)序列的操作棒坏。
直接添加額外選項，設(shè)定為true 即可REMOVE_DUPLICATES=true遭笋。

創(chuàng)建序列索引

這一步的目的坝冕，是方便我們可以快速地訪問基因組上任意位置的比對情況，有助于我們隨時了解數(shù)據(jù)瓦呼。

samtools index E_coli_K12.sorted.markdup.bam

局部重比對

• 這一步與下面的BQSR 在gatk 分析E.coli K12 的本次實戰(zhàn)中都沒有應(yīng)用喂窟。原因有二。
  沒有局部重比對是因為GATK 4.0 中沒有相應(yīng)的功能項，但GATK4.0 也提供了GATK HaplotypeCaller 這個功能實現(xiàn)同樣的結(jié)果磨澡。
  而沒有BQSR 則是因為E.coli K12沒有那些必須的已知的變異集碗啄。

但我還是來了解一下重比對及BQSR 的執(zhí)行過程及它們的意義。

開始前的合并

由上圖所顯示的那樣稳摄，一般在局部重比對開始前稚字，還有一部merge操作。它的作用是將同個樣本的所有比對結(jié)果合并為一個大的bam文件厦酬〉瑁可以借助samtools 實現(xiàn)。
一般來說有些樣本測序深度會很高仗阅，其測序結(jié)果需要經(jīng)過多次測序（或者分布在不同的測序lane中）才能全部獲得昌讲。
這時候一般是先比對排序并去重后再進行合并。

samtools merge <out.bam> <in1.bam> [<in2.bam> ... <inN.bam>]

開始局部重比對

局部重比對的目的是將BWA 比對過程中所發(fā)現(xiàn)有潛在序列插入或者序列刪除（insertion和deletion霹菊，簡稱Indel）的區(qū)域進行重新校正剧蚣。這個過程往往還會把一些已知的Indel區(qū)域一并作為重比對的區(qū)域支竹。

那么為什么要進行校正這一步呢旋廷？

原因就在于參考基因組的序列特點以及bwa、bowtie 這類序列比對算法本身礼搁。這類在全局搜索序列最優(yōu)匹配的算法在存在Indel 區(qū)域的及其附近時往往會不準(zhǔn)確饶碘。【特別是當(dāng)一些存在長Indel馒吴、重復(fù)性序列的區(qū)域或者存在長串單一堿基（比如扎运，一長串的TTTT或者AAAAA等）的區(qū)域中更是如此∫粒】

另外豪治，在這些比對算法中，對于堿基錯配以及出現(xiàn)間隙的容忍度即罰分情況是不同的扯罐，而一般的序列比對算法則更傾向于堿基錯配的情況负拟。這就會導(dǎo)致基因組上本來是非常大的一段Indel 的地方，被錯誤的變?yōu)殄e配內(nèi)容與短的gap歹河。（關(guān)于雙序列比對罰分規(guī)則的影響掩浙，可以參見：https://www.yuque.com/mugpeng/bioinfo/fgk5gz#b95f85b7

這必然會讓我們檢測到很多錯誤的變異。而且秸歧，這種情況還會隨著所比對的read長度的增長（比如三代測序的Read厨姚，通常都有幾十kbp）而變得越加嚴(yán)重。

而進行序列重比對键菱，就使用了Smith-Waterman 算法谬墙，實現(xiàn)對全局比對結(jié)果的校正和調(diào)整，最大程度低地降低由全局比對算法的不足而帶來的錯誤。而且GATK的局部重比對模塊拭抬，除了應(yīng)用這個算法之外险耀，還會對這個區(qū)域中的read進行一次局部組裝，把它們連接成為長度更大的序列玖喘，這樣能夠更進一步提高局部重比對的準(zhǔn)確性甩牺。

比對先后的不同。（白色為空位累奈，黑色為讀取序列贬派，彩色為錯配）

GATK 提供了相關(guān)的功能。（推測這里作者使用的是3.x 版本的GATK）

• 首先是通過RealignerTargetCreator澎媒，定位出所有需要進行序列重比對的目標(biāo)區(qū)域搞乏。
  

java -jar GenomeAnalysisTK.jar \
 -T RealignerTargetCreator \
 -R human.fasta \
 -I sample_name.sorted.markdup.bam \
 -known 1000G_phase1.indels.b37.vcf \
 -known Mills_and_1000G_gold_standard.indels.b37.vcf \
 -o sample_name.IndelRealigner.intervals

• 接著通過IndelRealigner 對所有在第一步中找到的目標(biāo)區(qū)域運用算法進行序列重比對，最后得到新的結(jié)果戒努。也就是上圖的after 的結(jié)果请敦。

java -jar GenomeAnalysisTK.jar \
 -T IndelRealigner \
 -R human.fasta \
 -I sample_name.sorted.markdup.bam \
 -known 1000G_phase1.indels.b37.vcf \
 -known Mills_and_1000G_gold_standard.indels.b37.vcf \
 -o sample_name.sorted.markdup.realign.bam \
 --targetIntervals sample_name.IndelRealigner.intervals

需要特別說明的是，-R 參數(shù)所用的參考序列為.fasta 格式储玫。

這里還需要額外說明VCF文件侍筛。1000G_phase1.indels.b37.vcf， Mills_and_1000G_gold_standard.indels.b37.vcf分別來自于千人基因組項目和Mills項目撒穷，記錄了那些項目中檢測到的人群Indel區(qū)域匣椰。

候選的重比對，除了要在樣本自身的比對結(jié)果中尋找外端礼，還應(yīng)該把人群已知的Indel 區(qū)域也包含進來禽笑。這些文件可以在GATK bundle ftp 中下載，（ftp://ftp.broadinstitute.org/bundle/）該版本需和參考基因組版本一致蛤奥。

當(dāng)然我們之前也說過了佳镜，如果我們在后面的變異檢測使用GATK 的GATK HaplotypeCaller 模塊，則可以減少這個局部重比對過程凡桥。因為GATK HaplotypeCaller 模塊就包括了對潛在的變異區(qū)域進行相同的局部重比蟀伸。

堿基質(zhì)量重校正 BQSR

之前提到過，在上面的E_coli_K12基因組分析過程里唬血，由于沒有已知的變異集望蜡，無法進行這一步的操作。但這并不代表它不重要拷恨。

堿基質(zhì)量值是衡量測序出來的堿基準(zhǔn)確性的重要指標(biāo)脖律，而變異檢測也是極度依賴堿基質(zhì)量值。

而測定的堿基質(zhì)量值與真實結(jié)果（實際的堿基質(zhì)量）之間的偏差大小腕侄，也受到多方因素的影響小泉。包括物理和化學(xué)等對測序反應(yīng)的影響芦疏，甚至連儀器本身和周圍環(huán)境都是其重要的影響因素。

BQSR（Base Quality Score Recalibration）這一方法也就應(yīng)運而生微姊。其主要采用機器學(xué)習(xí)的方法來構(gòu)建測序堿基的錯誤率模型酸茴，并對這些堿基的質(zhì)量值進行相應(yīng)的調(diào)整。

通過對比不難發(fā)現(xiàn)兢交，經(jīng)過BQSR 校正后薪捍，橫軸的理論值 更加接近了縱軸真實值的數(shù)值，圖像也擬合為了一條直線配喳。

但是酪穿，理論值我們可以通過fastq 文件獲得，真實值是如何產(chǎn)生的晴裹？

"真實值"是采用統(tǒng)計學(xué)的方法被济，獲得極其接近的分布結(jié)果。下面是尋求“真實值”的思路涧团。

假如一個堿基的質(zhì)量報告為20只磷。根據(jù)公式Q = -10log(p_error)，則錯誤率也就是10**-2 為1%泌绣。換句話說钮追，意思就是如果有100個堿基的質(zhì)量值均為20，則其中將會有一個堿基測定結(jié)果是錯的赞别。

現(xiàn)在的問題就轉(zhuǎn)換成了如何找到這些錯誤的堿基畏陕。

通常來說人與人之間基因的整體差異是很小的，那么一般來說在一個群體中發(fā)現(xiàn)的已知變異仿滔，在某個人身上也很可能是同樣存在的。

因此我們可以嘗試排除掉所有的已知變異位點犹芹，接著計算每個報告中（測序結(jié)果）的質(zhì)量值和參考基因組比對有多少是不同的崎页，這些不同的就可以認定為是錯誤的堿基，這些不同的堿基的數(shù)目所占整體其所屬的質(zhì)量值下相同的堿基的總數(shù)目的比例腰埂，就是真實的堿基錯誤率飒焦。再通過Q = -10log(p_error)換算成質(zhì)量分數(shù)，就可以和測得的質(zhì)量數(shù)進行比較了屿笼。

不難發(fā)現(xiàn)牺荠，在上圖的原始質(zhì)量值中，大部分的質(zhì)量值都被高估了驴一。

和局部重比對一樣休雌，BQSR 也包含了兩個步驟。

1）通過BaseRecalibrator肝断，計算出所有需要進行重校正的read和特征值杈曲，并將這些信息輸出為一份校準(zhǔn)表文件（sample_name.recal_data.table）驰凛。

java -jar GenomeAnalysisTK.jar \
    -T BaseRecalibrator \
    -R human.fasta \
    -I sample_name.sorted.markdup.realign.bam \
    —knownSites 1000G_phase1.indels.b37.vcf \
    —knownSites Mills_and_1000G_gold_standard.indels.b37.vcf \
    —knownSites dbsnp_138.b37.vcf \
    -o sample_name.recal_data.table

2）然后通過PrintReads，利用第一步得到的校準(zhǔn)表文件（sample_name.recal_data.table）重新調(diào)整來BAM文件中的堿基質(zhì)量值担扑，并輸出一份新的BAM文件恰响。

java -jar /path/to/GenomeAnalysisTK.jar \
    -T PrintReads \
    -R human.fasta \
    -I sample_name.sorted.markdup.realign.bam \
    —BQSR sample_name.recal_data.table \
    -o sample_name.sorted.markdup.realign.BQSR.bam

這里必須注意，BQSR實際是為了校正測序過程中的系統(tǒng)性錯誤涌献，因此胚宦，在執(zhí)行的時候也是按照不同的測序文庫進行處理。而這時候@RG信息就顯得尤為重要了燕垃，BQSR 就是通過識別@RG中的ID信息间唉，從而分辨各個測序過程。

變異檢測

獲得最終的變異結(jié)果利术，是一般的全基因組分析的流程的目標(biāo)呈野。

變異檢測的內(nèi)容一般會包括：SNP、Indel印叁，CNV和SV等被冒。而這里使用GATK HaplotypeCaller模塊對樣本中的變異進行檢測，它也是目前最適合用于對二倍體基因組進行變異（SNP+Indel）檢測的算法轮蜕。

HaplotypeCaller 與一般的直接應(yīng)用貝葉斯推斷的算法有所不同昨悼，它會先推斷群體的單倍體組合情況，計算各個組合的幾率跃洛，然后根據(jù)這些信息再反推每個樣本的基因型組合率触。因此它不但特別適合應(yīng)用到群體的變異檢測中，而且還能夠依據(jù)群體的信息更好地計算每個個體的變異數(shù)據(jù)和它們的基因型組合汇竭。

如圖上所示葱蝗，變異檢測應(yīng)用HaplotypeCaller有兩種做法。對應(yīng)于但樣本與多樣本的情況细燎。

做法一

直接使用HaplotypeCaller两曼。這個方法適合單樣本，或者固定樣本的情況玻驻。但缺陷是悼凑，如果需要增加新的樣本就必須要重新運行HaplotypeCaller，也就是所謂N+1 難題璧瞬。

java -jar /path/to/GenomeAnalysisTK.jar \
 -T HaplotypeCaller \
 -R /path/to/human.fasta \
 -I sample_name.sorted.markdup.realign.BQSR.bam \
 -D /path/to/gatk/bundle/dbsnp_138.b37.vcf \
 -stand_call_conf 50 \ 
 -A QualByDepth \ 
 -A RMSMappingQuality \ 
 -A MappingQualityRankSumTest \ 
 -A ReadPosRankSumTest \ 
 -A FisherStrand \ 
 -A StrandOddsRatio \ 
 -A Coverage
 -o sample_name.HC.vcf

• 其中-D參數(shù)輸入的dbSNP 也同樣可以在GATK bundle 中找到户辫，這份文件匯集的是目前幾乎所有的公開人群變異數(shù)據(jù)集。

除外嗤锉，如果我們有多個.bam樣本輸入則可以繼續(xù)使用-I參數(shù)輸入渔欢。-I sample1.bam [-I sample2.bam ...]

一般來說，若是直接對人類基因組做變異檢測档冬，由于人類基因組數(shù)據(jù)過于龐大膘茎，往往整個過程會消耗大量的時間桃纯，通常需要幾十個小時甚至幾天。因此我們?yōu)榱颂岣咝士梢圆捎梅植紙?zhí)行的方法披坏。

我們可以將基因上不同的染色體理解為相互獨立的（結(jié)構(gòu)性變異除外）态坦，也就是說，我們可以對待不同的染色體棒拂，一條條獨立的來執(zhí)行變異檢測這個步驟伞梯，最后再把結(jié)果合并起來。這里就額外指定一個參數(shù)-L帚屉， 一般直接寫-L 1谜诫，但有的也會書寫-L chr1，這是由fasta 文件中的命名決定的攻旦。

做法二

第二個方法是在第一個方法之上喻旷，先將每個樣本各自生成一個gVCF，再對整個群體進行joint-genotype牢屋。

gVCF 全稱為genome VCF且预，是每個樣本用于變異檢測的中間文件，格式類似于VCF烙无，它將整個joint-genotype 過程中需要的所有信息都記錄在這里面锋谐。而從文件大小上和數(shù)據(jù)量上，都遠遠小于原來的BAM文件截酷。這樣一旦新增加樣本也不需要再重新去讀取所有人的BAM文件涮拗，只需要額外生成一份gVCF 并再進行一次joint-genotype就好。

開始變異檢測

在開始之前可以利用CreateSequenceDictonary 功能為E.coli K12的參考序列生成一個.dict 的文件迂苛。

gatk CreateSequenceDictionary -R E.coli_K12_MG1655.fa \
-O E.coli_K12_MG1655.dict

接著就開始變異檢測了三热。這里我們使用上面提到的做法二。

• 首先需要生成gVCF 的中間文件灾部。

gatk HaplotypeCaller -R E.coli_K12_MG1655.fa \
> --emit-ref-confidence GVCF \
> -I ~/20200524_wgs_practice/output/E.col/E_coli_K12.sorted.markdup.bam \
> -O ~/20200524_wgs_practice/output/E.col/E_coli_K12.g.vcf

這個過程一般會持續(xù)比較久的時間康铭，比如僅僅是E.coli_K12，就耗費了Elapsed time: 3.79 minutes赌髓。

其中--emit-ref-confidence GVCF 表示產(chǎn)生gVCF 的指令。

• 接著通過生成的gVCF檢測變異

gatk GenotypeGVCFs \
> -R ~/20200524_wgs_practice/input/E.col/fasta/E.coli_K12_MG1655.fa \
> -V E_coli_K12.g.vcf \
> -O E_coli_K12.vcf

這樣我們就獲得了E.coli K12 樣本的變異結(jié)果E_coli_K12.vcf催跪。

接著我們就可以用tabix為它構(gòu)建索引锁蠕，方便以后的分析。

tabix -p vcf

提示需要先將其壓縮懊蒸。

[tabix] was bgzip used to compress this file? E_coli_K12.vcf

壓縮一下

bgzip -f E_coli_K12.vcf

再重新執(zhí)行tabix 即可荣倾。

變異檢測質(zhì)控和過濾（VQSR）

和原始數(shù)據(jù)的質(zhì)控一樣，對變異檢測的結(jié)果進行質(zhì)控和過濾骑丸，也是同樣的重要舌仍。VQSR妒貌，Variant Quality Score Recalibration，是通過構(gòu)建GMM模型對好和壞的變異進行區(qū)分铸豁，從而實現(xiàn)對變異的質(zhì)控灌曙。它通過機器學(xué)習(xí)的方式，利用多個不同的數(shù)據(jù)特征訓(xùn)練一個模型以對變異數(shù)據(jù)進行質(zhì)控节芥。使用VQSR 需要具備兩個條件在刺。

條件1:已知變異集

首先需要一個已知的變異集，將作為訓(xùn)練質(zhì)控模型的數(shù)集⊥纺鳎現(xiàn)在有Hapmap蚣驼、OMNI，1000G和dbsnp等這些國際性項目的數(shù)據(jù)相艇，可以作為高質(zhì)量的變異集颖杏。GATK的bundle主要就是對這四個數(shù)據(jù)集做了精心的處理和選擇，然后把它們作為VQSR時的真集位點坛芽。

需要注意的是留储，VQSR并不是用這些變異集里的數(shù)據(jù)來訓(xùn)練的，而是用我們自己的變異數(shù)據(jù)靡馁。這些國際項目的數(shù)據(jù)集只是為我們提供了位點欲鹏，告訴模型群體中哪些位點存在著變異，如果在其他人的數(shù)據(jù)里能觀察到落入這個集合中的變異位點臭墨，那么這些被已知集包括的變異就有很大的可能是正確的赔嚎。

因此我們可以從數(shù)據(jù)中篩選出那些和國際項目的數(shù)據(jù)集 也就是真集相同的變異，并把它們當(dāng)作真實的變異胧弛。接著尤误，進行VQSR的時候，程序就可以用這個篩選出來的數(shù)據(jù)作為真集數(shù)據(jù)來訓(xùn)練结缚，并構(gòu)造模型了损晤。

VQSR 的原理

VQSR 的核心就在于構(gòu)建一個區(qū)分“好”與“壞”變異的分類器。

該分類器是通過GMM 模型來構(gòu)建的红竭。這個模型在構(gòu)造時尤勋，并不是使用所有的數(shù)據(jù)來進行構(gòu)造的。而是挑出一部分茵宪，那些可以和已知的變異集合位點（國際項目的參考數(shù)據(jù)集最冰，Hapmap、千人數(shù)據(jù)等）對應(yīng)的數(shù)據(jù)稀火，為它們分配不同的高可信度權(quán)重暖哨。基于群體遺傳的原理凰狞，已知的變異會被認為是更加靠譜的變異篇裁，因此在初始化的時候先把它們當(dāng)作是正確的變異沛慢。

利用這個初始變異集 ，我們可以按照特征訓(xùn)練一個區(qū)分好變異的GMM 模型达布，接著對全部數(shù)據(jù)進行打分团甲，再把那些評分低的挑選出來，構(gòu)成一個最不正確的變異集合往枣，并用這個集合構(gòu)建一個區(qū)分壞變異的GMM伐庭。

最后再同時利用好變異與壞變異的GMM 對變異數(shù)據(jù)進行打分，通過該變異的評分分冈，就可以評判這個變異的質(zhì)量值了圾另。

條件2: 檢測結(jié)果有足夠多的變異

檢測結(jié)果必須有足夠多的變異信息，不然VQSR在進行模型訓(xùn)練的時候會因為可用的變異位點數(shù)目不足而無法進行雕沉。

由于條件1的限制（已知的變異集）集乔，包括我們上面的E.coli K12 在內(nèi)，很多的非人的物種在完成變異檢測后是無法進行VQSR 的坡椒。

而條件2則會導(dǎo)致某些捕獲測序（Panel測序數(shù)據(jù)扰路，甚至外顯子測序）的數(shù)據(jù)，由于變異位點的不夠倔叼，也無法只使用VQSR汗唱。

硬過濾是什么

當(dāng)無法使用VQSR時，就不得不考慮硬過濾的方法了丈攒。

所謂硬過濾哩罪，就是人為的設(shè)定一個或者若干個指標(biāo)閾值（數(shù)據(jù)特征值），然后把所有不滿足閾值的變異位點采用一刀切的方法切除巡验。

我們該如何執(zhí)行硬過濾呢际插？或者說，該如何選擇指標(biāo)用以標(biāo)準(zhǔn)來評定變異呢显设？

我們可以直接使用gatk_vqsr 中精心挑選的指標(biāo)框弛。這些指標(biāo)一共有6個（這些指標(biāo)都在VCF文件的INFO域中，如果不是GATK得到的變異捕捂，可能會有所不同瑟枫，但知道它們的含義之后我們也是可以通過自己計算來確定的）
包括：

1.QualByDepth（QD）
2.FisherStrand (FS)
3.StrandOddsRatio (SOR)
4.RMSMappingQuality (MQ)
5.MappingQualityRankSumTest (MQRankSum)
6.ReadPosRankSumTest (ReadPosRankSum)

那么我們該如何設(shè)定這些閾值呢？

QualByDepth（QD）

QD 是變異質(zhì)量值（Quality）除以覆蓋深度（Depth）得到的比值指攒。而這里的變異質(zhì)量值就是VCF中QUAL的值力奋。QUAL 用來用來衡量變異的可靠程度；而Depth則表示這個位點上所有的含有變異堿基的樣本的覆蓋深度之和幽七。通俗來說，這個數(shù)值可以由每個含變異的樣本的覆蓋深度累加獲得溅呢。

其中澡屡，對于樣本中的G/T 非0 樣本（GT=1/1表示純合變異猿挚；GT=0/1表示雜合變異；GT=0/0表示非變異）驶鹉，計算它們的加和绩蜻。并按照定義，QD值即為質(zhì)量值QUAL 除以含有變異樣本的深度之和Depth室埋。

QD這個值描述的實際上就是單位深度的變異質(zhì)量值办绝，一般來說QD 越高變異的可信度也就越高。

但是QD 的數(shù)值是如何確定的呢姚淆？可以通過VQSR質(zhì)控先后的數(shù)據(jù)來進行解釋孕蝉。（這里選用的是NA12878 的數(shù)據(jù)）

我們可以先把所有變異位點的QD值都提取出來，然后畫一個密度分布圖腌逢。通過這個分布圖降淮，不難了解到，QD的范圍主要集中在0~40之間搏讶。而兩個峰值分別為QD = 12 以及QD = 32佳鳖。

其中第一個峰（QD=12）代表雜合，而第二峰（QD=32）代表純合媒惕。因為對于純合變異來說系吩，貢獻給質(zhì)量值的read 是雜合變異的兩倍，因而同樣深度的情況下妒蔚，QD 會更大一些穿挨。

接著我們可以畫出VQSR 后的所有可信變異（綠色），與不可信變異（紅色）的QD分布圖面睛。

由圖可見絮蒿，大部分的不可信變異都集中在了左側(cè)的低QD區(qū)域，而恰好紅色區(qū)域的波峰為QD = 2 處叁鉴。

但我們也不難發(fā)現(xiàn)土涝，在高的QD 區(qū)域其實也存在小部分的不可信變異（紅）分布；而同樣在低QD 區(qū)域也存在小部分的可信變異（綠）幌墓。而這樣的結(jié)果也是硬過濾的弊端但壮。它不會像VQSR那樣，通過多個不同維度的數(shù)據(jù)構(gòu)造合適的高維分類模型常侣，從而能夠更準(zhǔn)確地區(qū)分好和壞的變異蜡饵，而僅僅是一刀切。

FisherStrand (FS)

FS 是通過Fisher檢驗的p-value轉(zhuǎn)換而來的值胳施，它是用來描述測序或者比對時對于只含有變異的read以及只含有參考序列堿基的read是否存在著明顯的正負鏈特異性（Strand bias溯祸，或者說是差異性）

這個差異反映了測序過程是否隨機，或者說比對算法是否在基因組的某些區(qū)域存在一定的選擇傾向。如果測序過程是隨機的焦辅，而且比對不存在問題博杖，則無論read 是否存在變異，再或者該變異來自于基因組的正鏈或負鏈筷登，只要是真實的剃根，這些read 都應(yīng)該是比較均勻的，也就是說FS 的值應(yīng)該接近于零前方，即不出現(xiàn)鏈特異的比對結(jié)果狈醉。

除此之外，在VCF文件的INFO區(qū)域內(nèi)惠险，有些時候除了FS苗傅，還可能會有SB 或SOR。它們實際上也是從不同層面對鏈特異描述的指證莺匠。SB 表示原始各鏈比對數(shù)目金吗，而FS 是對結(jié)果進行了Fisher檢驗后轉(zhuǎn)換而來的結(jié)果；SOR則和FS 類似趣竣，只是使用了不同的統(tǒng)計方法來計算差異程度摇庙，它更適合高覆蓋度的數(shù)據(jù)。

而FS 硬過濾數(shù)值的確定遥缕，和QD 的確定過程是類似的卫袒。

由圖可見，F(xiàn)S 的數(shù)值范圍非常廣泛单匣，而大多數(shù)的結(jié)果還是集中在了100以內(nèi)夕凝。

StrandOddsRatio（SOR）

和FS 一樣仲吏，SOR 也是對鏈特異的一種描述不铆。由于常常FS 在硬過濾的時候不是很能滿足過濾的需求（如上所述的那樣）蝌焚，而且很多時候read 在外顯子區(qū)域末端的覆蓋存在著一定的鏈特異（是正常的），往往只有一個方向的read狂男，這個時候如果該區(qū)域中有變異位點的話综看，則FS 會給出很差的分值，而SOR 就可以起到補充校正的作用了岖食。

SOR 計算使用的是symmetric odds ratio test，數(shù)據(jù)表現(xiàn)為一個2×2的列聯(lián)表舞吭。

考慮的其實就是ALT和REF這兩個堿基的read覆蓋方向的比例是否有偏泡垃，如果完全無偏，那么應(yīng)該等于1羡鸥。

SOR = (ref_fwd/ref_rev) / (alt_fwd/alt_rev)

可以看到蔑穴，大部分的SOR 分布集中在0～3 之間。

可信變異大部分集中在0～3 之間惧浴，因此可以將閾值選擇為3存和。

RMSMappingQuality（MQ）

MQ 表示的是所有比對至該位點上的read的比對質(zhì)量值的平方和的平均值的平方根。

這個數(shù)值相比平均值衷旅，能夠更好的反應(yīng)整體的比對質(zhì)量值的離散程度捐腿。如果我們的比對工具用的是bwa mem，那么按照它的算法柿顶，對于一個好的變異位點茄袖，我們可以預(yù)期，它的MQ值將等于60嘁锯。

由圖顯示宪祥，在40以下，不可信變異幾乎完全沒有了家乘，因此可以將閾值設(shè)定為40蝗羊。

除此之外，對于MappingQualityRankSumTest（MQRankSum） 與ReadPOSRankSumTest（ReadPOSRankSum）的設(shè)定原理仁锯，和其他四個是完全一樣的耀找， 都是通過VQSR 密度分布圖，確定了硬過濾的閾值（保留大多數(shù)的可信變異扑馁，剔除大多數(shù)的不可信變異）涯呻。

執(zhí)行硬過濾

gatk 提供了VariantFiltration 模塊來進行硬過濾。但需要注意的是腻要，過濾的時候需要分SNP 與Indel 兩個不同的變異類型來進行复罐。因為這兩個變異類型的一些閾值是不同的，需要區(qū)別對待雄家。

為SNP進行硬過濾

首先要選擇出SNP

gatk SelectVariants \
    -select-type SNP \
    -V ../output/E.coli/E_coli_K12.vcf.gz \
    -O ../output/E.coli/E_coli_K12.snp.vcf.gz

接著對這些SNP 進行硬過濾效诅。

gatk VariantFiltration \
    -V ../output/E.coli/E_coli_K12.snp.vcf.gz \
    --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
    --filter-name "PASS" \
    -O ../output/E.coli/E_coli_K12.snp.filter.vcf.gz

為Indel 進行硬過濾

選擇出Indel

gatk SelectVariants \
    -select-type INDEL \
    -V ../output/E.coli/E_coli_K12.vcf.gz \
    -O ../output/E.coli/E_coli_K12.indel.vcf.gz

對Indel 進行過濾

gatk VariantFiltration \
    -V ../output/E.coli/E_coli_K12.indel.vcf.gz \
    --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
    --filter-name "PASS" \
    -O ../output/E.coli/E_coli_K12.indel.filter.vcf.gz

合并過濾結(jié)果

上面關(guān)于Indel 與SNP 的過濾所使用的閾值都是借用NA12878（來自GIAB）的高深度數(shù)據(jù)進行計算獲得的，但如果序列來源不是哺乳動物，就最好選定新的合適的閾值乱投。

gatk MergeVcfs \
    -I ../output/E.coli/E_coli_K12.snp.filter.vcf.gz \
    -I ../output/E.coli/E_coli_K12.indel.filter.vcf.gz \
    -O ../output/E.coli/E_coli_K12.filter.vcf.gz

Ti/Tv 的范圍

可以將Ti/Tv 值理解為物種在與自然相互作用和演化過程中在基因組上留下來的一個統(tǒng)計標(biāo)記咽笼，在物種中這個值具備一定的穩(wěn)定性。

Ti 是轉(zhuǎn)換 Transition戚炫， 指的是嘌呤轉(zhuǎn)嘌呤剑刑，嘧啶轉(zhuǎn)嘧啶，A-G双肤，C-T施掏。
Tv是顛換 Transversion，指的是嘌呤轉(zhuǎn)嘧啶茅糜，嘧啶轉(zhuǎn)嘌呤七芭，A-C，A-T蔑赘，G-T或 G-C狸驳。

因此在完成變異質(zhì)控之后本讥，我們最好再看一下這些變異位點Ti/Tv的值是多少相速，以此來進一步確定結(jié)果的可靠程度胸竞。

另外鳄炉，在哺乳動物基因組上C->T的轉(zhuǎn)換比較多验毡，這是因為基因組上的胞嘧啶C在甲基化的修飾下容易發(fā)生C->T的轉(zhuǎn)變陵且。

一般來說卦碾，如果沒有選擇壓力图云，Ti/Tv將等于0.5物喷，因為從概率上來說卤材，Tv 發(fā)生的概率是Ti 的兩倍。而事實上峦失，對于不同的物種扇丛，有不同的差異。對于人來說尉辑，全基因組正常的Ti/Tv在2.1左右帆精，而外顯子中則在3.0，新發(fā)生的變異則在1.5左右隧魄。（Tv 更容易發(fā)生）

Ti/Tv 實際上是一個被動指標(biāo)卓练，因此我們只能夠用它來檢驗結(jié)果的可靠程度，而不能夠?qū)⑵渥鳛橐粋€變異質(zhì)控的過濾標(biāo)準(zhǔn)购啄。

最后

其實我們都可以把各個步驟的命令寫稱shell 腳本保存起來襟企，這樣就像程序一樣，我們可以合并或者修改腳本中的某些參數(shù)狮含，達到一個快速的流程化的全基因組分析顽悼。