中 華 大 學

中 華 大 學

碩 士 論 文

使用字尾樹發展引子挑選演算法

Development of the Primers Selection Algorithm Using Suffix Trees

系 所 別：資訊工程學系碩士班 學號姓名：M09302055 謝智偉 指導教授：侯玉松 博士

中 華 民 國 九 十 六 年 二 月

使用字尾樹發展引子挑選演算法

研究生：謝智偉 指導教授：侯玉松 博士 中華大學資訊工程研究所

摘 要

DNA 晶片的應用廣泛，有鑑於 DNA 微陣列上的探針與欲檢測的檢體作雜交，

在檢測過程中需將欲檢測的檢體作標記，並使用聚合酶鏈結反應複製檢體，所以 引子的挑選左右了實驗結果的成敗，並直接反應至實驗成本上，故藉由資訊工具 在實驗前挑選較好的引子，以增進實驗成功率和降低實驗用成本。

本論文使用 Gusfield 團隊改進的 Ukkonen 字尾樹演算法，來發展引子挑選 演算法，改善以往使用雜湊演算法挑選引子導致引子長度不足的問題，更導入引 子挑選準則：引子的特異性，準確的挑選出符合實驗室所用的引子。

關鍵字：引子、字尾樹、特異性、生物資訊

Development of the Primers Selection Algorithm Using Suffix Trees

Student：Hsieh Chih－Wei Advisor：Dr. Hou Yu－Sung Institute of Computer Science and Information Engineering

Chung Hua University

ABSTRACT

The application of the biochip is extensive so far, respecting the target sequences hybridized the probe on the DNA microarray, need to label the target sequence and use polymerase chain reaction technique to amplify the target sequences when doing measuring, so selecting primers will make a success or failure of experimental result and react to the experimental costs directly, so our research method for selecting the better primers before the experiment with the computer science information tools and predict whether the experimental of result is ideal or not.

Our research use the suffix tree algorithm which created by Ukkonen and further improved by Gusfield team for develop the primers selection algorithm, we reform the problem of the primer length which is using the hash algorithm to select the candidate primers, and further consider the better primers selection criterion make the primers to fit biology experiment with the specificity of primer.

Keyword：primer, suffix tree, specificity, bioinformatic

誌 謝

經過多方的研究探討本論文才能順利完成，首先要感謝指導教授侯玉松的諄 諄教誨，指導我做研究的方法和應持有的精神，並在我怠惰時當頭棒喝，督促我 前進，以到達最終目標。

在研究期間感謝我的恩師教導我許多撰寫程式的技巧，使得我的程式功力更 上層樓，也傳授了許多報告、閱讀期刊論文和撰寫論文方面的心得；更感謝其它 生物資訊研究所的老師常幫我解開生物方面的疑惑和課業上的問題。

此外，我要感謝我的父母，他們是我背後的支柱，我前進的動力，也感謝同 儕間的互相砥礪使我能突破一些困境。

最後感謝所有曾幫助過我的人，謝謝你們。

目 錄

中文摘要 ………Ⅰ 英文摘要 ………Ⅱ 誌謝 ………Ⅲ 目錄 ………Ⅳ 圖目錄 ………Ⅵ 表目錄 ………Ⅶ

第一章 簡介………1

第二章 文獻回顧………3

2-1 微陣列系統實驗流程及其優點………3

2-2 探針的製備………4

2-2.1 聚合酶鏈鎖反應………4

2-2.2 隨機引子………5

2-2.3 最小引子集合問題………6

2-3 常用引子設計準則 ………11

第三章 字尾樹演算法 ………12

3-1 字尾樹演算法簡介 ………12

3-2 字尾樹演算法的基本定義 ………13

3-3 在字尾樹裡搜尋一段子字串 ………16

3-4 線性時間建構字尾樹 ………18

3-5 動態加入方法建構多字串字尾樹 ………26

3-6 減少字尾樹所佔空間 ………31

第四章 引子挑選演算法 ………32

4-1 利用字尾樹特性解決最小引子集合問題 ………32

4-2 基本挑選引子的方法 ………34

4-2.1 使用者輸入 ………36

4-2.2 如何以貪婪演算挑選引子方法 ………37

4-3 導入常用引子設計準則 ………40

4-3.1 更進一步改善挑選子的方法 ………40

4-3.2 引子特異性辨識 ………43

4-3.2.1 方針一：正區和負區的比重………43

4-3.2.2 方針二：以涵蓋正區數為主要挑選………45

第五章 執行效能 ………46

5-1 實驗數據與硬體規格 ………46

5-2 實驗結果與評估 ………47

5-2.1 實驗(1)：S-Primer 適用性………47

5-2.2 實驗(2)：不考慮任何挑選引子準則………48

5-2.3 實驗(3)：考慮特異性辨識對引子選取的影響………48

5-3 效能分析 ………49

第六章 結論 ………50

6-1 研究成果 ………50

6-2 未來發展方向 ………50

參考文獻 ………51

圖 目 錄

第二章

圖 2.1：DNA 晶片實驗流程 ………3

圖 2.2：聚合酶鏈鎖反應 ………5

圖 2.3：使用寡胸腺嘧啶引子的聚合酶鏈鎖反應 ………6

圖 2.4：Genome-directed primers 製作流程 ………8

圖 2.5：FASTA 檔案格式 ………8

圖 2.6：以雜湊演算法改善 FindGDPs 範例………10

第三章 圖 3.1：字串 xabxac 的字尾樹 T ………13

圖 3.2：字串 xabxa 的字尾樹 T………14

圖 3.3：字串 xabxa$的字尾樹 T………15

圖 3.4：字串 xabxa$字尾樹中的子樹 T＇的子字串 xa 包含到的葉節點 ………16

圖 3.5：字串 ababc$建構字尾樹 T 的詳細流程 ………25

圖 3.6：字串 ababc$和 bbc#以動態加入方式建構字尾樹詳細流程 ………30

圖 3.7：改善字尾樹 T 所佔的空間………31

第四章 圖 4.1：挑選候選引子的示意圖………33

圖 4.2：基本挑選引子的演算法流程圖………35

圖 4.3：候選引子和其開放讀碼區序列的資料結構示意圖………37

圖 4.4：貪婪演算法挑選候選引子步驟(1)(2)………38

圖 4.5：貪婪演算法挑選候選引子步驟(3)(4)………39

圖 4.6：改善挑選引子演算法的流程圖………41

圖 4.7：計算方針一特異性值範例………43

圖 4.8：候選引子在字尾樹涵蓋到的正區和負區示意圖………44

圖 4.9：方針二取涵蓋正區多數 ORF 為候選引子………45

表 目 錄

第五章

表一：微生物基因組內含開放讀碼區序列的數量對照表 ………46

表二：硬體規格表 ………46

表三：實驗(1)以方針一挑選引子………47

表四：實驗(1)以方針二挑選引子………47

表五：實驗(2)結果比較………48

表六：本論文實驗(3)以方針一挑選引子結果………48

表七：本論文實驗(3)以方針二挑選引子結果………49

表八：文獻[24]實驗(3)結果………49

表九：FindGDPs 實驗(3)結果………49

第一章 簡介

由於生物晶片(biochip)較傳統生化檢測快速、低廉且可大量生產等特性，

近年來吸引越來越多的專家學者，投入此領域的研究。生物晶片依製程的不同可 分作兩大主流，一是微陣列晶片(microarray)[36-37]，另一個則是微處理型晶 片(lab-on-a-chip)[38]，由於微陣列晶片在微小面積的基質上種植高密度的生 物探針(probe)，做為大量篩檢及平行分析的工具。微陣列晶片具備快速、方便、

經濟、省時等特性，適用於大量基因表達、篩檢、及比對等研究，並可以應用在 病原體基因檢測、基因表現比較、基因突變分析、基因序列分析、及新藥物開發 等領域[19]，因此微陣列晶片是目前較為廣泛應用的生物晶片；而微處理型晶片 是將實驗室中的實驗流程移至晶片中，可用來處理生物樣品、進行生物特性反應 及成份分析等，其可減少人為操作所造成的危險和污染，應用層面廣泛，是未來 生物晶片發展的方向。

DNA 微陣列係屬微陣列晶片的一種應用，依探針製作方式的不同可以分成兩 種主要的平台：(1)cDNA(complementary DNA)微陣列 、(2)寡核苷酸

(oligonucleotide)微陣列。cDNA 微陣列是使用 DNA 經過反轉錄生成的序列，即 cDNA，以精密儀器快速的將其殖株在高密度的一塊小點陣載玻片上，而玻片上的 點裡都各自有代表一個不同功能的基因，這些玻片點裡的片段序列則稱為探針。

一個晶片上可以點上一萬個以上的探針，探針長度可以數百至數千個鹼基對 (base pairs, bp)，受測樣本為 cDNA 的標的物(target)。寡核苷酸微陣列則是 利用光罩法合成寡核酸，將 DNA 的核苷酸序列，逐層打上，而其玻片上的探針長 度約為 20~25bp。兩者的應用皆是利用核苷酸序列的互補鹼基配對原理，探針序 列會與對應的基因序列雜交。

