為什麼我們需要在某個物種當中搜尋一段DNA序列呢？

前幾天我拿了些例子說明了怎麼使用一些簡單的步驟，來檢視我們拿到的這個樣品裡頭有哪些細菌種類，這算是這個問題的一種應用。在生物資訊領域，為了回答各種不同的問題，我們常常都得透過序列比對來達成這件事情。比如說，當我確信某些序列就是出自某個物種了，但是我想知道這些序列會比對到該物種基因體上的哪些位置，可能是RNA-Seq產出的數據，所以想知道基因的表現量；也可能是ChIP-Seq產出的數據，想知道目前研究的蛋白質會跟哪些DNA片段作用。所以，當我們手邊有序列資料時，常常都是得對各種參考序列做比對的。今天想分享的，不是我前幾天寫的那些使用BLAST之類的本地端序列比對工具來做這件事，今天我要分享的是怎麼使用線上的比對工具，但是是透過Julia這個語言工具來做。

又是BioServices.jl

今天我打算使用BioServices.jl裡頭的GGGenome這個module。這個GGGenome module其實是採用了日本的一個研究中心Database Center for Life Science（簡稱DBCLS）開發的一個線上超快DNA序列搜索工具GGGenome所提供的REST API實作而成的。那我們馬上來看看在Julia底下要怎麼用這個工具呢？

• 基本用法

using BioServices.GGGenome
gggsearch(query::AbstractString; 
        db="hg19", k=0, strand=nothing,
        format="html", timeout=5,
        output=nothing, show_url=false)

這裡頭，query::AbstractString指的是一段核苷酸序列（由ACGT組成，大小寫都無所謂）。db指的是要搜尋的目標資料庫，目前DBCLS官方提供的可供搜尋資料庫都列在這裡了（也可透過gggdbs()這個函數列出）。strand指的是我們要搜尋的序列有沒有分正股跟反股（還記得以前生物課提到DNA時都會講到雙股螺旋吧？）。k指的是在這次的搜尋裡頭，我們可以允許多少核苷酸比對不上（mismatch/gaps），超過這個數字的一律算是比對不上。format則是指比對結果的輸出格式，目前BioServices.GGGenome允許了html|txt|csv|bed|gff|json這幾種格式，如果不指定的話，就會以html格式輸出。output這個選項也很重要，它指定了序列比對結果的輸出形式，預設是返回HTTP.Messages.Response類別的物件，但也可轉換成文字。由於我們給的query序列可能會比對到目標資料庫中的很多地方，如果我們設定成toString，它會全都輸出；但設定成extractTopHit，則會輸出比對結果最好的那一個，不過目前這個選項只有在format設定成txt時才有效。最後這個show_url::Bool則是指要不要顯示這個函數會生成怎樣的網頁連結來讓GGGenome接收我們的request之後返回我們要的結果。

• 幾個簡單的範例

  • 我想在hg19這個人類基因體的資料庫當中搜尋TTCATTGACAACATT這段序列，不管是比對到正股還是反股，同時我只想看可以完全比對上的結果，結果以BED格式返回：

  res = gggsearch("TTCATTGACAACATT", format="bed", output="toString")
  print(res)
  

  • 在小鼠基因體mm10這個資料庫裡搜尋TTCATTGACAACATTGCGT這段序列，允許2個mismatchs或gap只看可以比對上正股的結果，結果以txt返回：

  res = gggsearch("TTCATTGACAACATTGCGT", db="mm10", k=2, strand="+", format="txt", output="toString")
  print(res)
  

希望這個工具對您有幫助。