wes定義： 全外顯子組測序畦幢，是利用目標(biāo)序列捕獲技術(shù)贰谣， 將全基因組編碼基因外顯子區(qū)域的DNA捕獲并富集后剃执，進行高通量測序的基因組分析方法

通過測序可以找到大量的單核苷酸多態(tài)性位點(Single Nucleotide Polymorphism, SNP)，插入缺失位點(Insertion/Deletion, InDel)放典、雜合性缺失(Loss of Heterozygosity, LOH)、拷貝數(shù)變異(Copy Number Variation, CNV)【拷貝數(shù)變異只針對腫瘤】

配合比較基因組學(xué)分析基茵、群體遺傳學(xué)分析奋构、進化分析和計算生物學(xué)分析方法既可以用于深入探索疾病基因組的奧秘。

優(yōu)點：但在同樣的測序通量下拱层，外顯子組測序針對目標(biāo)區(qū)段的測序深度以及覆蓋度都高于全基因組重測序弥臼， 因此該技術(shù)對于發(fā)生在外顯子區(qū)的常見和罕見的染色體變異具有很高的檢測靈敏度。

wes流程首先要建立一個環(huán)境變量

然后在這個環(huán)境變量下安裝會用到得軟件根灯，下載比對的參考基因組和數(shù)據(jù)庫

先安裝conda径缅，用鏡像下載快一些

Miniconda是conda輕量級掺栅，我們處理數(shù)據(jù)就夠用了

conda命令詳解：http://www.reibang.com/p/8bb429659b34

Miniconda 安裝包可以到?https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/?下載。

conda create -n wes python=2??? #建立wes環(huán)境變量

conda? install? bwa? ?

conda install sra-tools

conda install samtools

conda install? ?bcftools? ?vcftools? ?snpeff

conda install? ?multiqc? ?qualimap

conda install??gatk? ?#這里需要了解gatk只是conda是不夠的

提到gatk就不得不說一下java纳猪，生物信息很多軟件都是用java寫的氧卧，所以需要在linux上配置java運行環(huán)境。平臺上的java 是1.7x的版本氏堤，但是一些生物信息軟件（比如picard）需要更高級版本的java沙绝。

安裝java

什么是GATK:?stands for?GenomeAnalysisToolkit

It is a collection of command-line tools for analyzing high-throughput sequencing data with a primary focus on variant discovery. The tools can be used individually or chained together into complete workflows. We provide end-to-end workflows, called?GATK Best Practices

gatk4.0包含了picard工具，全部是基于illumina數(shù)據(jù)格式

下載地址：http://www.broadinstitute.org/gsa/wiki/index.php/Downloading_the_GATK

這個是已經(jīng)編譯好的版本鼠锈，GATK是一個java程序闪檬，解開編譯好的壓縮包，我們可以看到下面有幾個文件：

$ ls

AnalyzeCovariates.jar?GenomeAnalysisTK.jar?resources

兩個后綴名為jar的文件和一個名為resources的文件夾购笆，這兩個.jar文件就是后面我們要用到j(luò)ava程序粗悯，主要要用的就是GenomeAnalysisTK.jar這個程序，resources文件夾里面有一些example的數(shù)據(jù)同欠，可以用來測試样傍。

在GATK使用過程中，有些步驟需要用到已知變異信息铺遂，對于這些已知變異铭乾，GATK只提供了人類的已知變異信息，可以在GATK的FTP站點下載（GATK resource bundle）

https://software.broadinstitute.org/gatk/download/bundle

GATK在進行BQSR和VQSR的過程中會使用到R軟件繪制一些圖娃循，因此炕檩，在運行GATK之前最好先檢查一下是否正確安裝了R和所需要的包，所需要的包大概包括ggplot2捌斧、gplots笛质、bitops、caTools捞蚂、colorspace妇押、gdata、gsalib姓迅、reshape敲霍、RColorBrewer等。如果畫圖時出現(xiàn)錯誤丁存，會提示需要安裝的包的名稱肩杈。

二、GATK的使用流程

GATK最佳使用方案：共3大步驟解寝。原始數(shù)據(jù)的處理—變異檢測—初步分析扩然。

1.原始數(shù)據(jù)比對和過濾，就是illumina下機數(shù)據(jù)為fastq格式聋伦，推薦使用bwa進行mapping

BWA夫偶，即Burrows-Wheeler-Alignment Tool界睁。