核糖核酸雜交(nucleotide hybridization)技術，利用序列與序列的互補特 性，將取自檢體之 DNA 或 RNA，以複製酵素放大並加入攜帶放射線或螢光標記的 核糖核酸，將具有標記之 DNA 加熱成為單股，再將單股 DNA 與附著於特製薄膜上 的 DNA 探針雜交成為雙股螺旋。洗去未形成雙股螺旋的單股 DNA，呈標記反應即 表示檢體含有與該探針 DNA 同源之 DNA 或 RNA。

核糖核酸雜交最困難之處，必須找到具專一性的 DNA 片段作為探針，即此 DNA 片段須在欲檢測目標序列上出現，而不會在其他非目標序列中出現。本論文 基於此困難點，進行相關研究。

聚合酶鏈鎖反應(polymerase chain reaction, PCR)的發明是基因研究上的

一大突破，由於 PCR 可以大量複製 DNA 且放大其訊號，使得生化檢測或基因定序 過程中，能夠有充足的樣本進行實驗。然而將生物晶片技術和聚合酶鏈鎖反應技 術合併，即為現今扮演基因工程等研究的重大功臣，並發展出聚合酶鏈鎖反應晶 片。

PCR 需一小段 DAN 片段做為引子(primer)，藉著特殊且專一的引子設計，PCR 可將剪制酶辨識序列加入引子，使聚合酶鏈鎖反應產物攜帶該剪制酶的辨識序 列，使無適當剪制酶辨識序列無法進行的基因工程實驗，得以順利進行。因引子 須與欲複製的 DNA 序列雜交，所以引子的專一性是考慮的因素，故透過設計一個 引子，可避免引子與實驗中的非目標序列黏合。

引子為經設計的小片段 DNA 序列，通常可以是幾個到二十多個鹼基長度的核 苷酸序列(oligo nucleotide)，可以和 cDNA 結合，以進行 PCR 或序列定序 (sequencing)。

文獻[11、15-17、21]表示想要設計一個引子，需考慮相當多的因素，在設 計或挑選引子的過程中，常因引用的生物特性和引子挑選準則的不同，而有不同 的結果出現，目前大多數用以挑選引子的工具，幾乎都具備以下引子挑選準則：

(1)引子長度、(2)引子特異性。所以能將挑選出的候選引子發揮其最大功效，更 能符合實驗的要求，且降低成本並增進實驗結果，是本論文主要研究的課題。

第二章 文獻回顧

2-1 微陣列系統實驗流程及其優點

一般標準的微陣列[36-37]系統實驗流程如圖 2.1 所示，在實驗樣本蒐集完 畢以後，接下來就是進行樣本標定，將樣本的 mRNA 經過反轉錄成 cDNA 或是試管 外轉錄(in vitro transcription)時，以螢光物質 Cy3、Cy5 或生物素(biotin) 進行標定。為了區分微陣列上的 cDNA 以及樣本標定後的 cDNA，定義微陣列點裡 的 DNA 片段為探針[19]，而欲進行實驗的樣本則稱為標的物(target)[33]；在進 行檢測之前，必須針對標的物的核苷酸序列進行處理，利用聚合酶鏈鎖反應放大 標的物，標定好的標的物與探針進行雜交，以及後續清洗處理後，再以適當的偵 測方法獲得雜交後的影像。接下來的訊號分析以及群組分類，則必須仰賴適當的 軟體以及分類運算，將微陣列所得到的大量基因訊號，轉換成有生物意義的訊 息，以便提供生物學家進行更進一步的研究。

圖 2.1：DNA 晶片實驗流程

微陣列系統最主要的優點，就是能在短時間之內得到大量基因表現的結果，

透過這些實驗，許多先前未被發現的訊息傳導路徑的基因，也獲得證實，如 MYC、

WT1 和 p53 等。

2-2 探針的製備

傳統的雜交反應，在薄膜上的核酸分子稱為標的物；加入雜交反應袋內與薄 膜行雜交反應之標記核酸稱為探針。而在基因晶片雜交反應中，於薄膜或其他承 載物上的核酸，因每一點代表一組基因，故稱其為探針；加入與之雜交的核酸，

因是由整個 mRNA 經標記作用後而來，則稱之為標的物。

在 cDNA 微陣列實驗中發生交互雜交的情況，可能是一個嚴重的問題，為了 避免交互雜交(cross hybridization)的情況發生[26]，探針挑選變得相當重要。

2-2.1 聚合酶鏈鎖反應

聚合酶鏈鎖反應是在 1983 年由 Dr. Kary B. Mullis 提出的方法，在 1985 年他與 Saiki 等人，正式發表第一篇相關論文[35]。此後更因聚合酶鏈鎖反應的 技術，使得生物界的應用一日千里，Dr. Mullis 更因此獲頒諾貝爾桂冠的殊榮。

聚合酶鏈鎖反應的原理是利用 DNA 聚合酶(DNA polymerase)對特定基因做體 外或試管內大量合成的技術，其原理步驟如下：(如圖 2.2 所示)

(1) 先要在擴增的雙股 DNA 片段上，個別設計前置引子(forward primer)及反置引子(reverse primer)。

(2) 變性作用(denaturation)：將溫度提高使雙股 DNA 變成單股。

(3) 引子接合作用(annealing)：將上述引子對與已變性的單股 DNA 配對結合。

(4) 延長作用(extension)：利用 DNA 聚合酶分別以兩股目標 DNA 做 為模板(template)，引子接合處為起點，合成新的 DNA 股。

聚合酶鏈鎖反應則是重覆上述的，變性作用、配對結合與延長作用，週而復 始的進行 DNA 量的放大動作，使得產物得以倍數的速率急速地增加。

圖 2.2：聚合酶鏈鎖反應 2-2.2 隨機引子

在較早之前，微陣列主要是使用隨機引子製備探針，然而DNA聚合酶要進行 DNA聚合反應時，一定要有引子[11、15-17、21]，因為聚合酶可以引子所提供的 3＇端－OH為起始點，進行延伸(extension)的工作。引子一般可用機器合成出 來，要進行特定DNA片段之複製時，可以根據已知DNA序列設計核酸引子，也可以 直接採用隨機引子(random primer)。所謂的隨機引子即其序列為隨機序列，不 經特別的設計，目前應用最有效果者，是以機器合成的6~8個bp的寡核苷酸混合 物。反應中若有任何一條隨機引子能結合到模板DNA上，即可引導DNA聚合反應之 進行。

雖然以隨機引子方式挑選引子，在實驗前不需先知道基因組的任何分子信 息，且在生物應用方面具有廣泛性和通用性，但其隨機後的引子長度短，因此會 有引子特異性不足的情況產生，又因真核細胞的mRNA的3＇端皆含有多腺苷酸 (poly A)，故以隨機引子方法效率不一定好；所以文獻[2]提出，當確定mRNA的3＇

端為多腺苷酸時，則可改用寡胸腺嘧啶(oligo-dT)做為引子，即為一段全為寡胸 腺嘧啶的DNA片段；使用寡胸腺嘧啶作為引子的聚合酶鏈鎖反應過程如圖2.3所 示。

圖2.3:使用寡胸腺嘧啶引子的聚合酶鏈鎖反應 2-2.3 最小引子集合問題

文獻[17]提出最小引子集合(minimal primer set, MPS)問題，是一NP-Hard 問題，而許多研究針對此問題提出各自的解決方法，簡述如下：

在2-2.2節中很明顯的，使用寡胸腺嘧啶引子做聚合酶鏈鎖反應，只能針對 3＇端有多腺苷酸的序列，對於無多腺苷酸的生物序列而言是不適用的，所以在 2000年文獻[11]提出Genome-directed primers(GDPs)的方法，並設計了

FindGDPs演算法，主要是為了改善引子挑選的不適當，其考慮了引子的特異性，

將欲挑選引子的目標序列存入資料庫中，並將序列做反轉互補動作(reverse complementary)，且搭配滑動視窗挑選候選引子的方法，以此方法找出較具特異 性的引子集合，其詳細說明如下及圖2.4所示。

(1) 將目標序列存入資料庫，如圖2.4(a)。

目標序列可在NCBI的開放資料庫中取得，生物基因組內包含所有開放 讀碼區(open reading frame, ORF)，用為引子反轉錄回cDNA的序列。

(2) 取每條目標序列3＇端的序列長度500bp，進行序列的反轉互補動作，

如圖2.4(b)。

FindGDPs要求輸入一般通用的序列格式FASTA，所以要先將FASTA中的 所以序列讀出，並取每條序列3＇端的序列長度500bp，也就是取每條 序列最後面500bp長。由NCBI開放資料庫取得的FASTA序列格式，如圖 2.5。

(3) 將經步驟(2)處理過的序列，以長度為6~8的視窗在序列上滑動，若視 窗內字元個數為8時，則挑選並記錄在反轉錄引子表中。

