實用的 ncRNA 資料庫

自己準備 non-coding RNA (ncRNA) 的清單未嘗不可。讓我們來想像一下準備的步驟。一，從表現量的資料中過濾出可靠的轉錄物，二，找出這些轉錄物的座標，三，比對已知的蛋白質的對應的位置，四，分離出 protein-coding transcript 與 non-coding transcript。然而當中最麻煩的步驟可能就是找出所有已知蛋白質的資訊，因為這需要蒐集各大資料庫的蛋白質註記，還有令人白髮的諸多資料前置處理的程序。

建立自己的清單還是相當重要，不過若是資料來源相當的可靠(自己動手做也比不上人家的話)，或者只是想初步檢驗問題的潛力，那麼使用別人已建立好的資源可能會是幫助自己研究過程更有效率的辦法之一。

在 2008 年，John S.Mattick 的實驗室建立了 NRED 資料庫，全名是 Noncoding RNA Expression Database。這個資料庫厲害之處在於不僅展示自己實驗室得到的數據，還搜刮了 GNF 與 Allen Brain Atlas 兩個大型資料庫的資訊，前者擁有 44,775 筆人類、36,182 筆老鼠的 transcript 的資料，後者則藏有成鼠大腦裡大約兩萬筆 transcript 表現量的資料。看來已是相當豐富。

他們將這些蒐集而來的轉錄物和底下三個來源的蛋白質資料做比對： Refseq [1] , MGC [2] , UCSC known genes [3] 。除此之外，他們也利用 CRITICA 演算法 [4] 來預測轉錄物有多少 protein-coding 的潛在能力。剩下的轉錄物，若沒有比對到已知的 ncRNA 就會被歸類為 'Unknown'，比對到已知的 ncRNA 的話則歸類為 'noncoding'。

NRED 擁有很貼心的介面，設計了許多過濾器(filter)。表現量： P-values, B-statistics, q-values。 fold change：M-values。表現量的強度： A-values, Affymetrix present/absent calls。另外還提供了許多有助於研究的資訊，如 transcript 有無 spliced，有沒有 PhastCons 的 conservation 分數，有沒有 RNAz 預測的二級結構等等。也提供了關於 sense, cis-antisense 等資訊。

決定要拿這個資料庫裡的資料來玩一陣子了。


[1] Pruitt, K.D., Tatusova, T. and Maglott, D.R. (2005) NCBI Reference Sequence (RefSeq): a curated non-redundant sequence database of genomes, transcripts and proteins. Nucleic Acids Res, 33, D501-504.
[2] Gerhard, D.S., Wagner, L., Feingold, E.A., Shenmen, C.M., Grouse, L.H., Schuler, G., Klein, S.L., Old, S., Rasooly, R., Good, P. et al. (2004) The status, quality, and expansion of the NIH full-length cDNA project: the Mammalian Gene Collection (MGC). Genome Res, 14, 2121-2127.
[3] Hsu, F., Kent, W.J., Clawson, H., Kuhn, R.M., Diekhans, M. and Haussler, D. (2006) The UCSC Known Genes. Bioinformatics, 22, 1036-1046.
[4] Badger, J.H. and Olsen, G.J. (1999) CRITICA: coding region identification tool invoking comparative analysis. Mol Biol Evol, 16, 512-524.

閱讀更多 »

淺論分辨蛋白質編碼基因和 ncRNA 的困難

人類擁有超過三億個鹼基對，然而卻只有 2.3% 的部份能夠轉錄出 mRNA。這讓人們對於那些最終未能轉成蛋白質的 DNA 區域感到好奇，它們的存在有什麼作用呢？甚至於有些科學家還以 junk DNA 形容這些區域。直到 2005 一篇發表在 Science 的研究，指出了其實細胞內有大量不在預期中的轉錄產物，而這些轉錄產物絕大多數似乎都不會成為蛋白質，了解到這項發現的重要性，於是人們開始以「暗物質」(Dark matter ) 等引人遐想的名稱來形容這些物質 [1]。

毫無意外地，自那時候起便開始有了許多針對這些物質的研究。也因此出現了許多新的類別，比如： microRNA, PIWI-associated RNAs, endogenous small interfering RNAs ...等。而之後科學家在透過大規模地實驗發現長度較長的 (超過 200 nt) noncoding RNA (ncRNA) 似乎是這些「暗物質」的主要成份，研究結果也暗示這些長的轉錄產物可能是之後短的、具功能的 ncRNAs 的 precursor [2]。

然而這類的研究也並非總是一帆風順。有許多標記上的問題需要面對。怎麼說呢？在當時，老鼠的基因體中有 33% 的轉錄產物被標記為 noncoding，然而換成了人類的基因體，卻只有 7% 的轉錄產物被標記為 ncRNAs。 26% 是很大的差異，而這樣的差異也凸顯了將長 ncRNA 與 mRNA 清楚分類的困難度。