BWA 是一種能夠?qū)⒉町惗容^小的序列比對到一個較大的參考基因組上的軟件包。它由三個不同的算法

1.BWA-backtrack: 是用來比對 Illumina 的序列的兵拢，reads 長度最長能到 100bp翻斟。

2.BWA-SW: 用于比對 long-read ，支持的長度為 70bp-1Mbp说铃；同時支持剪接性比對访惜。

3.BWA-MEM: 都支持較長的read長度，同時都支持剪接性比對（split alignments)截汪，但是BWA-MEM是更新的算法疾牲，也更快植捎，更準(zhǔn)確衙解，且 BWA-MEM 對于 70bp-100bp 的 Illumina 數(shù)據(jù)來說，效果也更好些焰枢。選擇這個算法蚓峦，該算法先使用 MEM(maximal exact matches) 進行 seeding alignments，再使用 SW(affine-gap Smith-Waterman) 算法進行 seeds 的延伸济锄。

BWA–MEM 算法執(zhí)行局部比對和剪接性暑椰。可能會出現(xiàn) query 序列的多個不同的部位出現(xiàn)各自的最優(yōu)匹配荐绝，導(dǎo)致 reads 有多個最佳匹配位點一汽。這對 long reads 的比對時比較重要的結(jié)果。但是卻會和 Picard 的 markDuplicates 程序部兼容低滩。

特別說明：

如果 mates.fq 缺省召夹，且參數(shù) –p 未設(shè)定，那么 reads.fq 被認(rèn)為是 single-end;如果 mates.fq 存在恕沫，且參數(shù) –p 未設(shè)定监憎，那么 mem 命令會認(rèn)為 read.fq 和 mates.fq 中的 i-th reads 組成一個read對 (a read pair)，這個模式是常用的 paired-end mode婶溯。

如果參數(shù) –p 被設(shè)定鲸阔，那么， mem 命令會認(rèn)為 read.fq 中的 第 2i-th 和 第 (2i + 1)-th 的 reads 組成一個 read 對 （a read pair）迄委，這種方式也被成為交錯式的（interleaved paired-end)褐筛。 在這種情況下，即使有 mates.fq叙身，也會被忽略死讹。

bwa下載有2種方式，一種是我們前邊介紹的conda曲梗，另一種可以直接下載二進制文件然后進行編譯

$wget http://superb-sea2.dl.sourceforge.net/project/bio-bwa/bwa-0.7.12.tar.bz2

$tar jxf bwa-0.7.12.tar.bz2

$cd bwa-0.7.12

$make

$echo'PATH=$PATH:/your/bwa/path'>> ~/.bashrc

$source ~/.bashrc

應(yīng)用bwa算法之前赞警，需要利用 BWA 的 index 命令妓忍，構(gòu)建出參考基因組的 FM-index

bwa index -a??bwtsw? hg38.fa? ?#構(gòu)建參考基因組索引

構(gòu)建索引時需要注意的問題：bwa構(gòu)建索引有兩種算法，兩種算法都是基于BWT的愧旦，這兩種算法通過參數(shù)-a is 和-a bwtsw進行選擇世剖。 其中-a bwtsw對于短的參考序列是不工作的，必須要大于等于10Mb笤虫；-a is是默認(rèn)參數(shù)旁瘫，這個參數(shù)不適用于大的參考序列，必須要小于等于2G琼蚯。

最后生成文件：hg19.fa.amb酬凳、hg19.fa.ann、hg19.fa.bwt遭庶、hg19.fa.pac和hg19.fa.sa宁仔。

bwa常用參數(shù)：

-t? ?INT? 線程數(shù)，默認(rèn)是1峦睡。

-M? ? ? ? 將 shorter split hits 標(biāo)記為次優(yōu)翎苫，以兼容 Picard’s markDuplicates 軟件。

-p? ? ? ? 若無此參數(shù)：輸入文件只有1個榨了，則進行單端比對煎谍；若輸入文件有2個，則作為paired?

? ? ? ? ? ?reads進行比對龙屉。若加入此參數(shù)：則僅以第1個文件作為輸入(輸入的文件若有2個呐粘，則忽

? ? ? ? ? 略之)，該文件必須是read1.fq和read2.fa進行reads交叉的數(shù)據(jù)转捕。

-R? ? ? STR? 完整的read group的頭部作岖，可以用'\t'作為分隔符， 在輸出的SAM文件中被解釋為制

