Query Linked Data所使用的技術SPARQL:過去歷史和1.1版本

SPARQL的歷史

在RDF於1999年變成標準時,一直到2004年,有很多種查詢用的語言圍繞者商業或是學術需求被開發出來,W3C 於是成立了一個小組RDF Data Access Working Group(RDF-DAWG)來搜集和整理相關的需求,在2004年底左右發表了第一個SPQRQL相關的標準草搞,同時期原本就開始出來用來查詢RDF的工具也陸續支持SPARQL的標準,在2009年W3C RDF-DAWG小組重新改制成SPARQL working group。

在SPARQL 1.0標準裡面,主要有三份文件
1. SPARQL Query Language for RDF
* 描述query的語法
2. SPARQL Protocol for RDF
* 描述成是在執行query時後,要回傳的格式
1. SPARQL Query Results XML format
* 回傳的資料支持xml格式

在最近的SPARQL 1.1版本再新增了8個文件:

1. Overview文件
2. Federate Query如何跨資料檔案的查詢,為了可以支持在distributed 環境下
3. 新增了update的概念,這部分是1.1和1.0最重要的差異,從此可以新增、刪除、取代的功能
4. 關於sevice的功能敘述,主要是為了發展SPARQL client的開發者所寫
5. 支持query結果以JSON格式回傳

screenshot.png
6. 支持query結果以CSV/TSV格式回傳

screenshot.png
7. 支持Graph Store HTTP protocol
8. 支持entailment regimen

PBS(Portable Batch System)系統效調細節筆記:使用叢集係統來做百萬級抗原變異點的預測

之前有測試幾種加速運算的方法,用python裡的線程搭配外掛程式來躲避GIL的效能限制等,但不同種類的架構能應付的問題(高I/O, 高運算…),這邊嘗試使用手邊另一個叢集系統來做計算,利用PBS中分發任務的能力來調用底下的電腦群們來做運算,過程也是一堆的大坑,主要的框架是用shell script寫成的,因為手邊這個叢集電腦python是舊版的,然後還得確定下面的電腦們每個環境都要乾淨,所以就先試者用最簡單的方式來加速計算,最後效調的結果可以在一周內分析完3425328個蛋白序列,雖然還是不夠快。(不確定是否是最佳做法,可能偏“土炮”)

整個過程有幾個重點:
1. 使用shell指令來處理大檔案的技巧
2. PBS系統的理解

第一部分在shell script的撰寫上,相對於高階語言如python,就少掉很多如線城池、共享變數等等用來暫時存放運算資料的方法,這時候就會需要比較複雜的代碼來實現想要的邏輯。技術上可以先把script用最小單位來撰寫,接者在測試成功後,慢慢合成一個大的,因為當牽涉越多層的呼叫使用,需多未知的東西狀態出現,尤其是想要用&的方式來加速(因為目前這叢集系統無法安裝linux parell指令)。

這邊可以堆間使用awk和sed的功能來做不雜一點的計算和資料搬移(他們很輕且通常很快),尤其搭配awk裡面的條件語法和函數,計算速度可以因此稍微提生一點,且能用“一行”來解決在高階語法可能要很多行代碼才搞定的字串處理問題。

第二部分關於PBS語法和功能的理解,雖然每個叢集系統多少會不同(就跟每個linux管理員管理的想法不同一樣,會對自己的系統加入不同的觀念),但指令的說明是差不多的。

三大關於PBS的指令分別為:pbsnodes, qstat, qsub

pbsnodes: 用來看目前叢集系統中,各個節點的資源和任務條波狀況

screenshot.pngqstat:用來看目前任務的運行狀況和單個任務所使用的資源

screenshot.png

qsub:提交任務到PBS系統上

screenshot.png

 

 

還有意外發現使用parent-child的shell process中,當程式的架構互叫的方式可能會造成無法預期的意外,比如簡單的file descriptor的指定,在parent和child process中分開執行,且一對多的時候運行的邏輯很怪,這部份可能要更深入底層才會理解這部分的坑。

R空間視覺化的學習資源(GIS or Spatial Visualization)

R在空間視覺化上在最近有很多的提升,尤其是整合如Google map api, leaflet等含有多種美麗地圖風格的package,這邊整理一下相關的學習資源。

 

R語言官方針對空間資料的文檔

CRAN Spatial Data Task View
裡面也有各種處理圖資的package可以使用


Roger Bivand. Applied Spatial Data Analysis with R. (2008). Springer
Chris Brunsdon, Lex Comber. An Introduction to R for Spatial Analysis and Mapping. (2015). SAGE

文章
R Spatial Cheatsheet
Zev Ross. Manipulating and mapping US Census data in R using the acs, tigris and leaflet packages.

Gene Transcription Regulation Database(GTRD):轉錄因子調控資料庫