底下，筆者就會針對六種不同的策略，剖析科學界在區分 protein-coding 與 noncoding RNAs 上所面臨的混沌與挑戰。


策略一： Open reading frame 長度
有鑑於超過 95% 的蛋白質擁有超過 100 aa 的長度，所以FANTOM consortium 原本使用 300 nt，也就是 100 codons 的長度作為門檻，但誰都知道這樣是很粗糙的辦法。一些已知的具功能的 ncRNA 如 H19, Xist, Mirg, Gtl2 都會因此被剔除掉。另一方面，也有大約近 3700 個蛋白質會小於這麼門檻。


策略二： ORF conservation
科學界普遍存在的一種假設是，若 ORF 缺乏 conseravtion，則其「沒有功能」的可能性很高。一些常用的計算 ORF conservation 的工具有 BLASTX. Pfam, rsCDS, SUPERFAMILY。其它還有一些加入了跨物種資訊的計算工具，包括 CSTminer, CRITICA 等。這種方法的限制在於現有蛋白質資料庫的「完備程度」和「準確程度」。舉例來說， Xist 在 1992 年被確認為具功能性的 ncRNA 之前，被錯誤地歸類在蛋白質當中，時間長達十五年之久。另外，以 conservation 為依據也可能造成誤判，理由是一些 ncRNAs 可能是從 protein-coding 基因演化而來，因此可能留有 mRNA 的遺跡。比如說， Xist RNA 基因便是一個 protein-coding 基因 pseudogenization 後的結果 [3]。


