閱讀時(shí)間：全文共3節(jié)，約6400字，約11分鐘

關(guān)鍵詞：參考序列、比對軟件、SAM文件

獲得人類基因組全外顯子組離線數(shù)據(jù)的fastq文件后，如何進(jìn)行后續(xù)的突變檢測？ 首先要做的是將測序獲得的讀數(shù)與人類基因組參考序列進(jìn)行比對。

人類基因組參考序列是如何獲得的？

隨著人類基因組計(jì)劃（Human，HGP）的進(jìn)展，Human于2001年首次公布了人類基因組序列草案。2003年，人類基因組計(jì)劃宣布完成。 HGP樣本的來源是大量匿名歐洲捐獻(xiàn)者（采集血液或精子），然后選擇少量樣本提取DNA。 為了保護(hù)捐贈(zèng)者的身份，研究人員并不知道隨后對誰的DNA進(jìn)行了測序（PMID：）。

2004年發(fā)布了更準(zhǔn)確的基因組序列，版本號(hào)是NCBI Human Bulid 35，2004年發(fā)表的一篇文章（PMID：）詳細(xì)描述了該版本基因組的組裝過程。 Bulid 35版本的基因組序列包括2.85個(gè)核酸，只有341個(gè)缺口（缺口主要是由重復(fù)片段造成），覆蓋了99%的基因組序列。

隨后，GRC（）在2009年基于NCBI Human Bulid 35版本發(fā)布了（hg19）版本，總長度為3,137,144,693； 并于2013年發(fā)布了一個(gè)版本，總長度為3,209,286,105。

目前，很多關(guān)于人類基因組的數(shù)據(jù)庫，如dbSNP等數(shù)據(jù)庫，都已經(jīng)更新到最新版本。 至于具體分析過程中使用哪個(gè)版本的參考基因組序列，您可以根據(jù)自己的需求進(jìn)行選擇。

各版本的人類基因組詳細(xì)信息請參見：

各版本詳細(xì)下載地址請參見：

＃人類）

如何使用BWA進(jìn)行序列比對？

人類基因組參考序列的來源、詳細(xì)信息和下載地址都已經(jīng)知道了，那么我們看看用什么軟件或算法將測序數(shù)據(jù)比對到31億堿基序列呢？

目前使用最廣泛的軟件非BWA莫屬。 BWA軟件除了最常見的平臺(tái)測序數(shù)據(jù)外，還可以用于SOLiD、454、reads、reads等。

Heng Li發(fā)布的BWA（-）軟件目前包含三種算法，分別是2009 BWA-算法（BWA-ALN；PMID：）、2010 BWA-SW算法（PMID：）和2013 BWA-MEM算法（arXiv： 1303.[q-bio.GN]）：

1）BWA-ALN算法可用于讀長小于或等于100 bp的測序數(shù)據(jù)；

2）BWA-SW算法可用于讀長在70 bp到1 Mb之間的測序數(shù)據(jù)；

3）BWA-MEM算法可用于讀長在70 bp到1 Mb之間的測序數(shù)據(jù)； 與BWA-SW相比，該算法更快、更準(zhǔn)確； 與BWA-ALN算法相比，它可以比較70 bp和100 bp之間的reads。 ，該算法也具有較好的性能。

在使用這三種算法之前，需要為參考序列構(gòu)建FM索引（全文索引空間）。 FM-index是一種基于全文壓縮和索引構(gòu)建的算法。

BWA 指數(shù) hg19.fa

建立FM-index后，您可以選擇三種算法之一進(jìn)行比較：

-end --- aln算法：

bwa aln hg19.fa 讀取.fastq> .sai

bwa samse hg19.fa .sai 讀取.fastq > aln.sam

-結(jié)束---bwasw算法：

bwa bwasw hg19.fa .fq > aln.sam

-結(jié)束---內(nèi)存算法：

bwa mem hg19.fa 讀取.fastq> aln.sam

-end --- aln算法：

bwa aln hg19.fa .fq> aln1.sai

bwa aln hg19.fa .fq> aln2.sai

bwa 樣本 hg19.fa aln1.sai aln2.sai .fq .fq > aln.sam

-結(jié)束---bwasw算法：

bwa bwasw hg19.fa .fq .fa> aln.sam

-結(jié)束---內(nèi)存算法：

bwa mem hg19.fa read1.fq read2.fq > aln.sam

