取火之路: R

顯示具有 R 標籤的文章。顯示所有文章

處理 R 的 sparse matrix 問題

張貼者： Unknown 2012年4月11日星期三

遇到 sparse matrix 是很棘手的事情。原因是一不小心它會使檔案變得肥大，佔盡了記憶體資源，同時拖垮了 script 的效率。

真的是很棘手的敵人。

幸虧已經有許多人研究、提出了一些解決辦法。比如說 R 的 Matrix 套件。

了解格式

比較推薦使用的是 MatrixMarket 格式。

5  5  8
1  1  1.000e+00
2  2  1.050e+01
3  3  1.500e-02

第一行記錄的分別是 rows, columns, entries (列、行、總數)。接著每一列表示每一個 entry。可加入註解，行首以 % 標示，放在檔案的開頭處。更詳細的說明可以見此網頁。

建立 sparse matrix

可參考底下的範例，試著在 R 重複一遍(鍵入 ?sparseMatrix 查更多說明)。

i <- c(1,3:8); j <- c(2,9,6:10); x <- 7 * (1:7)
(A <- sparseMatrix(i, j, x = x))

輸出結果

8 x 10 sparse Matrix of class "dgCMatrix"
                             
[1,] . 7 . . .  .  .  .  .  .
[2,] . . . . .  .  .  .  .  .
[3,] . . . . .  .  .  . 14  .
[4,] . . . . . 21  .  .  .  .
[5,] . . . . .  . 28  .  .  .
[6,] . . . . .  .  . 35  .  .
[7,] . . . . .  .  .  . 42  .
[8,] . . . . .  .  .  .  . 49

讀、寫外部 sparse matrix

我們可以藉著 read.table() 讀入 MatrixMarket 格式的檔案，再利用 sparseMatrix() 製作 sparse matrix。筆者相信應該還有更好的方法，讀者可以自行嘗試其他指令。底下的範例亦僅供參考。

讀

dd <- read.table("D:/test.mtx", header=FALSE)
mm <- sparseMatrix(i = dd$i, j = dd$j, x = dd$x, dims=c(row_dim,col_dim))

寫
如果你已經擁有了一個 sparse matrix，使用 writeMM() 可以幫您輸出成為外部檔案

data(KNex)
writeMM(KNex$mm, "mmMM.mtx")

這部份可鍵入 ?readMM 看更多說明。

用 Jensen-Shannon divergence 計算相似度

張貼者： Unknown 2012年4月5日星期四

0 意見

今天認識到了 Jensen-Shannon divergence，它似乎是一個常用於計算兩個機率分佈間之相似程度的方法。

where

and

is Kullback–Leibler divergence

根據維基百科的條目，它又被稱為 information radius (IRad) 或 total divergence to the average。儘管是筆者自己的猜測，但我想後者名字的由來應該與 M = 1/2(P+Q) 這個平均值脫不了干係。

2009 年 PNAS 刊了一篇充分運用 Jensen-Shannon divergence 方法以比較多種基因體的文章。對於想運用 Jensen-Shannon divergence 以計算 DNA 序列間相似程度的話，裡頭的公式與步驟應已經夠詳細與清楚了。

標題：Alignment-free genome comparison with featurefrequency profiles (FFP) and optimal resolutions

另外，筆者發現 R 有一個名為 textcat 的套件，裡頭含有可以計算 Jensen-Shannon divergence 的函式 ─ textcat_xdist()。這個函式的優點在於，給定 n-gram profiles，它還可以計算 Kullback-Leibler I-divergence, Kullback-Leibler J-divergence, the sum of the absolute differences in n-gram log frequencies, 以及 cosine dissimilarity 等等。

雖然有現成的 R 套件可以使用，但也不需要高興的太早。依照筆者這幾天的經驗，資料量太大(一千個維度以上) 的情況就開始可能會有「記憶體不足」的錯誤發生，差不多的意思就是這函式報廢了。

所以最好還是從頭了解，弄清楚每個算式，別擔心把手弄髒。如此一來最壞的情況下也有辦法實作演算法，把任務完成。

提醒

若要得到距離，記得將 Jensen-Shannon divergence 開平方根