反轉錄引子表中的所有引子長度皆為6~8bp，即為候選引子，除了記錄 挑選到的引子，也一併儲存每條候選引子涵蓋的開放讀碼區的資訊，

以便下一步驟的再次篩選。

(4) 將反轉錄引子表中的序列與資料庫中的目標序列進行比對，出現最多 次的候選引子即為最好的引子，並刪除此引子涵蓋的目標序列，重覆 步驟(4)，直至目標序列皆被涵蓋到為止。

步驟(4)主要利用貪婪演算法(greedy method)來挑選涵蓋最多的引 子，依序更新資料庫中的目標序列，再以反轉錄引子表中涵蓋目標序 列的數量，依序的取出記在標示為較好的引子表中。

ATGATC ATGTCA

(a)：步驟1

CACTT．．．TAATGCCGT

(b)：步驟2

ACGTTA…G

500bps DB

ATGACA ATGCCT

．．．

Sequences Table

ctagcgt tcggagacgaccgc

圖2.5：FASTA檔案格式 (c)：步驟3

(d)：步驟4

圖2.4 接上頁：Genome-directed primers製作流程 Reverse Primers Table

ctagcagt

>ref|NC_004463.1|:484-1023

ATGTCTCCACTAGAAGCAATCACCCGTTGGCTGGACTCAC AAGATTTGGGCATGCAACTTGAGATCGCCGGGTTACCTTT CTAACACACGAGGCTCTCAGTAGTTGGAGCGCCGCCGAG CCAGCGTTCCTGTTGCTGGGCTTGGGTCGCCTCATCGAGG CGCTGCTTAA

>ref|NC_004463.1|:3089-3682

AATGAATTGGCCGCGACGTATAACACGCTCGTAGCTAGTC CTAACACACGAGGCTCTCAGTAGTTGGAGCGCCGCCGAG AAGCTCCGAGGCCTGAGGTGGGCTTGGGTCGCCTCTTTA tagcagtt

agcagttc etc…

agctagctgagttgcgagatc gatcgagcatcgactcagcat cagctacgggcatagcatcag agcaaacgactacgactacg actacgatcagcatcagctac

…

Primer # Covered Best Primers Table

ctagcagt 686 ctagcagt 1206 etc…

8

c tagcgtt cggagacgaccgc

8

ct agcgttc ggagacgaccgc

Reverser Primers Table ctagcagt

tagcagtt agcagttc

8

雖然 ，在 cDNA 微陣列的表現上也優於隨 機引子，但由於

作法，將會在第 引子的特異性

足的；而其引子長度短

度也是必須考慮的一個因素。

在 2003 間，主要將步驟

雜湊候選引子表 4 以貪婪

演算法挑選最好的引子。其引子長度最長可取到

18~25bp有些差距。文獻[4]執 長度，q為目標序列總數。

圖 2

(a)：以滑動視窗挑選候選引子 FindGDPs 演算法相較隨機引子為佳

FindGDPs 演算法的執行效率太差，若以龐大數量的序列，已其 2、第 3 和第 4 步驟花費相當的時間和空間，然而雖然其已考慮

，但以生物方面的引子挑選準則而言，只考慮引子的特異是顯然不

，實際實驗的引子長度最好能介於 18~25bp，故引子的長

年文獻[16]利用雜湊演算法(hash)，大幅減少FindGDPs的執行時 2 挑選出的候選引子，以雜湊函式(hashing function)建立一個

，藉由雜湊函式可以很快的找到候選引子，以方便步驟

12bp，仍離實務上的引子長度 行時間為O(pq²)，p為引子

.6 為FindGDPs以雜湊演算法挑選候選引子並建立雜湊候選引子表的範例。

滑動視窗長度 = 2 【後選引子表】

候選引子： 涵蓋目標序列數：

A T 3 T T 1 G A 2 T G 1 G T 1 5’ 3’

1. T A C C G C A T 2. A C C G G A T C 3. C C C C C A A C 4. T C G A T A T C

引子搜尋範圍 = 3

(b)：雜湊候選引子表

圖 2.6 接上頁 ndGDPs 範例

圖 2.6(a) TCCCCC 和

GATTGACT，以滑動視窗長度為 2，在目標序列 3＇端長度為 3 搜尋，即引子搜尋 範圍

有所有可能引子的陣列中，

候選引子陣列每個候選引子均指向其涵蓋目標序列的編號。

基於引子長度未達實際標準長度 18~25bp，在 2003 年文獻[17]便提出改善 子長度的 U-primer 演算法，其作法可取最長引子長度 17bp，不但可以找出最 小引

執行時間、引子長度和引子特異性，都有 幅的改善，但卻未見將實務上的引子挑選準則，導入其演算法中，有鑑於此，

在 2

最小引子集合的發展，由 候選引子

AA

3 1

： AT

：以雜湊演算法改善 Fi

中有四條目標序列，分別是 ATGCCAT、GATGGCCA、GT 為 3，取出的候選引子儲存在候選引子表中。

圖 2.6(b)將圖 2.6(a)的候選引子表經存入一具

引

子集合，更增加了引子長度因而提升了引子特異性，而 U-primer 演算允許 1~2 個鹼基的錯誤黏合(mismatch)。

以前三篇文獻來看，在挑選引子的 大

004 年，文獻[21]依舊使用雜湊方法來挑選引子和建立候選引子表，也利用 貪婪演算法挑選最小引子集合，不允許引子有錯誤黏合情況，且較具創新的加入 部份的引子挑選準則：（a）不考慮內部互補引子、(b)GC 鹼基的含量值、(c)終 端限制，也就是引子終端不允許有 G 鹼基或 C 鹼基。

綜合以上四篇的文獻，我們瞭解到挑選引子及解決 TT

TG 1 GA 2 GT 1

： GG

涵蓋序列編號 1 2 4

4 3 1 2 3

短引子、低特異性、較差的執行時間、無生物條件限制，慢慢的發展至 2004 年 雜湊

合實驗用的引子是需被依標準設定，才能達到預期的實驗結果，由文獻[23]

常 在以下簡述之。

pin loop)的二級結構或引子間彼

(4) 須不能切割到基因序列，且必須在設計之

能 (5)

鹼基對，但須注意添加鹼基對後必須查看其氨基酸 (6)

(7) 序列中

黏合至

故本論文參考了以上的引子設計準則，將以引子長度及特異性等準則，作為 引子挑選演算法的生物條件。

方法改善引子長度、縮短執行時間、提高時間、加入引子挑選準則，但以龐 大的生物基因而言，雜湊的執行時間和空間都不盡理想，所以本論文將以目前最 新的技術，以字尾樹演算法(suffix tree algorithm)取代傳統的雜湊演算法，

更進一步的提高挑選引子的執行時間、空間和正確性，更加的加入以往文獻的優 點和生物條件，以便設計出一套引子挑選演算法。

2-3 常用引子設計準則 符

得知 用的引子設計準則，

(1) 最適當的引子長度應介於 18~25bp 個核苷酸。