? ? ? ? ?表符TAB. read group 的ID瓜富，會被添加到輸出文件的每一個read的頭部鳍咱。

-T? INT? 當(dāng)比對的分值比 INT 小時，不輸出該比對結(jié)果与柑，這個參數(shù)只影響輸出的結(jié)果谤辜，不影響比

? ? ? ? ? 對的過程。

-a? ? ? 將所有的比對結(jié)果都輸出价捧，包括 single-end和 unpaired paired-end的 reads丑念，但是這些比

? ? ? ? ? 對的結(jié)果會被標(biāo)記為次優(yōu)。

比對命令如下：

$ bwa mem ref.fareads.fq> mem-se.sam

$ bwa mem ref.faread1.fqread2.fq> mem-pe.sam

$ bwa mem -t4-M -R"\@RG\tID:{library}\tLB:{library}\tPL:Illumina\tPU:{sample}\tSM:{sample}\"? ref.fa read1.fastq read2.fastq > mem-pe.sam 2> ./mem-pe.log

2.尋找SA coordinates结蟋，確定reads在參考基因組上的坐標(biāo)脯倚。

使用軟件：BWA

使用命令：?bwa aln

bwa aln -n 0.04 -o 1 -k 2 -t 20 reference.fa leftRead.fastq > leftRead.sai

bwa aln -n 0.04 -o 1 -k 2 -t 20 reference.fa? rightRead.fastq > rightRead.sai?

參數(shù)說明：?-n ：NUM Maximum edit distance if the value is INT, or the fraction of missing alignments given 2% uniform base error rate if FLOAT. In the latter case, the maximum edit distance is automatically chosen for different read lengths. [0.04]

-o：Maximum number of gap opens

-k：Maximum edit distance in the seed

-t：Number of threads (multi-threading mode)

將SA coordinates輸出文件轉(zhuǎn)換為sam格式

使用軟件：BWA

使用命令：bwa sample

bwa sample -a 1000 -n 3 -N 10 reference.fa leftRead.sai leftRead.fastq > leftRead.sam

bwa sample -a 1000 -n 3 -N 10 reference.fa? rightRead.sai? rightread.fastq > rightRead.sam

參數(shù)說明：?-a:Maximum insert size for a read pair to be considered being mapped properly. Since 0.4.5, this option is only used when there are not enough good alignment to infer the distribution of insert sizes.

-n:Maximum number of alignments to output in the XA tag for reads paired properly. If a read has more than INT hits, the XA tag will not be written.

-N:Maximum number of alignments to output in the XA tag for disconcordant read pairs (excluding singletons). If a read has more than INT hits, the XA tag will not be written.

將sam文件轉(zhuǎn)換為bam文件

使用軟件：samtools

使用命令：samtools view

samtools view -Sb leftRead.sam > leftRead.bam

samtools view -Sb rightRead.sam > rightRead.bam

參數(shù)說明：?-S:?input is SAM

默認(rèn)下輸入是 BAM 文件，若是輸入是 SAM 文件，則最好加該參數(shù)推正，否則有時候會報錯恍涂。

-b:?output BAM

默認(rèn)下輸出是 SAM 格式文件，該參數(shù)設(shè)置輸出 BAM 格式

但是為什么要把SAM文件轉(zhuǎn)化成BAM文件植榕，兩者都存儲了什么再沧，是什么樣的？

bwa尊残、Bowtie2是現(xiàn)下最流行的短序列比對軟件炒瘸，SAM(Sequence Alignment/Map)格式是一種通用的比對格式，用來存儲reads到參考序列的比對信息寝衫。SAM是一種序列比對格式標(biāo)準(zhǔn)顷扩，由sanger制定，是以TAB為分割符的文本格式慰毅。主要應(yīng)用于測序序列mapping到基因組上的結(jié)果表示隘截。

SAM文件官網(wǎng)：https://en.wikipedia.org/wiki/SAM_%28file_format%29

SAM做的事情：

1. 非常多序列（read)，mapping到多個參考基因組（reference）上事富；

2. 同一條序列技俐，分多段（segment）比對到參考基因組上乘陪；

3. 無限量的统台，結(jié)構(gòu)化信息表示，包括錯配啡邑、刪除贱勃、插入等比對信息；

行：除注釋外谤逼，每一行是一個read贵扰。

@HD，說明符合標(biāo)準(zhǔn)的版本流部、對比序列的排列順序戚绕；2@SQ，參考序列說明枝冀；3@RG舞丛，比對上的序列（read）說明；4@PG果漾，使用的程序說明球切；5@CO，任意的說明信息绒障。