The Jensen-Shannon divergence is not a distance (as it does not obey the triangle inequality), but its square root is. ─ "What Makes a Query Difficult?" Carmel et al., ACM, 2006.

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

張貼者： Unknown 2012年3月27日星期二

0 意見

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

值得注目的 affinity propagation clustering

張貼者： Unknown 2012年3月26日星期一

0 意見

近來有個問題在腦中逐漸成形，出現了將 DNA 序列做分群的需求，因此著手尋找有用的套件。結果我找到了 APCluster，一個基於 affinity propagation (AP) clustering 的 R 套件。什麼是 affinity propagation clustering？它很強大嗎？查閱了一下原始出處，才發現這個演算法竟然是發表於 2007 年的《Science》。

標題："Clustering by Passing Messages Between Data Points"

《科學》裡的核心概念

在開始談 AP clustering 之前，先聊聊知名且常被使用的 K-centers clustering 演算法。K-centers 演算法會從一些隨機取得的 "exemplars" (被當做 center 的 data points)開始，並且迭代調整更新每個集合以減少整體的 suqared error。這個方法對於初始選擇的 exemplars 很敏感，所以通常需要更換初始選擇並且重跑許多次以獲得好的解答。然而這樣做通常也只有在 clusters 的數目少且恰好初始時選的種子 (seeds) 夠好才會趨近好的解答。

他們的方法中的一個創舉是，「同時考慮全部的資料點為可能的 exemplars」！藉由點與點之間的遞迴地交換資訊 ─ "responsibilities" 與 "availavilities"，讓各個節點能夠知道誰是最適合自己的 exemplar 以及自己有多大程度適合當 exemplar。

真是妙極了。

若想要知道較多實作的細節，不妨查閱原始論文。如果您偏好分析程式碼， AP cluster 也有 Matlab 的實作版本。在程式裡頭的註解檔，有交待了演算法的整體概念，比如：

Each cluster is represented by a data point called a cluster center, and the method searches for clusters so as to maximize a fitness function called net similarity. The method is iterative and stops after maxits iterations.

For N data points, there may be as many as N^2-N pair-wise similarities (note that the similarity of data point i to k need not be equal to the similarity of data point k to i). These may be passed to APCLUSTER in an NxN matrix s, where s(i,k) is the similarity of point i to point k.

APCluster：R 套件

AP clustering 演算法也被實作成了 R 的套件，發表在 2011 年的 Bioinformatics。

標題：APCluster: an R package for affinity propagation clustering.

其實這篇論文才是我的主角。想了想後，發現序列不能以有序數值表示，代表內建的 negDistMat( ) 在分群序列上頭是沒有意義的。若參考作者們在分析蛋白質 coiled coil 所採取的方式，還得詳細了解 Carsten et al. 在 2011 的研究工作才行 (似乎不脫離 BLAST pair-wise alignment 的計算分數)。

程式擁有其他內建的函式，能將資料視覺化

參考文件
APCluster 的參考文件有很多，底下僅列舉一些：

APCluster - An R Package for Affinity Propagation Clustering
http://www.bioinf.jku.at/software/apcluster/
R Graphical Manual : Affinity Propagation Package
http://rgm2.lab.nig.ac.jp/RGM2/func.php?rd_id=apcluster:apcluster-package
apcluster in R CRAN
http://cran.r-project.org/web/packages/apcluster/index.html

研究演化的 R 套件：RPHAST

張貼者： Unknown 2012年3月23日星期五

0 意見

筆者會注意到 PHAST 的原因首先是因為身邊同事使用 PhastCons 的輸出資料進行分析，然後接著因為自己手上新題目的需要而尋找適當的工具，結果找到了 RPHAST ─ 一個 R 版本的 PHAST。

RPHAST 的官方下載網頁：http://compgen.bscb.cornell.edu/rphast/
原始文件：PHAST and RPHAST: phylogenetic analysis with space/time models. Brief Bioinform (2011)
http://bib.oxfordjournals.org/content/12/1/41

沒用過倒算了，一試用才發現這實在是太好用了。與同樣頂頂大名的 PAML 相比，RPHAST 可說像是有視窗介面的 windows 3.1 ，而 PAML 則是醜陋又 bugy 的 DOS。RPHAST 提供了足夠充分的說明文件，而且安裝起來快速容易，真是揪感心啊。