(2) 引子序列應避免形成 U 型迴圈(hair 此形成二聚物(dimer)。

(3) 引子 G 鹼基和 C 鹼基的含量值，依照生物體內含量值設定為佳。

引子所選用的限制酶切點必

限制酶序列的 5＇段前補 3 個核苷酸小片段(如 GCG)，以便使限制酶 完全反應切割。

引子的基因片段序列，必須計算其基因碼經接合入質體後是否仍正確 無誤，是否須添加

碼(codon)必須不能為停止碼(stop codon)UAA、UAG、UGA。

引子解鏈溫度，理想引子反應的環境溫度應介於 55~72°C。

引子特異性，主要考慮除了預計引子黏合處之外，為避免引子 含有其他與 DNA 模板中高度相似的互補區域，否則將導致引子 多處。

第三章 字尾樹演算法

字串的精確比對問題(Exact Strings Matching problem)從 1980 年開始就 入研究此領域，一直到了 1990 年因人類基因的漸漸 的被

複雜的問題上，

線 時間(linear time)建構出字尾樹的演算法，然而這個先趨的 Weiner 演算法，

稱之

不同於 Weiner 位置樹，更有效 率的建構字尾樹，和較佳的空間配置。過了多年後，直至最近，1997 年因 Ukko

及 3-1 字尾樹演算法簡介

有許多學者和專家陸續的投

解碼後，字串的精確比對便成了一項熱門且急於突破的課題，至今許多研究 和應用，尤以生物序列處理方面都參考，並使用字尾樹演算法來解決問題。

在此之前不乏有強力的演算法問世，像是 Boyer-Moore[28]和

Kunth-Morris-Pratt[29]，因建構字尾樹僅需線性時間即可完成，並可以解決上 述兩種演算法所不能解的字串精確比對問題，更可以廣泛運用在

而在 Gusfield[32]的著作裡提到，字尾樹被公認為是目前效率較佳和應用範圍 涵蓋較廣的一種演算法，並使用字尾樹演算法能夠有效的解決二十餘種應用問 題，其問題主要包括字串精確比對(Exact Strings matching)、子字串問題(The Substrings Problem)、最長共同子字串(Longest Common Substrings)、子字串 容錯問題(k-mismatch Problem)、字尾陣列(Suffix Array) … 等問題。

字尾樹演算法的概念在 1973 年由 Weiner[30]所提出的，當時是第一位以 性

為位置樹(position tree)，但因位置樹的建構相當的複雜，和之前 1970 年代的兩篇文獻過於困難，所以使得在早期字串處理方面並沒有造成太大的迴 響。1970 年代有些宣稱也能在線性時間內建構字尾樹，由於其使用的字母集 (collection alphabet)和計算模型(computational model)的不同，則建構所花 費的時間和空間均會被侷限於字母集的大小。

幾年後，於 1976 年 McCreight[27]提出一個

nen[31]參照 1976 年的 McCrreight 演算法，提出了一個擁有其所有優點，

並且突破了學者對字尾樹建構複雜的刻板觀念，所以在 Ukkonen 演算法提出 後，許多的字串精確比對都陸續的被解開，然而投入在字尾樹這塊領域的學者和 專家也漸漸變多，這樣一來更能使得字尾樹加速的演進並能處理更多有關生物 基因序列方面的課題。

3-2 字尾樹演算法的基本定義

為了描述字尾樹適用於多樣化的字串，我們使用生物序列 S 長度為 m 做為欲 來源字串也可為欲搜尋字串(pattern)、文字 (Te

以字串 S、長度為 m 產生的字尾樹 T，是一顆由根(root)開始建構 向樹(directed tree)，且包涵編號從 1 ~ m 的葉子。除了根節點 外，

舉例 生的字尾樹 T。其根節點至編號為 1

的葉節點路徑則為字串 S = xabxac，而根節點至編號為 5 的葉節點路徑則為字 串 S

建構成字尾樹的字串，字尾樹的 xt) … 等字母集。

定義(一)：

的有

每一個內部節點(internal node)至少有兩個子節點(children)，

其節點與節點聯結的邊線(edge)，都會被標記上一串 S 字串的子字串 (substring)，且每段邊線上的子字串絕不會與其他邊線上的子字串雷 同，也就是每一個子字串均是獨立的。字尾樹最主要的一項特點，當由 根節點走訪至葉節點

i

i

i

來說，圖 3.1 為字串 xabxac 所產 的子字串，也就是一個字尾為 ac。

圖 3.1：字串 xabxac 的後置樹 T

但如果刪除字串 S 若以此字串 S 建構字 尾樹則其產生的葉節點只有 3 個，因為字尾 xa(字尾 4)隱含在 xabxa(字尾 1) 裡，

的最後一個字元 c，即 S = xabxa，

則無法由根節點走訪到編號為 4 的葉節點，也就是如果 S 字串的某個字尾出 現在另一個字尾的字首(prefix)時則字尾樹無法針對每一個字尾都可以產生一 個葉節點，便會產生葉節點小於 m 的情況。

圖 3.2：字串 xabxa 的後置樹 T

為了解決此種情況，只要在字串 S 的最後加上一個獨特且唯一的字元，例如 S = xabxa$，便可以達到產生 m 個葉節點，且由根節點至每個葉節點的路徑均能 對應一個字尾，圖 3.3 為加入唯一字元＂$＂之字串 S = xabxa$ 的字尾樹 T。

圖 3.3：字串 xabxa$的字尾樹 T

3-3 在字尾樹裡搜尋一段子字串

在有了 3-2 節論述的基本字尾樹演算法概念後，我們首要探討的是如何以字

。舉例來說欲在字串 S = xabxa$中搜尋是 否包

節點 u，也就是第一次走訪至內部節點 u 並得知 P 為 S 的子字串，則再視 u 為根 節點

在字 尾樹演算法來解決字串精確比對問題

含著子字串 P = xa 則只要先將字串 S 建構成字尾樹 T，然後在字尾樹 T 中 從根節點往下走訪每個內部節點和外部節點，也就是當走訪過的邊線集合的路徑 等於子字串 P，即可得知此欲搜尋片段 P 為字串 S 的子字串，反之，當走訪整棵 字尾樹卻無路徑相對應於欲搜尋片段 P，即可得知在字串 S 裡無子字串 P 的存在。

若已在字尾樹 T 中搜尋到子字串 P = xa 時，且走訪到的路徑為根節點到內部

，而 u 節點以下即為字尾樹 T 的子樹 T＇，然後利用深先搜尋演算法

(depth-first search algorithm)走再一次的走訪子樹 T＇，便可得知子字串 P = xa 包含的葉節點，而子樹 T＇所包含的葉節點的編號就表示子字串 P = xa

串 S 中的起始位置；以圖 3.4 來看子字串 P = xa 包含的葉節點編號有 1 和 4 號，

故可得知當由字尾樹 T 的根節點開始走訪至葉節點編號為 1 和 4 號時，其邊線路 徑皆包含子字串 P = xa。

圖 3.4：字串 xabxa$字尾樹中的子樹 T＇的子字串 xa 包含到的葉節點

以一般用來解決精確字串比對的演算法來說，若要在長度為 m 的字串 S 中搜

尋長度為 n 的子字串 P，搜尋子字串 P 在字尾樹 T 中全部出現的次數需花費 O(m+n) 的搜尋時間，然而只要導入字尾樹演算法的觀念，僅需先花費 O(m)的字尾樹建 構時間，將字串 S 建構成字尾樹 T，此動作即稱為前處理(preprocessing)，之 後只要花費 O(n+k)的搜尋時間，做子字串 P 在字串 S 中出現全部次數的精確比 對動作，而 k 為子字串 P 在字串 S 中出現的全部次數。然而 k 的取得則可利用深 先搜尋演算法走訪子樹 T＇，便可得知子字串 P 所包含的葉節點，所以 k 值就等 於子字串所包含的葉節點個數，k 值也就是子字串 P 在字串 S 的不同位置(index) 所出現的累加次數。如果只需找出字串 S 是否包含子字串 P 的問題，以字尾樹演 算法解決僅需花費 O(n)的搜尋時間，因為只需發現字串 S 是否包含子字串 P 所 以 k = 1；所以字尾樹演算法可節省下許多時間且高效率的完成精確比對動作。

3-4 線性時間建構字尾樹

至今最被推崇的字尾樹演算法有三個，依提出的順序分別為，Weiner、

文將重點放在擁有 McCreight 演算法的優點且較為 易懂

只有 3-2 節的基本定義是不夠的，若 以 節論述的基本定義建構字尾樹則會花費O(m³)的建構時間，所以在此導入 字尾

定義一字串 S = xα，其 x 為單一個字元，α則為字串且α可為空

，如字尾樹 T 中有一內部節點 v 其路徑為 xα，假定字尾樹 T 中存 在著

Ukkonen所使用的斷頭指標觀念與McCreight相同，但兩者不同處在於後者是 以Headi和Tai (i…m)分別的計算是否達到斷頭指標的條件，而前者即本論文使 用的

McCreight 和 Ukkonen，本論

的 Ukkonen 演算法來建構字尾樹。

為確保建構字尾樹的時間為線性時間，

3-2

樹的另一項關鍵的觀念：斷頭指標(suffix link)。

定義(二)：

字串

另一個內部節點 u 其路徑為α，則存在由內部節點 u 指向內部節點 v 的斷頭指標，標記為(u → v)。在建構字尾樹的過程中，當依序加入 欲建構的字尾，若欲建構字尾的字首中有著字尾的子字串，便可利用斷 頭指標在已建構好的樹型中，快速的找到插入位置，即省去重覆的子字 串走訪。

li

演算法，則是以生長點(grow point)的觀念取代了後者的計算方式，並隨著 生長點在樹裡的移動逐漸將完整的樹建構出來，然而詳細的理論，在此便不贅 述，直接以一例來看加入斷頭指標後的字尾樹建構流程；圖 3.5 顯示以字串S = ababc$建構字尾樹的流程，空心圓代表生長點，虛線線段則表示成斷頭指標。

(a)：步驟 1 步驟 1：

插入字尾 1 為字串 ababc$；生長點開始走訪字尾 1 的第一個字元 a，因生長 點處於根節點位置，根節點無任何分支，所以無法往下走訪，即從生長點所在的 節點位置，長出一個新的分支和節點，其分支後的邊線記作 ababc$，節點編號 成 1。生長點不移動待在根節點。

(b)：步驟 2 步驟 2：

插入字尾 2 為字串 babc$；生長點開始走訪字尾 2 的第一個字元 b，因生長

點無法經由根節點的分支向下走訪，即從生長點所在的節點位置，長出一個新的 分支和節點，其分支後的邊線記作 ababc$，節點編號成 2。生長點不移動待在根 節點。

(c)：步驟 3

步驟 3：

插入字尾 3 為字串 abc$；生長點開始走訪字尾 3 的第一個字元 a，因可從根 節點往下經由邊線 ababc$走訪，所以生長點往邊線 ababc$移動，並停在邊線上。

若生長點停在邊線而非節點，則標記 t 為走訪過的字元數，所以 t = 1。

(d)：步驟 4

步驟 4：

插入字尾 4 為字串 bc$；生長點開始走訪字尾 3 的第一個字元 b，因為已在 邊線 ababc$上走過一個字元 a，t = 1，故將字尾 4 的第一個字元 b 與所在邊線 的第二個字元 b 比較，相同則代表可繼續往下走訪，並累記 t 值，故 t 標記為 t

= 2。

(e)：步驟 5-1

步驟 5-1：

插入字尾 5 為字串 c$；生長點開始走訪字尾 5 的第一個字元 c，因已在邊線 ababc$走過兩個字 ab，t = 2；故以字尾 5 的第一個字元 c 與所在邊線的第三個 字尾 a 相比較，不同則由生長點所處位置長出新的分支和節點，並更新連接到生 長點的邊線資訊；新節點編號成 3；生長點繼續向上走訪。

(f)：步驟 5-2

步驟 5-2：

生長點向上經由邊線 ab 走訪到根節點時，此時便要利用剛走訪到的邊線資 訊，將邊線 ab 的第一個字元除去，並往下走訪字元 b；所以生長點此時停在邊 線 babc$，因生長點的移動，符合了斷頭指標的定義，所以將之前生長點所在節 點位置連結到現在生長點的節點位置，即為斷頭指標；因走訪了一個字元，則 t

= 1；再將字尾 5 的第一個字元 c 和所在邊線的第二個字元 a 相比較，不同則由 生長點所處位置長出新的分支和節點，並更新連接到生長點的邊線資訊；新節點 編號成 4；生長點繼續向上走訪。

(g)：步驟 5-3

步驟 5-3：

生長點向上經由邊線 b 走訪到根節點時，此時便要利用剛走訪到的邊線資 訊，將邊線 b 的第一個字元除去，因將字原除去後變為需走訪空字串，此情況生 長點則不移動，停在跟節點；因生長點的移動，符合了斷頭指標的定義，所以將 之前生長點所在節點位置連結到現在生長點的節點位置，即為斷頭指標；此時生 長點須從根節點再次走訪字尾 5 的第一個字元 c，因無分支可走訪，所以由生長 的所處的位置處長出新的分支和節點，並降邊線記作 c，新節點邊號成 5；因生 長點最後停在根節點，即代表字尾 5 的走訪已結束。

(h)：步驟 6

步驟 6：

插入字尾 6 為字串$；生長點開始走訪字尾 6 的第一個字元$，因根節點無任 何分支可走訪，即生長點長出新的分支和節點，並將邊線標記為$，節點編號成 6。

圖 3.5 接上頁：字串 ababc$建構字尾樹 T 詳細流程

3-5 動態加入方法建構多字串字尾樹

一般來說，在多字串中尋找單一子字串的精確比對是必須的，有些研究方法 [11、14、24]是先將多字串連結成一條長字串，每個字串與字串間插入一個唯一 的特殊字元，再建立字尾樹，若字串眾多時，則插入的唯一字元便成了一項待解 決的問題，而且以此法做欲建構字串的前處理，必須在前處理方面，花費空間來 儲存以多條字串連結成的長字串；而本論文採用更加有效率的方式，以動態加入 欲建構字串，至正在建構的字尾樹中，如此一來，每當加入一欲建構字串，其生 長點一定從根節點開始，在現有的樹型裡建構，如遇到唯一字元，如＂$＂便自 動的視之為唯一，所以用不著擔心唯一字元不夠用的窘境，又藉著已產生的斷頭 指標，做快速子字串位置定位，便可以免除字串的前處理動作，且高效率的完成 多字串的字尾樹建構。圖 3.6 即為兩字串 S1 = ababc$和 S2 = bbc#以動態加入 方式建構字尾樹的流程。

(a)：步驟 1

步驟 1：

加入字串 S2 的字尾 1 到 3-4 節已建構好的字串 S1 = ababc$的字尾樹 T 中；

生長點從樹根往下走訪字串 S2 的字尾 1 中的第一個字元 b，生長點停在內部節 點上。

(b)：步驟 2

步驟 2：

加入字串 S2 的字尾 2；生長點往下走訪字串 S2 的字尾 2 中的第一個字元 b，

因無分支可走訪，所以由生長點所在節點長出新的分支和節點，並標記邊線為 bc#，新節點編號成 2:1；冒號前的數字代表為字串 S2 產生的葉節點，冒號後的 數字代表為字串 S2 產生的第一個葉節點；生長點繼續往上走訪，因經過斷頭指 標經過邊線 b，所以此時生長點停在根節點處，因已有斷頭指標指向走訪路徑，

所以無須再次建立相同的斷頭指標；生長點繼續往下走訪字串 S2 的字尾 2 中的 第一個字元 b，所以生長點停在內部節點上。

(c)：步驟 3

步驟 3：

加入字串 S2 的字尾 3；生長點往下走訪字串 S2 的字尾 3 中的第一個字元 c，

可走向邊線 c，因停在邊線上所以標記 t = 1；生長點停止走訪。

(d)：步驟 4 步驟 4：

加入字串 S2 的字尾 4；生長點開始往下走訪字串 S2 的字尾 4 中的第一個字 元#，因已在邊線 c$走過一個字元 c，t = 1，故以字串 S2 的字尾 4 中的第一個 字元#與邊線 c$的第二個字元$相比較，不同，即由生長點所處位置長出新的分 支和節點，並更新連接生長點的邊線資訊，新節點編號成 2:2；生長點繼續往上 走訪，經由邊線 c 和斷頭指標到達根節點，此時生長點須往下走訪剛剛經過的邊 線字串 c，並將之前生長點的位置連結到目前生長點的位置，即為斷頭指標；因 生長點走已在 c$訪過一個字元 c，所以 t = 1；故以字串 S2 的字尾 4 中的第一 個字元#與邊線 c$的第二個字元$相比較，不同即由生長點所處位置長出新的分 支和節點，並更新連接生長點的邊線資訊，新節點編號成 2:3；生長點繼續往上 走訪到根結點，因生長點移動則產生一個斷頭指標；生長點再次向下走訪字串 S2 的字尾 4 中的第一個字元#，因無任何分支可走，所以由生長點所在的節點長 出新的分支和節點，將新邊線標記為#，新節點邊號成 2:4；即完成 S2 加入以字 串 S1 建構成的字尾樹 T 中。

圖 3.6 接上頁：字串 ababc$和 bbc#以動態加入方式建構字尾樹詳細流程

3-6 減少字尾樹所佔空間

雖然字尾樹演算法所佔的空間已被證明為線性空間(linear space)，當欲建 構字串數量龐大且長度很長時，相對的字尾樹的樹型必定相當龐大，其全部的節 點數更會增加。為減少字尾樹的空間運用問題，本論文便參照 Gusfield 學者與 其團隊所提出的一種改善空間使用的方法，其說明如下：

以字尾樹的資料結構角度來看，因每個節點都必須儲存一些資訊，例如：子 節點，父節點，與父節點間的邊線路徑字串…等，所以若能減少一些節點，卻不 影響字尾樹的資料結構，更能有效的搜尋片段字串 P，則可達到減少字尾樹所佔 的空間量。

故改善空間的使用度的方法便是將原本的葉節點，依葉節點的邊線路徑字串 分成兩種類型的葉節點，分別稱作 Leaf 和 Leave；Leaf 為一葉節點且邊線路徑 字串除了唯一字元外，也包含至少一個其他字元、Leave 為一葉節點且邊線路徑 字串只有包含唯一字元，然後再將 Leave 隱藏在其父節點裡，如此一來，便能減 少一些葉節點所佔據的空間，然而被隱藏 Leave 的父節點僅須記錄此葉節點是由 那一條字串產生。以圖 3.6 的(d)為例，以此方法改善字尾樹所佔的空間，便可 由圖 3.7 明顯的看出原有 10 個葉節點減少成 4 個葉節點。

圖 3.7：改善字尾樹 T 所佔的空間

第四章 引子挑選演算法

4-1 利用字尾樹特性解決最小引子集合問題

由文獻[11、16-17、21]得知，最小引子的挑選是一個 NP-Hard 的問題，且 引子長度以 18-25bps 為佳，而挑選出一個完全精確的引子，其中不容許有錯字 存在，適合使用字尾樹解決最小引子集合問題。

然而最小引子集合的挑選，最主要是要選出一個高效率的引子，也就是在微 生物基因組內涵蓋最多的開放讀碼區(ORF)，在眾多的開放讀碼區序列中，可以 挑選的引子有很多，而這些首先符合欲挑選的引子長度的引子，便稱為候選引子 (candidate primer)。

在 3-3 節題到字尾樹可找出子字串的共同字串，這特性正好可用來挑選候選 引子，在字尾樹 T 內走訪欲挑選的引子長度，而挑選到的子字串，其子樹旗下的 葉節點編號，也就是此候選引子所涵蓋到的開放讀碼區，而更以此特性可解決最 小引子集合的引子長度受限的問題。

以字串 S1 = ababc$和字串 S2 = bbc#為例，欲挑選引子長度為 2 的候選引 子，則圖 4.1 為挑選候選引子的示意圖。

圖 4.1：挑選候選引子的示意圖

由圖 4.1 可以看出由根節點向下走訪到路徑長度大於等於 2 的節點，所包含 的子字串為 bc、babc、bbc 和 ab，這些符合引子長度的子字串，即為候選引子。

當由字尾樹 T 中得到候選引子後，便可得知每個候選引子涵蓋那一條開放讀 碼區的序列，再以每個候選引子涵蓋序列數量的多寡，做貪婪演算法，進而再一 次的挑選出能覆蓋所有開放讀碼區的 3＇端序列。針對最小引子集合的詳細做法 和流程均在下節開始論述。

4-2 基本挑選引子的方法

挑選引子最主要的兩步驟，已在 4-1 節稍有提到，即為候選引子的挑選和貪 婪演算法的篩選，在做挑選引子兩大步驟之前，必須要將字串做前處理，使得挑 選後得到的引子與文獻[21]所提出的引子相同，而圖 4.2 為基本挑選引子的演算 法流程圖。

開始

1. 使用者輸入：

輸入 FASTA 檔名

圖 4.2：基本挑選引子的演算法流程圖 S-Primer

檔案格式

否 將 FASTA 序列格式轉成 S-Primer 的序列格式 是

2. 將欲建構字串 轉為逆互補序列

3. 序列集合建構 成字尾樹

4. 使用者輸入：

欲挑選引子長度

5. 經後深先搜尋 算法挑選出符合條 件的候選引子

6. 經貪婪演算法 挑選出引子

結束

4-2.1 使用者輸入

為了解決文獻[24]遇到空間上的不足和執行效率不甚快速的問題，本論文引 用 Gusfield 團隊所開發的精確字串比對程式，稱作 strmat，其主體架構在 Ukkonen 的字尾樹演算法並加以改進，而 3-6 節則為 strmat 改進 Ukkonen 字尾 樹演算法的方式與應用。本論文沿用其所有執行的優點，並針對挑選引子的問 題，更進一步的改進程式架構，在此稱本論文所撰寫的程式為 S-Primer。

一般來說，相較以往的文獻來看，大多都由 NCBI 網站的資料庫取得實驗用 的目標生物基因組，而 FASTA 序列格式取得容易且易於處理，所以許多處理生物 序列的相關程式，均以 FASTA 序列格式為首選，例如文獻[11]中的 GDPFinder 便是使用 FASTA 序列格式。但 strmat 是以自定的格式讀檔，所以當使用者讀入 FASTA 序列格式，須經過一道轉檔步驟，才能有效的建構字尾樹。

由文獻[15]隨機混合引子挑選法得知，引子主要是從序列的 5＇端開始挑 選，所以當使用者輸入字串或讀檔完成後，便需經過圖 4.2 的第二步驟，將欲建 構字串轉為逆互補序列(reverse complementary sequence)，此步驟即是將字串 反轉後再把相對應的去氧核糖核酸作轉換，也就是鹼基配對，AT 互補、GC 互補，

而此動作稱為字串的前處理。

字串的前處理完畢後，緊接著 strmat 會快速的將字尾樹建好，接著使用者 必需再次輸入欲挑選引子的長度，以文獻[14]中提到在實務上挑選的引子長度最 好能在 18-25bps 間，以程式的執行時間和效率來看，若欲挑選引子長度越長相 對的程式執行效率則會降低，所以輸入一個適當的引子長度是一件相當重要的 事。

4-2.2 如何以貪婪演算挑選引子方法

貪婪演算法是一種尋找最佳解的方法，尋找方式從某一起點開始，不斷的改 進該解答，直到無法改進為止，則最後得到的答案即為最佳解。貪婪演算法幾乎 可以解決大部份的最佳化問題，因最小引子集合問題為 NP-Hard，所以本論文則 利用貪婪演算法逐步的挑選引子，但挑選出的引子集合，無法稱作最佳解，只能 說是接近最佳解的近似值。

經過圖 4.2 基本挑選引子流程的第五個步驟挑出候選引子後，便可以得知每 個候選引子所涵蓋開放讀碼區的序列，然後以涵蓋開放讀碼區序列的多寡做貪婪 演算法，每一次貪婪法選取涵蓋最多開放讀碼區序列的候選引子為引子，同時更 新其他候選引子開放讀碼序列的數量和編號。

以圖 4.1 的字尾樹 T 為例，欲挑選引子長度為 2 的候選引子，則圖 4.3 為候 選引子和其開放讀碼區序列的資料結構示意圖。

head

圖 4.3：候選引子和其開放讀碼區序列的資料結構示意圖

bc

babc

bbc

1 2

1

2

ab 1

候選引子串列 包含候選引子的

開放讀碼區序列編號串列

有了候選引子的資料結構後，再導入貪婪演算法的觀念，以下為選取和更新 候選引子的資料結構步驟：(圖 4.4 和圖 4.5 為此四步驟的示意圖)

(1) 從候選引子串列的 head 節點開始每點走訪，在走訪的過程中可取得此 候選引涵蓋到開放讀碼區序列的數目，直到走訪完畢便可取得在候選引 子中涵蓋最多開放讀碼序列的候選引子節點。

(2) 將此候選引子節點獨立出，再走訪其開放讀碼區序列編號串列，當走訪 到每一個序列編號節點時，就同時將其他候選引子所涵蓋到同樣編號的 序列編號節點刪除，此動作稱為更新候選引子表。

(3) 然後走訪完獨立出的開放讀碼區序列編號並且更新其他候選引子後，則 此獨立出的候選引子便躍升為引子集合中的一員。

(4) 一直重覆步驟(1)、(2)、(3)，直到候選引子串列中的任一個節點已無 涵蓋任何的開放讀碼區序列。

步驟(1)、(2)：

head

babc

圖 4.4：貪婪演算法挑選候選引子步驟(1)(2)

bbc

ab

1

1 2

選取並獨立出候選引子節點

bc 1 2

步驟(3)、(4)：更新其他候選引子並選出引子

圖 4.5：貪婪演算法挑選候選引子步驟(3)(4)

babc

bbc

ab

NULL

NULL NULL head

因其他候選引子皆與 bc 同樣涵蓋編號 1 和 2 的 ORF 所以將相同編號的節點刪除

候選引子集合：bc

4-3 導入常用引子設計準則

在 4-1、4-2 節的論述中，我們還未加入一些挑選引子的準則，為了不讓挑 選出來的引子集合，在 PCR 反應中沾黏到非目標(off-target)序列，所以在以往 的實驗裡，會刻意挑選較長的引子來克服，這樣一來與目標序列接合的特異性勢 必有效的提高，但卻衍生出兩個問題：(1)引子長度越長時，所挑選的引子集合 則相對的增加，對生物實驗來說，其實驗成本必定增加。(2)長引子會降低 PCR 反應速率，故效率降低。

而本論文參考文獻[23]的一些常用引子設計準則和文獻[20]提到的特異性 辨識方法，更進一步的改善挑選引子演算法，將做到經由程式 S-Primer 挑選出 的引子更能夠符合實驗室的生物意義。

4-3.1 更進一步改善挑選引子的方法

加入 4-3 節提到的引子挑選準則，並加強引子在生物應用上的特異性，圖 4.6 即為改善挑選引子演算法的流程圖。

開始

1. 使用者輸入：

(a) 讀入 FASTA 檔案(讀檔)

(b) 欲取序列的正區、負區和忽略區(%) (4-3.2 節論述)

S-Primer 檔案格式

圖 4.6：改善挑選引子演算法的流程圖

將 FASTA 序列格式轉成 S-Primer 的序列格式 否

是

2. 將欲建構字串轉 為逆互補序列

4. 使用者輸入：

3. 序列集合建構成 字尾樹

(a) 欲挑選引子長度 (b) 引子特異性方針

5. 深先搜尋字尾樹 T 中的共同子字串

6. 符合條件的候選 引子集合

7. 經貪婪演算法挑 選出引子

引子

長度 是

是 特異性 否

結束 否

是

很明顯的可以看出使用者在基本挑選引子演算法流程中，僅需輸入欲挑選引 子的長度，在改善後的引子挑選演算法步驟 4 中，我們增加了一項更符合實驗室 應用的條件：(b)引子特異性辨識。而加入此項判斷條件即可增加候選引子的特 異性，也就是說不須如以往作法一般刻意加長引子長度來增加特異性，所以在第 5 個步驟便會略過特異性不足的候選引子，而繼續的挑選特異足夠的候選引子。

下節將詳述此項的特性。

4-3.2 引子特異性辨識

文獻[20]中提出特異性的辨識值，即引子在所有開放讀碼區前端區域的出現 數(記作 PN)，除以引子在所有開放讀碼區前端與後端區域的出現數(記作 TN)的 值，也就是特異性值等於 PN/TN，當求得的特異性值愈高則代表引子出現在開放 讀碼區前端區愈多，並使得 PCR 可複製較長的序列，所以特異性的門檻愈高，所 挑選的引子越接近實驗室中可用的引子。而引子在所有開放讀碼區後端區域的珠 現數，記作 NN。

改善挑選引子演算法流程中的第一個步驟，我們增加了一項條件(b)欲取序 列的正區、負區和忽略區，所謂的忽略區代表開放讀碼區後端區域，若引子函蓋 此區，並不會對實驗結果造成影響，正區即代表開放讀碼區前端區域，負區則代 表開放讀碼區介於正區和忽略區之間的區域。而特異性的辨識，本論文提出兩種 方針，分別以正區和負區的比重，和正區和負區涵蓋開放讀碼區的數量，更進一 步的由不同角度的特異性辨識達到越佳引子的挑選，以下詳述兩方針的辨識法。

4-3.2.1 方針一：正區和負區的比重

正區和負區集合的序列建構成字尾樹後，可容易的計算出每個候選引子在正 區和負區的出現數，由此資訊便可得知此候選引子的特異性。圖 4.7 為計算特異 性範例；圖 4.8 為候選引子在字尾樹涵蓋到的正區和負區示意圖。

候選引子 GT 涵蓋 3 條 ORF, 正區出現 2 次, 負區出現 1 次

5’ G A T A C A G T C C C A 3’

5’ G T C C T G G C A C A T 3’

5’ A C C G T C A T A C C G 3’

正區(＋) 負區(－)

忽略區

PN = 2, NN = 1 TN = PN + NN = 2 + 1 = 3

後選引子 GT 的特異性值 = 2/3 = 0.67

圖 4.7：計算方針一特異性值範例

特異性高 特異性低

－ ＋ －

＋ －

圖 4.8：候選引子在字尾樹涵蓋到的正區和負區示意圖

因考慮到引子的挑選準則和實驗室操作的設定，所以在方針一中加入區域權 重的觀念，當挑選候選引子時，考慮引子涵蓋到負區的因素，足以影響挑選的結 果，可將負區的權重加強，使得特異性更加嚴謹；反之，若較忽略引子涵蓋到負 區的因素，則可將負區的權重調降，使得較為突顯正區涵蓋開放讀碼區數量的因 素，又不失引子的特異性。

4-3.2.2 方針二：以涵蓋正區數為主要挑選

一個涵蓋越多開放讀碼區的引子，越利於實驗室中的操作，並且可提高實驗 的成功率和效率，若任兩個候選引子有相同的正區涵蓋開放讀碼區數量，則挑選 負區涵蓋開放讀碼區數量較少的引子，並將另一候選引子排除候選引子列；以此 方針二挑選候選引子，可避免以方針一挑選卻未挑選到正區涵蓋開放讀碼區數量 高的越佳引子。圖 4.9 為兩引子涵蓋相同正區數，不同負區涵蓋數，以方針二挑 選涵蓋正區數高者為候選引子。

5’ G A T A C A G T C C C A 3’

5’ G T C C T G G C A C A T 3’

正區(+) 負區(-)

忽略區

5’ A C C G T C A T A C C G 3’

引子 GT 涵蓋 3 條 ORF, 正區出現 2 次, 負區出現 1 次 引子 GG 涵蓋 3 條 ORF, 正區出現 3 次, 負區出現 0 次

候選引子：GG

圖 4.9：方針二取涵蓋正區多數 ORF 為候選引子

第五章 執行效能

5-1 實驗數據與硬體規格

本論文研究標的與文獻[24]相似，故以文獻[24]所以將其挑選特定的基因組 序列當成本論文實驗主要的數據，且以文獻[16]的 FindGDPs 演算法跑相同數 據，而這兩篇文獻執行出來的結果即為本論文的對照組，便可明顯的看出本論文 的挑選引子演算法改善效果。

本論文的實驗數據為五種不同的微生物基因組，分別為 Bradyrhizobium japonicum、Escherichia coli K12、Streptococcus pneumoniae TIGR4、

Haemophilus influenzae 和 Mycoplasma genitalium，表一為其微生物基因組內 含開放讀碼區序列的數量對照表，這五組微生物內含開放讀碼區序列的數量均不 盡相同，且差距甚大，而以大豆慢生根瘤菌(Bradyrhizobium japonicum)含有 8317 條開放讀碼區序列，其擁有多於一般微生物的開放讀碼區序列數量，故以 此微生物作為 S-Primer 是否能適用於眾多微生物的標的；且挑選第二列的微生 物作為與其他兩文獻實驗數據的比較來源，大腸桿菌(Escherichia coli K12) 含有 4254 條開放讀碼區序列，在微生物中算是相當龐大，若以此資料進行實驗 後的效能良好，相對的在實際生物實驗的應用上應是有貢獻的。本論文的實驗資 料由 NCBI 網站的資料庫中取得微生物基因組。

微生物名稱 開放讀碼區序列數量

Bradyrhizobium japonicum 8317 Escherichia coli K12 4254 Streptococcus pneumoniae TIGR4 2094 Haemophilus influenzae 1657 Mycoplasma genitalium 484

表一：微生物基因組內含開放讀碼區序列的數量對照表

本論文實驗使用的硬體配備。

硬體名稱 規格

中央處理器 CPU AMD 2800+

記憶體 Main Memory 1GB 硬碟 Hard Disk 40G

作業系統 Operation System Windows XP 表二：硬體規格表

5-2 實驗結果與評估

本論文由基本的引子挑選到實驗室使用的引子，所以將實驗循序漸進的分成 兩大項進行，(1)不考慮任何挑選引子準則、(2)考慮特異性辨識對引子選取的影 響，並考慮引子挑選準則，以一項實驗結果證明 S-Primer 可適用於多數的微生 物，並以兩項實驗的結果均和文獻[16、24]的實驗結果互相的比較與評估。

5-2.1 實驗(1)：S-Primer 適用性

以內含 8317 條開放讀碼區序列的大豆慢生根瘤菌，作為引子挑選演算法的 標的微生物，將其開放讀碼區序列全作逆互補序列動作，並取 30%的 3＇端正區 序列長度和 30%負區序列長度，設定 20%的 5＇負區忽略區，且以方針一和方針 二作為挑選候選引子的特異性辨識，為了與文獻[16、24]和隨機引子做比較，故 取引子長度為 6~8；實驗(1)結果如表三、表四所示。

表三為方針一：引子的特異性均需大於等於 0.5，且設定負區序列的權重為 1，作為挑選引子演算法的特異性辨識。

實驗(1)以方針一挑選引子數據

引子長度 引子數量 執行時間(秒)

6 80 106

7 181 114

8 378 127

表三：實驗(1)以方針一挑選引子 表四為方針二作為引子挑選演算法的特異性辨識。

實驗(1)以方針二挑選引子數據

引子長度 引子數量 執行時間(秒)

6 67 58

7 176 41

8 440 30

表四：實驗(1)以方針二挑選引子

實驗(1)在處理龐大數量的開放讀碼區序列資料，仍有不錯的執行時間，並 挑選出少量且符合特異性的引子；方針二平均挑選出的引子集合較少，執行時間 也較方針一快速許多，相對方針一在引子數量和執行時間上都略遜一籌，但以引 子長度 8 的引子挑選，方針一挑選出較少的引子集合，且優於方針二所挑選出的 引子集合數量，故可推出，若引子長度持續增長，方針一可能會挑選出相較方針 二的最小引子集合；方針二均有不錯的引子數量和執行時間，所以方針一和方針 一互有其優缺之處。然而實驗(1)也說明了，S-Primer 可適用於多數微生物開放

5-2.2 實驗(2)：不考慮任何挑選引子準則

我們將大腸桿菌內含的開放讀碼區序列全作逆互補序列動作，並取 3＇端正 區序列 30%的長度建構成字尾樹，且利用引子挑選演算法挑選出引子的數據與其 他兩文獻結果比較，如下表五。

引子長度 6 7 8

本論文引子數 58 124 245

文獻[24]引子數 55 122 250 FindGDPs 引子數 59 132 N/A

表五：實驗(2)結果比較

實驗(2)引子挑選個數與文獻[24]相近，僅在引子長度為 8 時優於文獻[24]

的時驗結果，所以可見當引子長度加長時，本論文可取得較少的引子數量，更優 於無法取得引子長度為 8 的文獻[16]的實驗結果。

5-2.3 實驗(3)：考慮特異性辨識對引子選取的影響

在 4-3.2 中提到引子的特異性的計算方法，所以我們將大腸桿菌內含的開放 讀碼區序列全做逆互補序列動作，取 5＇端忽略區序列 20%長度，並取 3＇端正 區和負區序列 30%的長度建構成字尾樹，以方針一和方針二作為挑選候選引子的 特異性辨識，方針一：設定每個引子的特異性需大於等於 0.5，負區序列權重為 1，且利用引子挑選演算法挑選出引子的數據與其他兩文獻結果比較，表六~表九 為各演算法挑選引子個數的結果，在不同長度的引子挑選下執行的時間和得到的 引子數量。

本論文實驗(3)以方針一挑選引子數據

引子長度 引子數量 最長引子長度 執行時間(秒)

6 85 11 40

7 183 11 45

8 379 22 58

表六：本論文實驗(3)以方針一挑選引子結果

本論文實驗(3)以方針二挑選引子數據

引子長度 引子數量 最長引子長度 執行時間(秒)

6 84 6 16

7 219 7 12

8 389 8 12

表七：本論文實驗(3)以方針二挑選引子結果

文獻[24]實驗(3)數據

引子長度 引子數量 最長引子長度 執行時間(秒)

6 61 11 43

7 133 11 43

8 254 19 54

表八：文獻[24]實驗(3)結果

FindGDPs 實驗(3)數據

引子長度 引子數量 特異性＜0.5 數量 執行時間(秒)

6 60 29 36

7 128 38 78

8 258 75 513

表九：FindGDPs 實驗(3)結果

經由表六、表七和表八、表九的互相比較，可以發現本論文在實驗(3)得到 的引子數量多於其他兩文獻的實驗結果，其結果不如人意，但以特異性角度看來 本論文挑選引子特異性的判別較為嚴謹，故造成挑選的引子數量較多。又以執行 時間來看，本論文的引子挑選演算法遠遠優於其他兩文獻的執行速度，所以整體 來看本論文實驗(3)的結果略微優於其他兩文獻的實驗結果。

5-3 效能分析

在 5-2 節中我們由淺入深的做了三項實驗，由實驗(1)證實了本論文的引子 挑選演算法，可適用於廣泛微生物的開放讀碼區序列，也比較解相同問題文獻的 實驗結果，在實驗(2)和實驗(3)讓我們證實在相同的微生物的情況下，越複雜的 條件依序的加入後，本論文演算法執行出的效能是良好的，以整體而言，從生物 實驗角度來看，本論文的引子挑選演算法在挑選引子集合和執行時間和空間上都 相當有貢獻。

第六章 結論

6-1 研究成果

字尾樹演算法至今雖還未廣泛運用在處理生物問題方面，然而本論文成功的 運用字尾樹在尋找共同子字串的特性和其線性時間和空間，且運用在 DNA 微陣列 的待測物標記上，再經由引子挑選演算法，所挑選出來的引子正好可以降低檢測 時的錯誤雜交發生；本論文針對生物意義方面加入了特異性辨識：方針一和方針 二的辨識方法，更增進引子符合 DNA 微陣列檢測的實務需求，又經實驗證實引子 挑選演算法，在執行時間和正確引子挑選的數量上均優於以往，又因字尾樹演算 法的改進，本論文的空間約為 24n~30n(n 為每個字元所佔位元)，而文獻[24]則 佔空間將近 150n 以上，明顯的，本論文所使用的演算法大幅的降低程式所佔的 記憶體空間。

6-2 未來發展方向

雖然本論文的演算法已大幅降低程式所佔據的記憶體空間，引用的

Gusfield 團隊改進的 Ukkonen 字尾樹演算法，但以文獻[25][27]來看，還有更 佳的演算法可將空間降至 10n 左右；因本論文尚未完全參考到所有的引子挑選準 則，故經引子挑選出來的引子是否真正符合實驗室實務上的標準，然而空間和確 切的生物意義還待未來繼續的研究。

參考文獻

[1] R.B. Wallace, J. Shaffer, R.F. Murphy, J. Bonner, T. Hirose, K., (1979). Itakura, Hybridization of synthetic oligodeoxyribonucleotides to phi chi 174 DNA: the effect of single base pair mismatch. Nucleic Acids Research, 6, 3543-3557 [2] K.J. Breslauer, R. Frank, H. Blocker, L.A. Marky, (1986). Predicting DNA

duplex stability from the base sequence. Proc Natl Acad Sci USA, 83, 3746-3750 [3] M.A. Innis, D.H. Gelfand, J.J. Sninsky, T.J. White, (1990). PCR Protocols : A

Guide to Methods and Ap-plications. Academic Press, New York [4] T. Lowe, J. Sharefkin, S.Q. Yang, C.W. Dieffenbach, (1990). A computer

program for selection of oli-gonucleotide primers for polymerase chain reac-tions. Nucleic Acids Research, 18, 1757-1761

[5] S. Kwok, D.E. Kellogg, N. McKinney, D. Spasic, L. Goda, C. Levenson, J.

Sninsky, (1990). Effects of primer-template mismatches on the polymerase chain reaction: Human immunodeficiency virus 1 model studies. Nucleic Acids

Research, 18, 999-1005

[6] P. O. Brown et al. (1995). Quantitative Monitoring of Gene Expression Pattern with a Complementary DNA Microarray. Science, 270, 467-470

[7] Gusfield Dan, (1997). Algorithms on strings, trees, and sequences : computer science and computational biology. Cambridge [England] New York Cambridge

University

[8] P. Li, K.C. Kupfer, C.J. Davies, D. Burbee, L.A. Ev-ans, H.R. Garner, (1997).

PRIMO: A primer design program that applies base quality statistics for auto-mated large-scale DNA sequencing. Genomics, 40, 476-485 [9] J. Jr. SantaLucia, (1998). A unified view of polymer, dumbbell, and

oligonucleotide DNA nearest-neighbor thermodynamics. Proc Natl Acad Sci

USA, 95, 1460-1465

[10] S. Haas, M. Vingron, A. Poustka, S. Wiemann, (1998). Primer design for large scale sequencing. Nucleic Acids Research, 26, 3006-3012

[11] A.M. Talaat, P. Hunter and S.A. Johnson, (2000). Genome-directed primers for selective labeling of bacterial transcripts for DNA microarray analysis. Nat.

Biotechnol., 18, 679-682

[12] T. Kampke, M. Kieninger, M. Mecklenburg, (2000) . Efficient primer design algorithms. Bioinformatics, 17, 214-225

[13] N. von Ahsen, C.T. Wittwer, E. Schutz, (2001) . Oli-gonucleotide melting temperature under PCR conditions: nearest-neighbor corrections for Mg2+, de-oxynucleotide triphosphate, and dimethyl sulfoxide concentrations with

comparison to alternative em-pirical formulas. Clinical Chemistry, 47, 1956-1961

[14] Jan-Jaap Wesselink and so on, (2003). Determining a unique defining

DNAsequence for yeast species using hashing techniques. Bioinformatics, 18, 1004-1010

[15] Bowtell, David, Sambrook, Joseph, DNA microarrays : a molecular cloning manual. Cold Spring Harbor Laboratory

[16] Robert J. Blick, Andrew T. Revel and Eric J. Hansen, (2003). FindGDPs:

identification of primers for labeling microbial transcriptomes for DNA microarray analysis. Bioinformatics, 19 ,1718-1719

[17] Hsieh,M.-H., Hsu,W.C., Chiu,S.K. and Tzeng,C.M., (2003). An efficient algorithm for minimal primer set selection. Bioinformatics, 19, 285-286 [18] S.H. Chen, C.Y. Lin, C.S. Cho, C.Z. Lo, C.A. Hsiung, (2003). Primer Design

Assistant (PDA): a web-based primer design tool. Nucleic Acids Research, 31, 3751-3754

[19] Steen Knudsen, (2004). Guide to analysis of DNA microarray data. Hoboken, N.J.Wiley-Lissc

[20] E.Barry Simpson, Shannon L. Ross, Sarah E.Marchetti, and John C.Kennell, (2004). Relaxed Primer Specificity Associated with Reverse Transcriptases Encoded by the pFOXC Retroplasmids of Fusarium oxysporum. Eukaryotic Cell, 3, 1589-1600

[21] Jiren Wang, Kuo-Bin Li and Wing-Kin Sung, (2004). G-PRIMER: greedy algorithm for selecting minimal primer set. Bioinformatics, 20, 2473-2475 [22] J.S. Wu, C. Lee, C.C. Wu, Y.L. Shiue, (2004). Primer design using genetic

algorithm. Bioinformatics, 20, 1710-171

[23] 陳容靜、潘人豪、陳永福、林一郎、詹永寬，(2005). 蛋白質功能區域研究 及其多重聚合酶連鎖反應引子設計，國際醫學資訊研討會

[24] 侯玉松、陳華, (2005). A study on Specific Primer selection Algorithms using Suffix Trees

[25] Stefan Kurtz, (1999). Reducing the Space Requirement of Suffix Trees. Softw.

Pract. Exper. 29(13), 1149–1171

[26] 許維中、江彥逸, Evaluation and analysis of cross hybridization effects of probes on cDNA Microarrays

[27] E. McGreight, (1976). A Space-Economical Suffix Tree Construction Algorithm.

Jouraal of the Amociation for Computing Machinery, 23(2), 262-272

Communications of the ACM, 20, 762-72

[29] D. Knuth, J. Morris and V. Pratt, (1997). Fast Pattern Matching Algorithm in