比較后得到的信息存儲(chǔ)在SAM文件中。 盡管 SAM 文件只有 11 列必填字段，但它包含的信息量非常豐富。 接下來我們看一下該文件中存儲(chǔ)了哪些信息。

SAM文件格式簡介

SAM文件格式是Heng Li在2009年提出的用于存儲(chǔ)比較結(jié)果信息的文本文件（用Tab鍵分隔），同時(shí)他發(fā)布了處理SAM文件的軟件（PMID：）。

SAM文件包括兩部分：注釋信息( )和比較結(jié)果部分( )。

注釋信息

這里主要介紹三個(gè)信息：@SQ、@PG、@HD。

@SQ

參考序列的描述

示例：@SQ SN:chr1 LN:

SN：人類參考基因組中的染色體編號(hào)，如chr1、chr2

LN：參考序列中序列的長度

@PG

所用程序的描述

示例： @PG ID:bwa VN:0.7.13-r1126 CL: bwa mem hg19.fa read1.fq read2.fq

ID：軟件名稱

VN：軟件版本號(hào)

CL：命令行

@高清

使用 SO 記錄對齊讀取的順序。

例子：

所以：

所以：

所以：

或者sam文件中沒有@HD這行信息：比較后的默認(rèn)排序順序與比較時(shí)輸入的fastq文件一致。

那么如何對比較后的文件進(jìn)行排序呢？ 使用軟件中的排序模塊。 該模塊需要輸入bam格式：

查看 -Sb aln.sam > aln.bam

下圖展示了不同排序方式下的sam文件：

我們先看一下fastq文件中的排序順序

fastq 文件：less -SN read1.fq

對比后默認(rèn)的sam文件與對比時(shí)輸入的fastq文件一致，如下圖：

-X aln.bam | 少-SN

注意：可以看到下面文件的第二列是一個(gè)字符串，與下面的截圖不一致。 這是因?yàn)椴榭磿r(shí)添加了-X參數(shù)，可以將第二列中的FLAG轉(zhuǎn)換為字符形式。 FLAG會(huì)在比較結(jié)果的解釋中指定，這里我們主要看FLAG一欄的最后一個(gè)數(shù)字。 “1”代表末端測序中的read1，“2”代表末端測序中的read2。

按查詢名稱排序的SAM文件，按（即比較結(jié)果部分的第一列）從小到大排序：sort -n aln.bam aln。

將排序后的sam文件按與參考序列對齊的位置（即對齊結(jié)果部分的第三列和第四列）從小到大排序：sort aln.bam aln。

綜上所述，注釋部分主要記錄了比較reads時(shí)使用的參考序列的基本信息、SAM文件使用的程序以及reads的排序規(guī)則。

接下來我們詳細(xì)解釋一下SAM文件中的比較結(jié)果部分。

比較結(jié)果

比對結(jié)果部分中的每一行代表一條read的比對信息（如果是雙端測序，則每行記錄單端reads的比對情況，該read在SAM/BAM文件中稱為Query），包括11個(gè)必填字段和其他可選字段由 Tab 鍵分隔。 共有 11 個(gè)必填字段，按固定順序排列，詳細(xì)信息如下：

上校

場地

類型

簡短的

例子

名稱

讀入fastq文件

:140::4:1223:25723:72631

旗幟

INT

FLAG 數(shù)（請參閱下面的詳細(xì)信息）

聚丙稀

名稱

此讀數(shù)與參考序列中的哪條染色體進(jìn)行比對？

chr1

銷售點(diǎn)

INT

此read比對到染色體上的位置，并且是從1開始的，即參考序列的第一個(gè)堿基位置編號(hào)為1（從0開始，參考序列的第一個(gè)堿基位置編號(hào)為0）

13198

MAPQ

INT

就是Phred-：MapQ = -10 log10(P)，比如MapQ=30，表示匹配到這個(gè)位置的概率是千分之一。 與MapQ=20相比，不是隨機(jī)事件，reads比較更準(zhǔn)確。 精確的。

48

雪茄

簡要對比信息表達(dá)（詳見下文）

150M

近端串?dāng)_

示例中，末端測序的read為read1，信息為與其配對的read2相比的染色體編號(hào)（如果末端測序的read為read2，則信息為與其配對的read1編號(hào)相比的染色體編號(hào)）：“= " 表示與讀到的一致； 如果沒有，則填寫相應(yīng)的染色體編號(hào)。

PNEXT

INT

另一個(gè)讀數(shù)在染色體上的位置與末端測序中的該讀數(shù)配對。

13019

特倫

INT

