標題：APE:Analyses of Phylogenetics and Evolution in R language
就以測量兩條 DNA 序列的距離來說好了，ape 能夠彈指間幫你計算完成，真是大快人心。可以泡杯咖啡，想想晚上要看什麼戲劇了。
由於筆者目前有比對 DNA 序列的需要，所以底下會簡單介紹如何利用 ape package 計算多條序列之間的相似程度 (similarities of DNA sequences)。
輸入資料 (讀取檔案)
比對序列前，至少要先有 input 吧。ape 擁有讀取多種資料格式的函式，並且轉化成自己專屬的資料型態。目前它支援的檔案格式有 "interleaved", "sequential", "clustal", 和 "fasta"。我想後兩者應該是很多生資學生熟悉的吧。
您可以在 R console 裡輸入
?read.dna 了解更多資訊和範例說明。為了方便解釋，筆者僅以 fasta format 為例：### ... 建立一個擁有三條 DNA 序列的 fasta 檔案， ### ... 取名為 exdna.txt cat("> No305", "NTTCGAAAAACACACCCACTACTAAAANTTATCAGTCACT", "> No304", "ATTCGAAAAACACACCCACTACTAAAAATTATCAACCACT", "> No306", "ATTCGAAAAACACACCCACTACTAAAAATTATCAATCACT", file = "exdna.txt", sep = "\n") ### ... 讀取 exdna.txt ex.dna <- read.dna("exdna.txt", format = "fasta")
計算距離
我們使用
dist.dna() 計算多個序列間的距離。 您可以鍵入 ? dist.dna 或者到此網頁看更多資訊。dist.dna() 支援的 model 非常多，有 "raw", "N", "TS", "TV", "JC69", "K80" (the default), "F81", "K81", "F84", "BH87", "T92", "TN93", "GG95", "logdet", "paralin", "indel" 以及 "indelblock"。不過我想最為人知的還是 1980 年 Kimura 所提出的 “Kimura's 2-parameters distance” model。在沒有特別指定 model 之下， K80 是
dist.dna() 的預設 model。### ... 計算距離 dist.dna(ex.dna) ### 結果 No305 No304 No304 0.05561282 No306 0.02703361 0.02703361
等等
筆者調整了範例中的 DNA 序列，修修改改後意外地發現若提供了太短的序列的話，會出現「錯誤在obj[i, ] <- sequ : 被替換的項目不是替換值長度的倍數」的警示訊息，然後就卡關了。原因是什麼呢？原來是
read.dna() 原始碼之中有最短長度的限制。筆者手邊的序列長度皆不到 10 bp，而 read.dna() 卻設有長度限制。要如何破解它呢？破解方法
先在 R console 裡輸入
read.dna 印出原始的函式，將程式碼複製，接著打開記事本或任何上手的文字編輯器，貼上 (不要保留 <environment: namespace:ape>)。找到這一行，
pat.base <- "[-AaCcGgTtUuMmRrWwSsYyKkVvHhDdBbNn?]{10}"
把後面的 10 換成你想要的任意長度，比如改成 7 好了：
pat.base <- "[-AaCcGgTtUuMmRrWwSsYyKkVvHhDdBbNn?]{7}"
為了避免破壞原來的函式，我們將修改過後的函式重新定義為自己的函式，比方說， 改名為
seyna.read.dna() 好了。在程式碼開頭處，找到：function (file, format = "interleaved", skip = 0, nlines = 0, comment.char = "#", seq.names = NULL, as.character = FALSE, as.matrix = NULL) {
我們只需要放上新定義的函式名稱，其餘維持不變就好了 (修改過的地方為紅色處)
seyna.read.dna <- function (file, format = "interleaved", skip = 0, nlines = 0, comment.char = "#", seq.names = NULL, as.character = FALSE, as.matrix = NULL) {
接著把筆記本上完整的程式碼複製貼上到 R console 裡，執行，讓 R 儲存此函式 。然後您就可以執行變更過長度限制的
seyna.read.dna()，讀取更短的 DNA 序列以進行比對囉# ... 用自行定義的函式讀取 DNA 序列檔 ex.dna <- seyna.read.dna("exdna.txt", format = "fasta") dist.dna(ex.dna)
更多
R 套件 Biostrings 是個強大的字串處理工具，也能做到類似的事情。只是其建立相似程度的模型並不包含基於 K80 等演化模型，多是 BLOSUM、PAM 或參數化的 nucleotideSubstitutionMatrix。然而使用者也可以自行定義 substitution matrix，以完成自己的任務。