比對結(jié)果部分（alignment section）吨凑，每一行表示一個片段（segment）的比對信息，包括11個必須的字段（mandatory fields）和一個可選的字段户辱，字段之間用tab分割鸵钝。必須的字段有11個糙臼，順序固定，不可用時恩商，根據(jù)字段定義弓摘，可以為’0’或者’*’，這是11個字段包括：

比對結(jié)果詳解：

SRR035022.2621862 163 16 59999 37 22S54M = 60102 179 CCAACCCAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCGACCCTCACCCTCACCC >AAA=>?AA>@@B@B?AABAB?AABAB?AAC@B?@AB@A?A>A@A?AAAAB??ABAB?79A?AAB;B?@?@<=8:8 XT:A:M XN:i:2 SM:i:37 AM:i:37 XM:i:0 XO:i:0 XG:i:0 RG:Z:SRR035022 NM:i:2 MD:Z:0N0N52 OQ:Z:CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCBCCCCCCBBCC@CCCCCCCCCCACCCCC;CCCBBC?CCCACCACA@

QNAME ????SRR035022.2621862

FLAG ????163

RNAME???? 16

POS ????59999

MAQ???? 37

CIGAR ????22S54M

MRNM =

MPOS???? 60102

ISIZE ????179

SEQ CCAACCCAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCGACCCTCACCCTCACCC

QUAL ????>AAA=>?AA>@@B@B?AABAB?AABAB?AAC@B?@AB@A?A>A@A?AAAAB??ABAB?79A?AAB;B?@?@<=8:8TAG XT:A:M

TAG ????XN:i:2

TAG ????SM:i:37

TAG???? AM:i:37

TAG???? XM:i:0

TAG ????XO:i:0

TAG ????XG:i:0

TAG ????RG:Z:SRR035022

TAG???? NM:i:2

TAG???? MD:Z:0N0N52

TAG OQ:Z:CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCBCCCCCCBBCC@CCCCCCCCCCACCCCC;CCCBBC?CCCACCACA

必須的字段有11個痕届，順序固定韧献，不可自行改動，根據(jù)字段定義研叫，可以為’0‘或者’*‘锤窑，這是11個字段包括：

QNAME，比對片段的（template）的編號嚷炉；

FLAG渊啰，位標(biāo)識几晤，template mapping情況的數(shù)字表示绕娘，每一個數(shù)字代表一種比對情況却音，這里的值是符合情況的數(shù)字相加總和重归；

（picard專門有一個工具解讀sam的flag:http://broadinstitute.github.io/picard/explain-flags.html）

1The read is one of a pair read是pair中的一條（read表示本條read灸拍，mate表示pair中的另一條read）

2The alignment is one end of a proper paired-end alignment pair一正一負(fù)完美的比對上

4The read has no reported alignments 這條read沒有比對上

8The read is one of a pair and has no reported alignments mate沒有比對上

16The alignment is to the reverse reference strand 這條read反向比對

32The other mate in the paired-end alignment is aligned to the reverse reference strand mate反向比對

