月份: 五月 2017

這邊主要是在記錄使用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)，這格式主要是由1000 Genomes Project所推動產生的，這計畫的目的希望用定序不同種族的正常人類來探討族群基因體間的異同。

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

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

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

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

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

這圖片中的背後的細胞圖標就是來自於Reactome Icon Library，頗佛心來著的。

使用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是因為在剖析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服務啟動

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

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

也可以直接使用命令行來讀取資料庫，因為其有良好的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 ;"}]}'

在做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的圖！

右邊的語法就是BioPAX用來描述一個路徑的方式，假如之前有理解Linked DATA/RDF/OWL等文件格式的話，應該一看就會理解。

BioPAX基本語彙組成

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

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

NCBI所掌管的Pubmed已經是全世界最大的醫療資料庫，近年來提供很多直接用程式串接的工具，但很重要的地方便是他們對於request的規範，仔細遵守下才可以避免被鎖IP的處置，下面是這個最常用的Pubmed E-utilities的相關規定：

主要有五大重點：
第一個：使用https://eutils.ncbi.nlm.nih.gov/entrez/eutils/作為base url
因為Pubmed內部管理負載平衡和使用效能的考慮，把不同服務的url都乾淨的分開，雖然你也可以直接用www.ncbi.nlm.gov開頭的url來爬，但這樣是會關注的。

第二個：在使用E-utility的URL request頻率規範
基本速限是每秒3個requests以下，在週末和平日的晚上九點到凌晨五點希望避開大量使用，假如是你開發的工具會很長期的使用，那麼可以跟NCBI註冊你的工具名稱和email，關於IP被禁的問題，可以禮貌地寄信給eutilities@ncbi.nlm.nih.gov，跟他們說明你的狀況，和IP。

第三個：增加每次request的內容量
盡量減少request的頻率，增加單次request的內容量，這樣可以比較有效的提高謝率，尤其是在使用EFetch這種下載全文資料或是蛋白質序列等的工具。

第四個：注意資料版權的問題
在使用E-utilities來開發產品的時候，要注意NCBI對於版權等資訊的規範，可以點選這邊來看細節，另外，關於在Pubmed查到的部分abstract內容也是受到美國或是其他國家著作權法的保護的。

第五個：注意特殊字符
在使用E-utilities時，會用完整的url來索取資源，這時候要注意特殊的字符和全半形的問題。

這次Eric Topol的One-on-One節目採訪到了Brian J. Druker，Brian Druker醫師是發展出第一個非常成功的標靶藥物Gleevec的人，也是Oregon Health & Science University裡面算是第一把手的科學家，最近在成功完成Niket創辦人的捐款挑戰Knight Cancer Challenge後，取得Phil Knight夫婦5億美元的捐款，加上對半跟超過1000名捐款者募集的5億美元，共十億美元的資金，依此成立了一個新的Knight Cancer Institute，預計希望推動一個新的計劃Precision Earily Detection of Cancer Program。在這筆資金下會招募約250-300個科學家和相關的醫師，以及至少25位國際頂尖的研究者。

http://c.brightcove.com/services/viewer/federated_f9?isVid=1&isUI=1

這邊有一份蠻有趣的關於 Precision Prevention and Early Detection Working Group Recommendation的資料，google Earily Detection 找到的，不知道有無關聯性，但提供一個不錯對於未來領域的想法。可能是NCI對於未來癌症規劃所做的計畫書。

開始整理一些使用owl或是rdf儲存的資料時，直接使用像是jena arq來做sparql，會發現效率非常的差，所以假如真的要將一套基於rdf/owl資料使用sparql查詢的應用華，必定要思考找尋一個好的database來專門解決效率的問題。

畢竟rdf標準發展也有10幾年了，其實與其相關連的資料庫其實比想像中的多，但要如何挑選適合的便是件頭痛的事情，從下面這張表可以一窺rdf相關資料庫的發展和相關的database system的演進。

 Apache’s Jena TDB system is stated as being faster, more scalable, and better supported than the Jena SDB, which is a non-native system relying on an RDBMS. TDB is, for instance, the system supporting persistence in the Fuseki SPARQL server. The architecture is built around three concepts, namely a node table, triples/quads indexes, and a prefixes table. The node table serves to store the dictionary and follows the two mappings approach presented in Chapter 4. Practically, the string-to-id and id-to-string operations are respectively implemented using B+trees and a sequential file. A large cache is dedicated to ensure fast data retrieval during query processing. Triples and quads indexes are stored in specialized structures and respectively store three and four identifiers from the node table. B+trees are used to persist these indexes. The system supports SPARQL update operations, which are handled using ACID transactions (with the serializable isolation level) through a WAL (write ahead logging) approach. This implies that write transactions are first written into a journal and then stored in the database when resources permit it. This approach presents the benefit of not requiring a locking solution for read transactions. Finally, Jena TDB supports a bulk-load solution that does not support transactions. The different features contained in Jena TDB, such as some security aspects as well as some APIs, make it a solution to consider in a production setting.

 

閱讀參考
1. Olivier Curé; Guillaume Blin. RDF Database Systems. (2014). Morgan Kaufmann

在設計Linked Data相關的schema時，最好能使用那些大家常用的，好處是能減少時間，提高開發效率，最重要的則是之後你所設計出來的database是可以跟其他用同樣schema儲存的資料做合併甚至提高一些inference的效果。

而Linked Open Vacabularies是一個搜集大家常用語料庫的地方，所以也可以在這邊找尋適合建立自己資料ontology的vacabularies。

這邊可以看到裡面大圈圈就是使用率很高的Vacabularies，像是foaf, vann, skos, dcterms, dce, cc等。

可以點選自己想要進一步了解的Vacabularies，他會有詳盡的URI、 Namespace和相關創建者的資訊，簡易乾淨的頁面讓你知道這個語料庫裡面的class, properties有幾種，比如下面這個很著名的foaf語料，是專門用來處理社交網絡的資訊關係。