根據(jù) RNAME+POS 和 RNEXT+PNEXT 計(jì)算這對 read 對應(yīng)的 DNA 文庫片段長度，如下例所示： POS – RNEXT + CIGAR 中的 M+I 編號(hào) = 13198-13019+150 = 329。則為什么是負(fù)329？ 因?yàn)閞ead1是在read2的下游對齊的。

-329

10

序列號(hào)

該讀段的堿基序列

……

11

質(zhì)量

read堿基序列對應(yīng)的質(zhì)量值

第二列是FLAG，F(xiàn)LAG是標(biāo)識(shí)符的總和，各種比較情況用不同的數(shù)字表示，每個(gè)數(shù)字代表一種比較情況。 將符合一定比較條件的查詢編號(hào)相加dnastar序列比對，得到的編號(hào)就是FLAG。 FLAG 詳細(xì)信息如下表所示：

標(biāo)志（十六進(jìn)制）

標(biāo)志（十進(jìn)制）

通過末端測序獲得的單端讀數(shù)

該標(biāo)志僅在 0x01 存在時(shí)才有意義：末端測序中的配對 read1 和 read2 與參考基因組上的適當(dāng)位置對齊。 詳細(xì)解釋見下圖1： 上圖---read1和read2對對齊到同一條染色體上； 下圖---read1和read2對齊到不同的染色體。

詳細(xì)說明1：

該讀取本身并未映射到參考基因組dnastar序列比對，請參見下面解釋 2 中的插圖。

該Flag僅在0x01存在時(shí)才有意義：與其配對的read沒有映射到參考基因組，參見下面詳細(xì)解釋2中的圖示。

詳細(xì)解釋2：

16

read本身與參考基因組的負(fù)鏈進(jìn)行比對，SAM文件中第10列的SEQ是fastq文件中堿基序列的反向互補(bǔ)。 詳細(xì)說明見3。

32

該標(biāo)志僅在 0x01 存在時(shí)才有意義：與其配對的讀取與參考基因組的負(fù)鏈對齊。

詳細(xì)解釋3：

:140::4:1101:11475:2381 的 read2 與參考基因組的負(fù)鏈對齊。 原始fastq文件中的序列是：

。 。 。 。 。 。

在SAM文件中，第十列SEQ為：

。 。 。 。 。 。

64

該標(biāo)志僅當(dāng) 0x01 存在時(shí)才有意義： -end read1

128

該標(biāo)志僅當(dāng) 0x01 存在時(shí)才有意義： -end read2

256

該比對信息不是read的最佳比對位置。 請參見解釋 4 中的插圖。

詳細(xì)說明4：

第512章

讀取失敗/

1024

該read是由PCR或()引起的read，即至少有一個(gè)其他read的堿基序列與該read一致。 PCR過程是DNA模板的復(fù)制。 如果測序數(shù)據(jù)量足夠大，就有可能檢測到具有相同堿基序列的reads。 光學(xué)重復(fù)簡單理解就是信號(hào)本身來自于測序、拍照獲取信號(hào)過程中的信號(hào)。 ，但被識(shí)別為兩個(gè)，導(dǎo)致出現(xiàn)具有相同堿基序列的兩個(gè)讀數(shù)。

FLAG如此之多，如何快速提取特定FLAG的reads呢？ 只需在查看后添加-f參數(shù)即可。 例如，要提取未映射到參考基因組的讀數(shù)：

-f參數(shù)指定FLAG對應(yīng)的十六進(jìn)制值：查看-X -f aln.bam | 少-SN

-f參數(shù)指定FLAG對應(yīng)的十進(jìn)制值：view -X -f 4 aln.bam | 少-SN

第六列CIGAR是一個(gè)簡短的比對信息表達(dá)式( )，它基于參考序列，用數(shù)字加字母表示比對結(jié)果，其中M-Match/; 我-; D-； S-軟夾；H-硬夾； P-； N-從 ，例子如下：

綜上所述，對于末端測序來說，SAM文件的比對結(jié)果部分中的每條信息不僅詳細(xì)說明了本次read的比較情況（如FLAG、RNAME、POS、CIGAR），還記錄了其匹配reads的比較情況。 情況（例如 FLAG、RNEXT、PNEXT）。

您對比對過程中涉及的人類基因組參考序列的來源、比對軟件、比對結(jié)果文件SAM/BAM是否有了更深入的了解呢？

下一期我們來談?wù)勅绾芜M(jìn)行突變檢測。

如有侵權(quán)請聯(lián)系刪除！