實用的演化分析工具：ape package of R

R 套件 ape 是筆者最近發現的一個非常好用的工具，它可以幫助你從事親緣關係 (phylogenetics) 和演化的分析研究。ape package 自從在 2004 年發表在 Bioinformatics 後，至今已有超過七百次的引用數，可以想見它的應用之廣泛與實用。

標題：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，以完成自己的任務。