在今年2017的NAR資料庫特輯中,GTRD資料庫整理人類和老鼠ChiP-seq資料並且彙整相關轉錄因子調控,使用四種負責處理ChiP-seq的peak caller程式,分別是MACS, SISSRs, GEM, PICS,並且將這些程式分析出的peak做clustering後,去掉重複的資訊。screenshot.png

目前使用起來不是很順手,每個要查看的轉錄因子都要一個個調整來看,有點不方便,似乎是使用biouml作為架構的,所以整體操作上頗有侷限性,當要查詢多個轉錄因子時會很麻煩,不過算是一個可以當作備用的轉錄因子調控資料庫,且是相對較新的,裡面用來估計peak左右可能影響的基因數量主要跟據來自於HOCOMOCO資料庫的資訊,這資料庫是同一個研究團隊在2013年的發表。screenshot.png

 

他本身可以直接下載它們整理好的interval檔案,資料欄位整理的很乾淨,並且提供那四種peak caller程式分析的資料可以挑選,或是最終其整合的版本。

screenshot.png

 

在進行變異辨認(Variant calling)前的序列品質校正(Base Quality Recalibration)流程

這邊主要是在記錄使用GATK中的RNAseq Variant Calling流程中很關鍵的一步,便是校正位點的品質資訊(Recalibrate base quality scores),這點為何重要呢?因為本質上我們就是在區分資料中單一位點是否有差異,其中最大的“障礙”,就是定序過程中的錯誤所造成的偽變異,如何去估算並且調整每一個位點下的品質數值是否是真的,就會決定我們變異辨認結果的準確程度。簡單來說,可以把Recalibrate Base Quality Score想成是Variant Calling的關鍵資料前處理步驟

關於GATK在Base Quality Recalibration的部分有詳盡的影片可以參考

實務面上重要的關鍵點:

理解在Read Group中會被用來矯正的是以lane和library為單位來思考,所以在進行Read Group的標誌的時候,會決定這一步的Recalibrate到底正不正確
(可以參考這邊GATK的文章)

在做Quality Score Recalibration主要有幾個步驟:
1. Add/Replace Read Group
把Read 資料貼上正確的Read Group(有的人會在alignment的時候就貼上,也可以在alignment後的bam檔來進行處理)

2. Markduplicates/CreatIndex(Picard)/Split’n’Trim(GATK) (調整Read alignmnet和建索引)

  1. 使用BaseRecalibrator(GATK),計算出矯正Quality Score的表

  2. 使用PrintReads(GATK),將上一步計算出的表來調整原本的Bam file,並且輸出矯正過後的Bam檔

閱讀參考:
1. GATK官方文檔和說明
https://software.broadinstitute.org/gatk/documentation/article.php?id=3891
http://gatkforums.broadinstitute.org/gatk/discussion/2801/howto-recalibrate-base-quality-scores-run-bqsr

VCF(Variant Call Format) 基因突變資料儲存格式

次世代定序技術出來後,緊接者是許多高通量資料的處理問題,從一開始單純看定序品質到後面用定序資料組成的基因體後,來看其變異資訊,此時就需要一個很好的儲存方式,這時候就有了VCF(Variant Call Format),這格式主要是由1000 Genomes Project所推動產生的,這計畫的目的希望用定序不同種族的正常人類來探討族群基因體間的異同。

一份VCF檔案大致長得像下面這樣
screenshot.png

基本上分成兩大部分:VCF header (描述所使用的工具、篩選標準、紀錄的細節名稱)和 body (變異資料本身)。screenshot.png

變異資訊基本上會有9行,分別紀錄CHR染色體號碼、POS位置、ID變異編號(是否在dbSNP中)、REF此位點在參考基因上的序列、ALT變異的序列、QUAL此序列的定序品質、FILTER是否通過所設定的篩選條件、INFO進一步的資訊、FORMAT格式資訊、SAMPLES各樣本在此位點的狀況。

會很驚訝於這樣的表示方法,其實就可以包含各種複雜的基因變異形式,如下圖所表現的樣子
screenshot.png

Reactome資料庫,可愛的視覺化小圖標(Reactome Icon Library)

之前就有稍微介紹一下Reactome,最近幾乎跟他形影不離,發現他不斷在提供一些超棒的小資源,其中就是可愛的生物小圖標,最近在練習畫複雜網絡的視覺化,他便提供了很好的素材。這邊下面的圖就是使用Cytoscape來視覺化我的項目的部分成果,整體空間的架構頗平衡的。
screenshot.png
這圖片中的背後的細胞圖標就是來自於Reactome Icon Library,頗佛心來著的。
screenshot.pngscreenshot.pngscreenshot.png

R資料處理小技巧:dplyr包中的summarise_each(funs()),自選函數來處理數值