若是需要分析基因體，官網裡的 vignette1.pdf 以內建的番茄第二號染色體為例，示範了讀取 alignments、features，以及用 PhyloFit, PhyloP 與 phastCons 評估 neutral model 的方法。簡單清楚。

僅僅 ~50 行的 R 程式碼，就能產出這張圖

USCS 的 Table Browser 裡也可以找到 phastCons table，只是很可惜地對 human genome 來說，只有提供 44 個物種的比對 (phastCons44wayPrimates)。但其實若只是想要 (ultra-) conserved region 的座標，在這裡拉拉選選就可以得到輸出，省去了自己執行程式的過程，也是很愜意啊。

PHAST，最早是 2002 年由 Cornell University 的助理教授 Adam Siepel 開啟的計畫，到 2006 年都陸續有相關著作，然而卻突然地中斷，直到 2010 年 Melissa J. Hubisz 的加入，在 2010 年和 2011 年大爆發在 Genome Research、 Briefings in Bioinfomatics 等期刊發表著作。若您好奇 PHAST/RPHAST 與其他工具的差別，可以見其官方網頁"Comparison with Other Packages"

利用 partykit 進行遞迴區分 (recursive partitioning)

張貼者： Unknown 2011年11月19日星期六

0 意見

遞迴區分 (Recursive partitioning) 是一種多變量分析的統計方法。遞迴區分，基於一些可二分的 (dichotomous)相依變數，試圖建立一個可以正確分類成員的決策樹。

雖然遞迴區分多用於醫療診斷測試，但其實它的應用程度相當廣泛。

相較於迴歸分析會建立一個公式，讓醫療照護者從中計算一個病患擁有疾病的機率；遞迴區分所建立的規則像是「如果病患身上有 x, y, 或 z, 那他可能會得了 q 疾病」

R 當中的 rpart 套件可以幫你做遞迴區分。 Torsten Hothorn 及 Achim Zeileis 進一步實作了擴充的工具組，稱作 partykit，可以幫助呈現、摘要以及視覺化所得模型的樹狀結構。

底下，我們會使用 HELP 資料，分類成員是否「無家可歸」(homeless)。

ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
library(rpart); library(partykit)
ds$sub = as.factor(ds$substance)
homeless.rpart = rpart(homeless ~ female + i1 + sub + sexrisk + mcs +
  pcs, method="class", data=ds)
plot(homeless.rpart)
text(homeless.rpart)

好極了，有點樣子了 (儘管有點陽春，但沒關係，待會我們會美化它)。我們可以用指令 printcp() 將分類樹的產物顯示出來。

printcp(homeless.rpart)

接著我們要利用 partykit 套件製作漂亮的圖形。注意我們使用了 plot.party() 和 as.party()，強迫 rpart 物件轉成適當的格式。

plot(as.party(homeless.rpart), type="simple")

你也可以將 type="simple" 移掉，看一看終端節點 (terminal node) 的分類情況。

plot(as.party(homeless.rpart))

想要知道更多的資訊和更多迷人的圖片，不妨參考此文件。文件裏頭有取出每個中間節點的 p-value 的範例。

相關資源:
Quick-R: http://www.statmethods.net/advstats/cart.html
消息來源: R-bloggers

你看過 R 食譜了沒?

張貼者： Unknown

0 意見

Yanchang Zhao 於 R-bloggers 上介紹了一個學習 R 的好網站： http://code.ca-net.org/R%20Cookbook。基於不藏私的理念，站主想將此分享給欲學習 R 的人們。

使用 R cookbook，你可以：

藉由 RSQLite、RMySQL、RdbiPgSQL、RODBC
您可以存取資料庫
讀取與寫入資料
日期/時間變數
繪圖
空間資訊 (Spatial data)
除錯 (Debugging)

筆者非常感激裏頭精簡實用的程式碼，往往能幫助筆者很快地製作出渴望的圖表。相信您也能從中獲益許多。

更多相關資源：
RDataMining: http://www.rdatamining.com
Twitter: http://www.twitter.com/RDataMining
Group on Linkedin: http://group2.rdatamining.com

photo credit: Βethan via photo pin cc