策略三：Structural approaches
利用 conserved predicted RNA 的二級結構作為判斷 ncRNA 的依據。常用的工具有 QRNA, RNAz, EvoFOLD 等。然而這種方法很容易造成許多誤判，而且，其實 mRNA 裡頭也擁有許多 conserved 二極結構(尤其是 3'UTR)。


策略四：實驗手段
實驗室利用 In vitro translation assays 來鑑定一個 ORF 有無轉譯成為一個蛋白質。另一種方法是評估一個轉錄產物是否與 polysomes 緊密連結(雖然這種方法仍不是很可靠)。實驗手段很可能是所有手段裡最為可靠的(直覺上)，然而在解讀這些實驗結果時我們仍須特別小心。原因是有些研究結果顯示， in vitro 中的肯定結果，在 in vivo 的情況下卻未能偵測到應該有的蛋白質。而另一方面，in vitro 實驗中的負面結果，我們也不能信心滿滿地百分之百地保證 in vivo 也能重複見到一樣的事情。


策略五：Artifact 過濾
舉例來說，FANTOM3 consortium 利用了很嚴謹的條件，包括 5'end 與 3'end 的 CAGE signal strength 和需吻合 polyadenylation 等，將原來的三萬多條 ncRNA 過濾到只剩三千六百多條 [4]。還有方法將那些會與另一股會產出 mRNA 重疊到的 ncRNA 通通丟棄。可想而知，若按照後者辦法，現在很熱門的 anti-sense transcript 等研究通通都無法做了。


策略六：綜合辦法
在 2006, 2007 年有兩個方法不約而同的使用了 SVM，透過學習 peptide length, amino acid composition, protein homologs, secondary structure, protein alignment information 等特徵，來預測 protein-coding 的機率值 [5,6]。這是一個將生物問題成功重塑為機器學習問題的好故事，也都刊上了很好的期刊。然而，光依賴(黑盒子)機器提供的預測數據，我們是無法對於核心問題 ─ protein-coding gene 與 ncRNA 之不同，有更清楚深入的了解的。


結論：
近來人們在多細胞和單細胞生物體內觀察到一種神奇的現象，那就是有些轉錄產物不僅可以以蛋白質的角色進行任務，也可以以 RNA 的角色執行工作。例子包括人類的 steroid receptor activator，果蠅的 Oskar mRNA，或者是 E.Coli 的 SgrS RNA。這些觀察都在挑戰著科學家之前對於轉錄產物的認知。RNA 事實上是可以有雙重角色的，一個是當它為蛋白質時擁有的，另一個是當它以 RNA 的角色出現時擁有的。這樣的結果衝擊人們熟悉的二分法，理論開始有了修補的需要。毫無疑問地，未來關於 ncRNA 的研究數量會更加速地成長，刺激人們的想像，開拓更多知識的疆土。眼前看來有無窮盡的寶藏即將開展。然而讓我們在此稍作歇息，拍拍肩上的灰塵，深呼吸，然後問問自己，你想出一個好問題了嗎？


[1] In the forests of RNA dark matter. Science, 2005
[2] RNA maps reveal new RNA classes and a possible function for pervasive transcription. Science, 2007
[3] The Xist RNA gene evolved in eutherians by pseudogenization of a protein-coding gene.Science, 2006
[4] The transcriptional landscape of the mammalian genome. Science, 2005
[5] Distinguishing protein-coding from non-coding RNAs through support vector machines. PLoS Genet, 2006
[6] CPC: Assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res, 2007

photo credit: NASA's Marshall Space Flight Center via photo pin cc

閱讀更多 »

處理 R 的 sparse matrix 問題

遇到 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 計算相似度

今天認識到了 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

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

近來有個問題在腦中逐漸成形，出現了將 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

筆者會注意到 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"

閱讀更多 »

另一則實驗：G-quadruplex 與轉錄因子 SP1 之結合

《核酸研究》 (Nucleic Acids Research) 又刊出了一個關於 G-quadruplex 的文章。標題是「透過體外實驗，非典型的 DNA 結構是轉錄因子 SP1 的結合特徵序列」(A non-canonical DNA structure is a binding motif for the transcription factor SP1 in vitro)   

他們的研究與前人的研究有很大的差別在於，他們所使用的G-quadruplex 候選序列是只會形成兩個平面 (tetrad) 的，而不像之前人們多是採用會形成三個以上平面的候選序列。然而從他們的研究結果中他們認為，這種只會形成兩個平面的 G-quadruplex 在真實的細胞環境下有可能存在。它們被討論的太少了，不應該被如此忽視。   

已有前人展示過鋅手指(zinc-finger)蛋白質和 G-quadruplex 結合的能力。好玩的是，SP1 也是一個鋅手指蛋白質。而透過 DNA 突變的分析實驗，他們展示了 SP1 蛋白質除了辨識雙股 DNA 上的結合序列之外，也會將 DNA 是形成特殊結構的 G-quadruplex 視為結合的另一個對象。   SP1 蛋白質其中一個最短的辨識序列是 5'-GGGCGG-3'，擁有相當多的鳥嘌呤(guanine)。他們的研究數據顯示 SP1 的結合位置有高達 77-87% 是和可能形成 G-quadruplex 的位置重疊的。   

更有甚者，之前科學家已發現 G-quadruplex 能夠限制 CpG 核苷酸的甲基化。Sp1 可能會透過和 G-quadruplex 結合，避免 DNA 序列後續的甲基化，達到調控基因的效果。   

「我們的研究展現了人們至今對於 SP1 蛋白質對於辨識結合上仍不了解的特質，未來可能被證明對於調控基因的表現至關重要」Balasubramanian 教授說道。  

原始文章： "A non-canonical DNA structure is a binding motif for the transcription factor SP1 in vitro"
Eun-Ang Raiber, Ramon Kranaster, Enid Lam, Mehran Nikan, and Shankar Balasubramanian Nucl. Acids Res. (2012) 40 (4): 1499-1508.


(本文同步刊載於生資櫥窗 BioWindow)

閱讀更多 »

測序儀大戰：GridION 和 MinION 出陣

如果真如新華網的新聞稿所言 (英研製出"USB"基因組測序儀)，能在串接的方式在 15 分鐘內以僅僅 5000 美元的代價定序完成人類的基因組，那未來定序儀的戰爭可會非常有趣了。儘管先前 Life technologies 已經率先宣成達到以 1000 美元定序人類基因組的目標 ("Ion Torrent claims to be first with $1K genome sequencer")。

但我不覺得 4% 的錯誤率與產品「一次性的特性」僅是「小小的瑕疵」。畢竟若要真能實際運用到學術研究上頭，相較於其他 I、L 公司不到 1% 的錯誤率，這個英國牛津納米孔技術公司的準確度就有點難以令人接受了。

閱讀更多 »

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

遞迴區分 (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 食譜了沒?

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

閱讀更多 »

美女科學家對於人類新生基因起源的觀點

2011 一月，
Aoife McLysaght 在都柏林的一場演講

在最新一期的 PLoS Genetics，愛爾蘭的演化生物學家 Aoife McLysaght 寫下了她對於人類新生基因 (De novo gene) 起源的觀點。(PLoS Genetics, 2011)

相較於其它對於此一問題的答案，比如 exon shuffling, tandem duplication, retrocopying, segmental duplication 以及 genome duplication。Aoife McLysaght 的答案與研究的方式顯然很與眾不同(也很直接)。

她試著從分析研究其他靈長類的非編碼區序列，找出最有可能的新生基因。

除了在科學上的表現，筆者(害羞地)認為 Aoife McLysaght 美的簡直就像是好萊塢的明星。不相信的話請看此影片(youtube)。

閱讀更多 »