有勁的基因資訊

        GATK ( Genome Analysis Toolkit)，一套用來分析DNA sequencing data，找尋genomic variation的工具，由Broad Institute開發。目前已被應用在幾個大計畫中，例如1000 genomes project、TCGA(The Cancer Genome Atlas)。而在今年七月也更新成第二版。

圖1. GATK流程圖

當我們將NGS 資料(reads data)利用各種方法去和參考序列(reference sequence)做比對之後，我們該如何表達比對之後的結果呢？這個問題的答案就是我們耳熟能詳的SAM檔案。SAM的縮寫是Sequence Alignment/Map，它是來自於Heng Li 等人在2009發表在Bioinformatics的文章。藉由標準的SAM檔案格式，我們得以描述每一次比對之後的結果。SAM是一個純文字檔案，可以用任何的文字編輯器開啟，其格式具有以下的優點：

1. Is flexible enough to store all the alignment information generated by various alignment programs;
2. Is simple enough to be easily generated by alignment programs or converted from existing alignment formats;
3. Is compact in file size;
4. Allows most of operations on the alignment to work on a stream without loading the whole alignment into memory;
5. Allows the file to be indexed by genomic position to efficiently retrieve all reads aligning to a locus.

簡單的說，透過SAM檔案的內容，我們可以很有彈性的去描述各種比對的狀況。此外我們也可以透過作者所提供的工具(SAMTools)來去抓取特定的區域，合併或排序比對的結果、甚至是可以根據不同的比對狀況來抓取相對應的序列資料….等等非常多的應用。

除了SAM檔案之外，我們也往往會看到BAM檔案。BAM其實就是2進位檔的SAM，因此具有檔案大小較小、存取速度較快的優點，但是缺點為無法用一般的文字編輯器開啟。由於以上特性，某些工具會僅針對BAM檔案去設計。遇到這種狀況，我們變需要先將SAM利用SAMTools去轉換成BAM，才能再往下進行相關的分析。

重複序列廣泛的出現在各物種的基因體序列中，這些序列對於生物資訊的分析上將會造成許多biases。

以人類的基因體序列為例，有將近51%的序列是屬於repeated sequence，請見圖一。

圖一： a.人類基因體序列中的repeated sequence分布情形。b. 23對染色體各自的repeated sequence分布情形。

在基因體分析上，對於研究單一個基因，可以透過它本身的基因註解，例如Gene Ontology (GO) term (可參考之前的文章，”Gene Ontology簡介”)；或是透過同源蛋白質的註解，來推測此基因的功能。然而，當研究的對象轉換為一群基因時，例如一群具有顯著表現量差異的基因時，則會得到太多基因註解的資料，造成分析困難。為了解決使問題，許多文獻透過Hypergeometric distribution，尋找具有顯著統計意義的基因註解。

首先，先簡介Hypergeometric test ，它描述了由有限個物件中抽出n個物件，成功抽出指定種類的物件的次數，並且抽出的物件不放回去。