使用dplyr的summarise_each提供自己寫函數用在每一行的方式,在做一些陣列資料處理還蠻方便的,簡單來說,可以用來把不同probe但屬於同一個基因的probe表現量合併起來,相關的函數,也可以使用summarise_all/summarise_at等等,靈活使用下頗為強大的。

tmp <- data %>%
        dplyr::select(-Date) %>%
        group_by(Id) %>%
        summarise_each(funs(sum)) %>%
        factor.summarise(.,2:67)

Neo4j簡介:"圖形"資料庫(Graph Database)

screenshot.png

會需要用到Neo4j是因為在剖析Reactome資料庫中的生物路徑資料時,發現他有提供Neo4J的資料庫備份檔案,可以直接下載後,匯入到Neo4j的資料庫,相對於Reactome提供的BioPAX檔案,置放在Neo4j中的資料schema比較偏向原生他們內部人員所設計的(在寄信跟Reactome內部人員討論時,他們也推薦我試試看,直接使用他們的Graph資料庫的備份),所以直接瀏覽和理解它內部的資料架構時,可以少轉一層,只能說有好有壞,相對於在處理BioPAX檔案格式時,直接使用Apache tdb/arq相關工具土炮來取得裡頭的檔案,Neo4J則提供蠻多操作選擇的,主要分成可以由命令行操作和網頁模式操作,而在讀取資料庫的時候,可以使用內建的Rest API直接去做資料的爬梳處理。

跟RDF/XML或是其他Triple store的資料庫,Neo4j比較偏向’單純’的圖形資料庫,也就是並沒有那麼重視資料間關係定義的“語意”,更注重在“資工面向”所在意的點,效能、操作介面、API接口、進階的Query方式,Neo4j有自己開發一套query的語言叫Cypher,相信在理解SPARQL後,其實在理解相關圖形資料庫的獨取值會顯得容易許多。

瀏覽器的使用介面

可以使用圖形介面來匯入資料庫,然後把Neo4j服務啟動
screenshot.png

可以直接登入Neo4j所開啟的服務,其有產生一個很好理解和管理的介面
screenshot.png

不得不說它內部的文檔非常豐富screenshot.png

也可以直接使用命令行來讀取資料庫,因為其有良好的API接口,當然可以用curl直接來對街,也能用Java, .NET, Javascript, Python, Ruby, PHP等寫好的函庫(本來有一個R的,但似乎原本負責人離職,後來就沒有支持了!)

curl -i -H "Content-Type:application/json" -X POST -u neo4j:neo4j http://localhost:7474/db/data/transaction/commit \
    -d '{"statements":[{"statement":"MATCH (p:Pathway)<-[:containedInPathway]-(re:Regulation) WHERE p.speciesName = \"Homo sapiens\" RETURN p.displayName,re.stId,re.displayName ;"}]}'

BioPAX(Biological Pathway Exchange)簡介:Pathway資料的交換標準格式

在做pathway analysis時,常會在各大資料庫中遊走,往往會希望能把自己關注的路徑資料下載下來仔細把玩,這時候會看到超多種的格式出現,像是在KEGG中儲存的KGML格式、SBML格式、SBGN格式等等,往往會讓大腦混亂,好在開始有計算生物社群的團體開始組織如何產生好的標準來交換生物路徑的資訊,那就是BioPax,一個Pathway資料交換的標準格式。

BioPAX(Biological Pathway Exchange)是一種用來描述生物路徑資料的語法,可以用來描述ㄒ訊息傳遞路徑、代謝路徑、分子和基因交互作用以及調控網絡(signaling pathway, metabolic pathway, molecular and genetic interactions and gene regulation networks).

本質上來說,BioPAX是奠基在RDF/OWL格式來發展的,只是使用符合生物路徑需求的描述語句(這部分可以參考前面寫RDF和OWL這兩種用來描述Linked Data的格式),BioPAX的建立也是在Open Biomedical Ontologies推動下所發起的,這是個在美國國家生物醫學數據模型中心補助下的計畫。

目前BioPAX在2010年公佈了level 3的標準,緊接者正在往level 4邁進。(其實level就是version的意思 XD)

那到底BioPAX長什麼樣子呢?可以看看下面這張來自The BioPAX community standard for pathway data sharing. (2010). Nature Biotechnology的圖!
螢幕快照 2017-05-17 上午11.15.30.png
右邊的語法就是BioPAX用來描述一個路徑的方式,假如之前有理解Linked DATA/RDF/OWL等文件格式的話,應該一看就會理解。

BioPAX基本語彙組成
screenshot.png

BioPAX基本的語彙可以看得出是以entity為主體,使用三大類的class來註解(為Linked Data/RDF的模式衍生的)

閱讀參考:
The BioPAX community standard for pathway data sharing. (2010). Nature Biotechnology