64The read is the first (#1) mate in a pair 這條read是read1

128The read is the second (#2) mate in a pair 這條read是read2

RNAME澎蛛，參考序列的編號啡彬，如果注釋中對SQ-SN進行了定義压真，這里必須和其保持一致杆煞，另外對于沒有mapping上的序列魏宽，這里是’*‘；

POS决乎，比對上的位置队询，注意是從1開始計數(shù)，沒有比對上构诚，此處為0蚌斩；

MAPQ，mappint的質(zhì)量范嘱；

CIGAR送膳，簡要比對信息表達式（Compact Idiosyncratic Gapped Alignment Report），其以參考序列為基礎(chǔ)彤侍，使用數(shù)字加字母表示比對結(jié)果肠缨，比如3S6M1P1I4M，前三個堿基被剪切去除了盏阶，然后6個比對上了晒奕，然后打開了一個缺口，有一個堿基插入，最后是4個比對上了脑慧，是按照順序的魄眉；

“M”表示 match或 mismatch；

“I”表示 insert闷袒；

“D”表示 deletion坑律；

“N”表示 skipped（跳過這段區(qū)域）；

“S”表示 soft clipping（被剪切的序列存在于序列中）囊骤；

“H”表示 hard clipping（被剪切的序列不存在于序列中）晃择；

“P”表示 padding；

“=”表示 match也物；

“X”表示 mismatch（錯配宫屠，位置是一一對應(yīng)的）；

RNEXT滑蚯，下一個片段比對上的參考序列的編號浪蹂，沒有另外的片段，這里是’*‘告材，同一個片段坤次，用’=‘；

PNEXT斥赋，下一個片段比對上的位置缰猴，如果不可用，此處為0灿渴；

TLEN洛波，Template的長度胰舆，最左邊得為正骚露，最右邊的為負(fù)，中間的不用定義正負(fù)缚窿，不分區(qū)段（single-segment)的比對上棘幸，或者不可用時，此處為0倦零；

SEQ误续，序列片段的序列信息，如果不存儲此類信息扫茅，此處為’*‘蹋嵌，注意CIGAR中M/I/S/=/X對應(yīng)數(shù)字的和要等于序列長度；

QUAL葫隙，序列的質(zhì)量信息栽烂，格式同F(xiàn)ASTQ一樣。read質(zhì)量的ASCII編碼。

可選字段（optional fields)腺办，格式如：TAG:TYPE:VALUE焰手，其中TAG有兩個大寫字母組成，每個TAG代表一類信息怀喉，每一行一個TAG只能出現(xiàn)一次书妻，TYPE表示TAG對應(yīng)值的類型，可以是字符串躬拢、整數(shù)躲履、字節(jié)、數(shù)組等聊闯。

2.bam和sam的區(qū)別與一致

sam是帶有比對信息的序列文件（即告訴你這個reads在染色體上的位置等）崇呵，用于儲存序列數(shù)據(jù)（SAM ?format is a generic format for storing large nucleotide sequence alignments. ）。

BAM is the compressed binary version of the Sequence Alignment/Map (SAM) format. 生物信息中的二進制文件主要是為了節(jié)約空間馅袁，計算機機可讀域慷。可以用samtools工具實現(xiàn)sam和bam文件之間的轉(zhuǎn)化汗销。

二者都是fastq文件經(jīng)過序列比對或者mapping后輸出的格式（其儲存的信息都是一致的）

在UCSC種對bam文件有詳細介紹

BAM is the compressed binary version of the?Sequence Alignment/Map (SAM)?format, a compact and index-able representation of nucleotide sequence alignments. Many?next-generation sequencing and analysis tools?work with SAM/BAM. For custom track display, the main advantage of indexed BAM over PSL and other human-readable alignment formats is that only the portions of the files needed to display a particular region are transferred to UCSC. This makes it possible to display alignments from files that are so large that the connection to UCSC would time out when attempting to upload the whole file to UCSC. Both the BAM file and its associated index file remain on your web-accessible server (http, https, or ftp), not on the UCSC server

Convert SAM to BAM using the samtools program:

```samtools view -S -b -o my.bam my.sam```

Sort and create an index for the BAM:

?????samtools sort my.bam my.sorted

? ? samtools index my.sorted.bam

The sort command appends?.bam?to?my.sorted, creating a BAM file of alignments ordered by leftmost position on the reference assembly. The index command generates a new file,?my.sorted.bam.bai, with which genomic coordinates can quickly be translated into file offsets in?my.sorted.bam

在UCSC瀏覽器中顯示：

Move both the BAM file and index file (my.sorted.bam?and?my.sorted.bam.bai) to an http, https, or ftp location.

If the file URL ends with .bam, you can paste the URL directly into the?custom track?management page, click submit and view in the Genome Browser. The track name will then be the name of the file. If you want to configure the track name and descriptions, you will need to create a?track line, as shown below in next step .

Construct a?custom track?using a single?track line. The most basic version of the "track" line will look something like this:

track type=bam name="My BAM" bigDataUrl=http://myorg.edu/mylab/my.sorted.bam

Again, in addition to?http://myorg.edu/mylab/my.sorted.bam, the associated index file?http://myorg.edu/mylab/my.sorted.bam.bai?must also be available at the same location.

Paste the custom track line into the text box in the?custom track management page, click submit and view in the Genome Browser

其他信息可查看UCSC官網(wǎng)：http://genome.ucsc.edu/goldenPath/help/bam.html

通過bwa軟件對fastqc文件比對結(jié)束犹褒，前邊要做質(zhì)控（fasqc和multiqc軟件）過濾（trim_galore軟件），下一步找變異弛针，下一篇我們梳理一下GATK找變異的流程