遺傳學上將複雜疾病(complex disease)定義為多基因性(polygenic)或是多因性(multifactorial)造成的疾病，可能由多個組合的基因變異和環境因子促進而成。已知這種特性的疾病包含氣喘、乳癌、慢性白血病、肺腺癌、部分精神疾病、第二型糖尿病等等一百多種。除此之外部分的疾病還會有相同的症狀。最近新聞上有一組美國研究團隊與台灣和泰國的醫院合作發現一種俱有和愛滋病類似症狀的遺傳疾病-成人免疫缺乏症候群(http://www.nejm.org/doi/full/10.1056/NEJMoa1111160)也被認為是多因性的疾病。患有此病的人都是亞洲人，或亞洲出生、住在其他地方的人，因此該疾病可能和遺傳或是環境有關。

複雜疾病常涉及多個基因(上千或是上百個基因)，使得研究這種疾病的病理並不容易有一個很好的對象(相對於由單一或是單一組合突變所引起的疾病)。因此當在比較患者和健康者之間的基因表現差異或是 SNＰ的差異時，常會發現患有相同疾病的基因差異數量不僅很多，並且還是不同基因的多種組合。在生物網路學方面常會將高通量的實驗所發現的表現量差異作分類，挑出俱有相同表現量變化的基因做聯結，因此可以構成一個以基因為節點的大型網路。這樣的網路和蛋白質交互作用(PPI)的網路有類似的特性，基因之間的聯結代表某些關聯，而彼此之間有高密度相互關聯的基因會被稱為基因模組(module)。一般認為致病基因可能就在這些表現量差異高密度相關的模組裡。

當我們利用NGS定序儀器所產出的資料格式為FASTQ，細節的說明可參閱之前的介紹(http://yourgene.pixnet.net/blog/post/84563506)，其後續的生物資訊分析步驟中，如定序的物種有已知的基因體參考序列，則可將定序出來的short reads針對已知的基因體參考序列做alignment。

傳統的alignment分析過程，最常使用的工具就是BLAST(Basic Local Alignment Search Tool)，不過由於NGS的技術所產出的資料特點是資料量龐大，序列片段都較短，針對NGS產出的資料如使用傳統的BLAST針對已知的基因體參考序列做alignment是非常沒有效率的。因此在近幾年針對NGS的資料所設計的alignment工具大量的被發展出來。

隨著越來越多的基因體被定序，如何有系統地且大規模的收集和整理基因的功能和所在的生物路徑，顯得越來越重要。例如，可以透過Enzyme Commission number(EC number)得以了解酵素的分類和所催化的化學反應。在生物資訊領域，Gene Ontology (簡稱為GO)常用來協助了解一群感興趣的基因產物。

Gene Ontology是由Gene Ontology Consortium於1998年所發起的計畫，主要目標為訂定對於所有基因功能之分類標準，起初整合了三個模式生物的資料庫，分別為FlyBase (果蠅)、Saccharomyces Genome Database (酵母菌)和the Mouse Genome Database (小鼠)，現在已經GO的範圍已經涵蓋到原核生物和其他真核生物。

GO資料庫主要包含三個分支：

BLAST 是生物學上常用的序列比對工具，隨著BLAST 改版至BLAST+ 後，許多指令的用法已經與之前不太相同，以下就對 BLAST+ 的指令做一個簡單的介紹。

安裝:

首先請到 ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/ 下載適合平台的 BLAST 程式，一般只需要下載tar.gz 檔，解壓縮後即可使用。

Bioperl是一個基於perl程式語言之上的工具箱，他提供數個不同套件來方便生物學家分析生物資訊的資料。當使用perl分析一些檔案時，透過Bioperl僅需要幾行的程式碼就可以擷取出檔案中的各項資訊，節省一些撰寫基本程式碼的時間。儘管如此，使用Bioperl的時候，開發者仍然要相當熟悉perl語言。越熟悉perl的進階用法以及其他perl模組的搭配，越可以用Bioperl作一些變化。在許多已發表的NGS相關分析軟體中，就可以看到Bioperl的蹤影。

功能

Bioperl中的許多功能皆使用物件導向的方式來提供，也就是在使用大多功能時皆需要先建立物件才能取得他的功能。在官方網站可以看到常見的基本用法有：

我們偶爾會聽到客戶方有這樣的問題：「在 illumina 其他平臺上定序的資料是否可以跟 Hiseq 2000 上的資料一起比較？」，這樣的問題其實很微妙。

首先，我們先探討不同廠牌平臺 (如 illumina與 Roche 454) 間的比較能力，由於 Roche 454 與illumina 平臺的樣品備製方法不同，定序原理不同，在序列取樣上的方式就不同，若直接進行比較，無法得知有差異之處究竟是來自樣本間的差異，還是定序方法上之差異。

但是同一廠牌，不同平臺間的資料是否能比較呢? 以illumina 為例，由於 Illumina 在製備樣品 library 時的方式固定，並且定序平臺使用的定序策略相近，理論上應該是可以被當作重複性的實驗，為了證明這個平臺相容性的問題，在最近 3 月 8 日的 Nature ISME 期刊中，illumia 與美國的一些知名大學教授合作，發表一篇 short communication 比較 Hiseq 以及 Miseq 的定序結果，探討跨越平臺是否會產生定序結果的差異性。

在這篇文章中所研究的主題為探討微生物的族群分布並且使用 Metagenomic 的方式來分析，材料來源分別取自土壤環境或是與寄主相關的口腔、表皮和排泄物，將全部 24 隻樣品平均上在 Hiseq2000 三條 lane 上以及 Miseq 一片 flow cell 上 (一片只有一條 lane)，再將這總共四條得到的 pair-end 資料分成 5′ 和 3′ 共 8 組。上圖 1 挑選出其中 4 組做為代表，比較 Hiseq 的 lane 與 lane 之間、5′ 與 3′ 之間和Hiseq 與 Miseq 之間定序的相關性。

GC rich的區域不易定序的原因，主要發生於以下兩個階段:

1.  PCR 階段

當我們在進行 de novo 定序時，一般而言，定序深度越深 (定序量越多)組裝效果會越好，就如同統計學中所述，抽樣的樣本數越多，其分佈會越接近母體之分佈。

不過，是不是只要一直增加定序量就能完整組出 de novo 的基因體序列呢? 目前來恐怕還是件很困難的事。主要的理由在於基因體中長片段重複序列造成組裝上之問題，由於長片段重複序列被打斷時會產生許多相似的序列，使得在組裝過程中無法判斷何種組裝結果是正確的。

以下就以非常簡化的例子來說明長片段重複序列組裝上的問題。

假設一個read只包含兩個base。

有時候會遇到使用者詢問，為什麼做不同長度的mate-paired 呢? 這是因為contig間的距離不同，為了將這些不同距離的contig 組裝起來，得到更完整的組裝資訊，所以才需要使用不同長度的mate-paired。

以下就用一個簡化的例子來說明不同 mate-paired 在組裝效果上的差異:

假設我們有三個contig，這三個contig在genome上的距離如下:

在處理NGS的資料時，一開始會碰到的資料型態即是Fastq的序列格式，可以想成是fasta格式+quality值。

在前幾週的blog我們介紹了fasta的序列格式，格式如下

>sequence_name

TTAATTGGTAAATAAATCTCCTAATAGCTTAGATNTTACCTTNNNNNNNNNNTAGTTTCTTGAGATTTG

 當我們需要記錄DNA、RNA和蛋白質的序列時，我們常用一種稱為FastA的檔案格式。其檔案範例如下：

>Annotation….

ATGCGGATCGATCGA

AAACCCTGA

Velvet 是 Daniel Zerbin 與 Ewan Birney 提出的一種應用在序列組裝 (assembly) 的演算法，這個演算法是利用一種特殊的序列索引 (De Bruijn graph) 來記錄各條互相配對 read 序列，並且在此索引上依照序列 depth 以及較長序列組裝路徑來處理定序誤差及 repeat 的情況。De Bruijn graph 索引結構的設計特別適合短序列並且索引結構本身並不容易受到 read 序列庫大小而增大，同時還可記錄序列上 repeat 的片段，因此此一演算法大多應用於 Solexa 所產生的序列庫。

Velvet 在組裝序列的過程中，先將 read 序列儲存於一張網絡結構 (De Bruijn graph) 以及一張查詢表格中(上圖節錄於原作論文中，是一張非常簡化的網絡)。此一網絡結構主要是將所有 read 片段以固定長度 (K-mer，像是一個 reading frame) 抓取序列下來，單一 read 序列由多個K-mer組成，並在網絡中以一個 base 位移相互並列在一起(也就是上圖傾斜排列的序列)。如果有遇到與其配對的 read，就會把配對的 K-mer 分離出(形成上圖中藍色的區塊)，並且將原來鄰近的 K-mer 以一個箭頭連接起來。將 read 序列轉變成上圖的過程中也會將各 read 在 De Bruijn graph 起始位置以及與其它序列配對的資訊記錄在一個表格中。

在談論eQTL之前，首先要介紹一下它的由來。

最早在研究遺傳與性狀時，大多數是根據孟德爾遺傳定律，認為性狀的不同主要是來自於等位基因不同組合所造成，當然，在某些性狀上，我們可以發現確實是如此，但是在一些性狀上卻不是這樣，例如黑貓和白貓交配的後代就可能是灰色、黑白斑，而不會只有黑或白兩種結果，於是有些人開始提出多基因調控一性狀的理論，會被多基因調控的性狀就稱為數量性狀 (Quantitative trait)，而參與調控的基因則稱為數量性狀基因座 (QTL, Quantitative trait loci)。

對於初學生物資訊的不知道 DNA 序列中的「N」是什麼是很正常的，對於有經驗的人，大多也只知道「N」是代表 A 或 T 或 G 或 C，其實這個符號是來自NC-IUB (Nomenclature Committee of International Union of Biochemistry) 所建議的編號，也稱為 NC-IUB codes 或 IUB codes。

IUB codes可分為兩類，一為 Standard Bases，一為 Mixed Bases (Wobble)。

Standard Bases:

Mixed Bases (Wobble):

FPKM 是 Fragments Per Kilobase of transcript per Million mapped reads 的縮寫，FPKM 與 RPKM 其實是幾乎一樣的東西，都是用來衡量 transcripts (mRNA) 表現量的一種指標，兩者主要的不同在於計算單位上的不同，RPKM 是以 Reads 為單位，而 FPKM 是以 paired-end reads 為單位，一對 paired-end reads 視為一個 Fragment，所以 FPKM 僅適用於 paired-end 的定序資料。
在 RPKM 中，只要 read 能 map 到 transcripts，即列入計算，在 FPKM 中，只有 paired-end reads 都能 map 到 transcript (在符合 quality 條件下)，才列入計算，若 paired-end reads 中只有一個 read 能 map 到 transcript，或是其中有某個 read 的 quality 太差，都不會列入計算。

FPKM 目前主要常見於 Cufflinks 這個 RNA-Seq 分析軟體的分析結果上。

N50是一個用來評估de novo assembly效果好壞的方法之一。當我們把contig或scaffold從大到小排列後，從最大的contig開始進行長度的累加，當累加長度達到全部contig或scaffold總長度的50%時，這時所加上的contig或scaffold長度，即為N50。如下圖所示：

而N50的50，即為50%的意思。同理，N25即代表25%；N75即代表75%。數字越大，則評估條件越嚴苛。