2017 24th Nucleic Acids Research 資料庫特輯

screenshot.png
每年Nucleic Acids Research期刊的資料庫專刊(database issue)都會整理整年新增加的生醫資料庫和有哪些資料庫有更新,能刊登在Nuclei Acids Research的資料庫都有一點的品質,也是找尋可以利用在自己相關領域的素材。

今年的Database issue總共包含152篇論文,有54個新增的資料庫,另外,98個資料庫有更新。這些資料庫的內容有跟基因體結構、基因表現量、基因調控資訊、蛋白質交互作用等。其中有三個資料庫被特別提出來說明其“breakthrough”的貢獻,分別是denovo-db, Monarch Initiative, Open Targets。

“Breakthrough"的新資料庫

denovo-db

screenshot.png

資料庫主要是收集人類相關的基因變異,總共收集40個研究、23098 trios,共32991個de novo variants,每個variant的資料除了基本的染色體位置、變異種類外,還有在轉錄和轉譯層面的資訊,還有嚴重度分數和頻率,目前驗證階段,還有來自於個人的phenotype資訊。整個資料庫的資料都可以被下載供研究使用。

Monarch Initiative

計畫(說資料庫的話,把他們的企圖講得太小了)主要是想大規模地將其他物種的資訊希望導入到人類疾病和相關機制的推測上,讓genotype-phenotype間的關係能有更多資訊供參考。
Open Targets

screenshot.png

這平台提供藥物標的開發的搜尋和視覺化工具,整合疾病跟潛在藥物標的之關係。

這三個平台都提供API串接資料,可以看出目前生醫界在處理資料的規格慢慢貼近軟體業界的水準

值得注意的資料庫

名稱 類別     介紹    論文
GTRD 轉錄因子調控 紀錄轉錄因子間的交互關係 link
TcoF-DB 轉錄因子調控 記錄轉錄因子和轉錄因子co-factor的交互關係 link
The Gene Ontology 基因功能註釋 基因功能註釋 link
PANTHER 蛋白質功能註釋 跨種蛋白質功能註釋 link
InterPro 蛋白質註釋 蛋白質功能分析,且由不同的蛋白domain link
Protein Ontology 蛋白質註釋 蛋白質功能詮釋 link
OntobeeOntobee link資料庫 有各種ontology結構的資料庫資訊 link
KEGG 代謝路徑 仲要代謝路徑資料庫 link
STRING 蛋白質交互 重要蛋白質交互作用資料庫 link
BioGRIDBioGRID 代謝路徑 重要代謝路徑資料庫 link
pathDIP 訊息傳遞路徑 將現有的訊息傳遞路徑加入但蛋白質間交互的關係來延昇|link
TissueNet 蛋白質互動 蛋白質間的交互關係 link
XTalkDBXTalkDB 訊息傳遞(crosstalk) 註釋不同訊息路徑間的關係 link
Exposome-Explorer 環境飲食基因體學 生物標記跟環境風險因子相關 link
Pharos 藥物蛋白體學 特定疾病相關的蛋白質標的 link
Bio-TDS 生醫資料文字探勘 用自然語言處理技術支持文本式搜尋 link

Gene Set Enrichment Analysis(GSEA)資料格式介紹:gct, cls, gmt

資料種類 內容 格式
基因表現量資料Expression dataset 主要為基因表現量資料 res, gct, pcl, txt
實驗設計資料Phenotype labels 主要為實驗設計,標記各樣本 cls
基因集Gene sets 要用來檢定的各個基因集(可以直接使用GSEA軟件的) gmt, gmx

注意事項:
1. 基因表現量資料中的基因名稱跟基因集都要使用HUGO中基因的名稱,假如使用的是microarray的資料,要注意使用的是哪個平台,要提供那平台的probe註釋(大部分affymetrix的都有內建資料,所以不用太擔心)

基因表現量資料
GCT:Gene Cluster Text
為tab分隔的資料,可以使用純文本編輯器或是excel,不太推薦使用excel,有時候會出現一些數字被轉換成日期的問題

基本原則:

  1. 第一行一定要有#1,2這字串。
  2. 第二行總共要有兩個數字,第一個數字代表這筆基因表現資料總共有幾行,第二行數字代表有幾個樣本。
  3. 接者便是基因表現量資料,第一列為基因名稱或是probe編號,第二行為這些基因或是probe的描述,再來就是所有樣本個別的表現量多寡。
    screenshot.png

實驗設計資料
CLS:Categorical/Continuous files format
screenshot.png
主要用三行表示
1. 第一行有三個數字,第一個數字代表樣本總數,第二個數字為總共有幾個類別,第三個數字都一定為1
2. 第二行要以#開頭,然後有兩個數字,緊接者為各個類別的名稱
3. 第三行則是各個樣本的類別種類

注意事項:
1. 樣本類別也可以是連續變數,比如時間序列的,可以用下面的方式來表示
#numeric
#IncreasingProfle
30 60 90 120 150

基因集資料
這邊介紹兩個主要基因集資料的格式,分別是GMT和CMX,兩個資料格式有兩個差別:
1. GMX基因集資料是以列方式儲存
2. GMT基因集剛好反過來,用行來儲存,適合儲存大於256個以上的基因集
GMX: Gene Matrix file format
screenshot.png
GMT: Gene Matrix file format
screenshot.png

論文閱讀:Tumor neoantigenes:building a framework for personalized cancer immunotherapy 精華簡報

 免疫療法在癌症治療上的角色越來越重要,在化學治療、標靶治療、手術治療、放射治療等等之後,癌症本質上基因變異的多樣性,在使用上面幾種療法下效果較差,而在免疫療法的策略下,讓自己的免疫細胞去殺死癌症,似乎可以解決上面幾個方法的死角,當然免疫療法也有他的局限性,他很大程度的需要依賴自己身體的免疫系統,另外,過度激活的免疫系統也會產生自體免疫的問題,當下最大的難關是很難預測病人的反應性,有的對治療反應性好的,癌細胞可能會戲劇性地消失,如何預測這件事情變得很重要,這篇論文的重點著重於免疫療法中其中一塊,針對個別癌症特異性抗原來當作下手目標。

MedRec:區塊鏈技術應用到病歷系統的實作

MedRecMedRec是由MIT Media Lab中研究生 Ariel Ekblaw, Asaf Azaria , Thiago Vieira在資深研究人員Andrew Lippman帶領下所提出基於區塊鏈在病歷資料上的應用,比較有趣的一點是他們加上一層讓研究者可以參與的方式近來,頗符合目前醫療資料複雜需求中的重要要角。而將區塊鏈應用在醫療資料處理上最重要的目的是提升所謂的interoperability,也就是如何讓單一個病人的資料原本散佈在各個醫療機構中,能有效互通,減少重複檢查的問題,另一方面,也開始逐間導向以病人為中心的儲存方式,因為許多健康資料的來源已不限於醫學中心,大量消費產品已經能讓病人有足夠多自主權或許各式各樣的資料,如23andMe等

This lack of interoperability costs 150,000 lives and $18.6 billion per year, according to the Premier Healthcare Alliance.

詳細的實現內容可以看這篇論文的細節,整體上這他們系統設計的原型有嘗試在跟Beth Israel Deaconess Medical Center內來試驗,希望能將醫院裡現有的資料匯入到MedRec中。

下面的圖表可以看出她們整個設計的架構,主要牽涉到的模組有:


1. smart contract infrastructure
– registrar contract
– patient-provider relationship contract
– summary contract
2. backend API library
– PyEthereum
– PyEthApp

3. Database getekeeper
4. EHR Manager
– ONC’s Blue Button design

screenshot.pngscreenshot.png

參考閱讀:
MedRec: Electronic Medical Records on the Blockchain

MedRec: Medical Data Management on the Blockchain

Asaph Azaria, Ariel Ekblaw, Thiago Vieira, (2016). Andrew Lippman. MedRec: Using Blockchain for Medical Data Access and Permission Management. 2016 2nd International Conference on Open and Big Data

R資料處理小技巧:關於missing data標記,dplyr包中的coalesce/na_if

dplyr中的函數內借鑒許多資料庫操作的函數而來,這邊要分享的一對函數便是這種關係,我們常見一個情況需要將資料中的NA取代成特定得值,或是將某些值當作NA處理,這邊在R dplyr中的coalesce和na_if便是很好的幫手,coalesce函數就是直接借來在SQL中同名的函數名稱在這邊使用,他的功能也是用來做NA取代,放入估計的函數,這邊便來釋放一些例子:

coalesce
coalesce接受vector,然後將第一個遇到的NA值取代成為特定的數值,要發揮它最大的效果可以跟purrr函數的map一起使用,便能大幅增加他的適用性。

# Use a single value to replace all missing values
x <- sample(c(1:5, NA, NA, NA))
coalesce(x, 0L)
# Or match together a complete vector from missing pieces
y <- c(1, 2, NA, NA, 5)
z <- c(NA, NA, 3, 4, 5)
coalesce(y, z)

na_if
na_if(x, y)剛好跟coalesce相反,輸入的vector x中,只要跟vector y中相等的就會被替換變成NA值。

na_if(1:5, 5:1)
x <- c(1, -1, 0, 10)
100 / x
100 / na_if(x, 0)
y <- c("abc", "def", "", "ghi")
na_if(y, "")

useR!2017 比利時布魯塞爾

screenshot.pngscreenshot.png

2017年的useR!會議將在美麗的比利時舉辦,聽這次的主辦人講,吸引大家來比利時唯一的東西便是生產全世界最多且最美味啤酒的地方了,今年會議的舉辦時間為7/4-7/7號,此次由比利時分析公司OpenAnalytics的Tobias Verbeke、Hasselt University統計系的Ziv Shkedy教授、英國Warwick University統計系的Heather Turner、OpenAnalytics的Matthias Verbeke作為籌辦委員。

今年會議的教學主題有下面這些:

Joaquin Vanschoren, Heidi Seibold and Bernd Bischl: OpenML: Connecting R to the Machine Learning Platform OpenML
Stephanie Kovalchik: Sports Analytics with R
Karline Soetaert, Thomas Petzoldtenvironmental modeling using R.
Hana Ševčíková: Introduction to parallel computing with R
Martyn Plummer: Introduction to Bayesian inference with JAGS
Gabor Csárdi: R package development with R-hub
Charlotte Wickham: purrr
Edzer Pebesma: Spatial data in R: new directions
Dirk Eddelbuettel: Extending R with C++: Motivation, Introduction and Examples
Colin Gillespie: Efficient R Programming
Toby Dylan Hocking, Rebecca Killick: Introduction to optimal changepoint detection algorithms
Francois Michonneau, Tracy K. Teal: Data Carpentry: Open and Reproducible Research with R
Arun Srinivasan: data.table for beginners
Signe M. Jensen, Christian Ritz: Dose-response analysis using R
Bhaskar V. Karambelkar: Geospatial visualization using R
Taylor Arnold, Lauren Tilton: Introduction to Natural Language Processing with R

藥物治療分類代碼系統ATC code

在爬梳藥物基因體相關資料時,很重要的便是去理解藥品的編碼系統,因為同一類藥品通常名字有很多種,不同國家的藥物編碼也不同,常混者商品名等等,使用編碼系統可以免去名稱混亂所造成的問題,且要是有一套國際使用的編碼系統變可以使得國間的資訊可以共通。目前世界常用的藥物分類編碼系統有兩種,一種是AHFS(American Hospital Formulary Service)/DI(Drug Information)和WHO(World Health Organization)/ATC(Anatomical Therapeutic Chemical Classification)。傳統上使用這套系統便能做很好的藥品項目管理、評估用藥品質、效益分析等等(因為有統一的代碼,所以各國資料可以彙整在一起,所以當要把基因資訊跟藥物彙整,當然會需要使用到藥物編碼系統)。

這邊來記錄一下ATC code的系統:

一個藥物的完整ATG code總共有七位由英文和數字組成,由左到右,不同位置代表的東西不同,共分成五組,分別代表一個藥物所主要作用的器官系統、治療哪類人體功能問題、藥物化學組成的成分、藥物動力學上的組別來表示:

第一層:第一位,為一個大寫字母,總共有14個字母,代表藥物主要作用的目標人體器官系統

A: Alimentary tract and metabolism
B: Blood and blood forming organs
C: Cardiovascular system
D: Dermatologicals
G: Genito-urinary system and sex hormone
H: System hormonal preparations, excluding sex hormones and insulins
J: Anti-infectives for systemic use
L: ANtineoplastic and immumodulating agents
M: Musculo-skeletal system
N: Nervous system
P: Antiparasitic products, insecticides and repellents
R: Respiratory system
S: Sensory organs
V: Various

第二層:第二三位,由兩個數字組成,代表藥物所屬於的治療功能類別

ex: C03 Diuretics

第三層:第四位,為一個大寫字母,代表藥物在藥物動力上的組別

ex: C03C High-ceiling diuretics

第四層:第五位,為一個大寫字母,藥物在化學特性上的組別

ex: C03CA Sulfonamides

第五層:第六七位,為兩個數字組成,藥物實際所屬於的分子特性組別

ex: C03CA01 Furosemide

GSEA(Gene Set Enrichment Analysis)分析方法和觀念介紹

screenshot.png

Gene Set Enrichment Analysis的方法是目前在pathway analysis方法中,歸類為functional score analysis裡的state of the art(可參考此篇),且其搭配的工具和gene sets都有蠻完整的支持,值得深入瞭解,並納為自己的口袋武器。實際上,Gene Set Enrichment Analysis有許多係分支的方法,但這邊主要談的是由Broad Institute裡的GSEA這款軟體。

基本觀念

所有的pathway analysis的觀念都是用一組特定的gene set簡化代表特定生物內的pathway或是功能,然後看這組gene set是否有在我們的資料中有特別顯著的表現。

不管是在做pathway analysis或是gene set enrichment analysis(其實大同小異,做時候的觀念差異不同而已),到後來會發現如何挑選gene set也變成是一種功夫,必須依靠對背景的生物問題有透徹的文獻研究才會能有較詳盡的分析結果。可以研讀Hammer lab中對於自制gene list的一些看法,會對觀念的釐清有頗大的幫助,其實在做pathway analysis的本身就是一種生物高通量資料探勘的過程。

要如何做GSEA呢?

(這邊以由Broad institue維護的GSEA使用來分享)

第一步理解GSEA所需要輸入的資料格式:*.gct, *.cls

GCT:Gene Cluster Text

screenshot.png

CLS: Categorical class

screenshot.png

第二步挑選或是設定一組基因列表S(關鍵在於priori defined)

這群基因可以代表某個metabolic pathway或是共有特定的GO term,想要測定這個特定的pathway或是go function group有沒有在整個L基因列表中有表現顯著,在GSEA論文中所使用的gene set內大小至少35個gene,不過很多人自製的gene set並不一定會超過這數量。(在GSEA的軟體裡能直接抓取MSigDB資料庫screenshot.png)screenshot.png

第三步製作一個根據基因在兩組phenotype中表現值差異來排序的基因列表

量測兩組樣本間某個基因表現量的差異,有許多作法,在GSEA官方軟體內所使用的主要有五個指標,用來表現基因在兩個表型間的差異,因為之後排序就是依照這些值,所以最好理解一下他們各自的算法,分別是:Signal2Noise(預設方法), tTest, Ratio_of_Classes, Diff_of_Classes, log2_Ratio_of_Classes。
下面則是他們的公式

Signal2Noise

\frac{\mu_{A} - \mu_{B}}{ \sigma_{A} + \sigma_{B} }

tTest

這方法至少每一組的樣本都要有三個會比較適合
\frac{\mu_{A} - \mu_{B}}{ \sqrt{\frac{\sigma^{2}_{A}}{n_{A}}+ \frac{\sigma^{2}_{B}}{n_{B}}}}

Ratio_of_Class

其實就是所謂的Fold Change,兩者差異越大,數值就會偏離1
\frac{\mu_{A}}{\mu_{B}}

Diff_of_Classes

直接用兩者的差值
\mu_{A} - \mu_{B}

log2_Ratio_of_Classes

將計算好的fold change取對
\log_{2}\frac{\mu_{A}}{\mu_{B}}

第四步計算Enrichment Score(ES)

整個GSEA分析中,最關鍵的觀念就是Enrichment Score的算法理解,在排序好的基因列表L中,我們所定義的S genes列表中的基因理論上會散佈在其中,此時由左到右,計算累計的統計量(running-sum statistic),假如遇到的基因是在S genes中的,那麼就加{P}_{hit}(S,i),假如不是,那麼就減{P}_{miss}(S,i) ,最後在過程中得到的最大值就是所謂的Enrichment Score
{P}_{hit}(S,i) = \sum_{{g}_{i}\in{S}}^{} \frac{|r_{j}^{p}|}{N_{R}} , where \quad N_{R} = \sum_{{g}_{i}\in{S}}|r_{j}|^{p}

{P}_{miss}(S,i) =\sum_{{g}_{i}\notin{S}}^{} \frac{1}{N-N_{H}}

第五步計算ES有無顯著

會發現這其實就是一種修改過的無母數統計做法,比較重要的是他比較的是當我們隨機把phenotype的標籤做置換,然後經過permute後,比較兩者有無差異的方式。

第六步多重比較的事後調整(adjustment)

當我們有很多組S Genes要看有無顯著時,這時候就有多重比較的問題,這邊主要使用False Discovery Rate來控制。

 其實假如直接下載使用Broad Institute維護的軟體GSEA話,只要做第一個步驟就行了,輸入資料後接下來就不用擔心了,這是理解這六步驟才能真的懂GSEA分析的邏輯。

下面展示一些使用GSEA分析後的輸出檔案,其實幫助很大,裡面內建的gene set非常多,也可以由這些gene set來建立一些生物感覺。

heat_map_245.png

screenshot.png

screenshot.png閱讀參考:
Aravind Subramanian, Pablo Tamayo, Vamsi K. Mootha et al.(2005). Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. PNAS, vol 12, 15545-15550

R語言套件:data.table,可處理超過100gb大小資料

之前做需要幾萬次重複的資料處理時,使用dplyr基本的data.frame操作語法,單一個運算單位包含了多組的pipe操作,大約五萬次運算,加上一些平行語法,後來發現盡然會讓R吃的資源爆表,仔細研究才發現dplyr在執行多組data.frame運算,搭配上pipe時,其實會進行資料的copy,幾層下來,很容易讓吃的ram大增,所以需要找一種在運算上較乾淨的套件,同時也能做複雜資料操作,而data.table就是為這件事情而生的,他可以操作100gb大小的data.frame處理,對於此量級資料的join, order, update, delete都可以有很高速的表現,相對於dplyr將資料操作抽象成比較貼近人思考的select,filter,mutate,summarize等,data.table的語法更貼近SQL一點。

整個Data.Table的語法概念可以由下面一張圖直接解釋
screenshot.png

Take DT, subset rows using i, then calculate j grouped by by

Data.Table懶人表

閱讀參考:
[心得] 資料整理套件介紹-第一章 data.

有趣的各種data.table語法跟其他套件的比較
plyr::ddply vs data.table::rbindlist
R data.table – collating data
R data.table – subsets
R data.table – sum by groups
unnest in R data.table

datacamp的data.table語法整理
data.table cheat sheet

Frequently asked questions
Introduction to data.table
Keys and fast binary search based subset
Reference semantics
Efficient reshaping using data.tables
Secondary indices and auto indexing

Health Level 7(HL7)研究筆記:五大基本架構

HL7基本上是一個試圖把整個醫療體系的行為資訊化的架構,所以非常龐大,要了解他的話可以分為五個部分來分別研究,這五個部分:Foundation, Specification Infrastructure, Implementation Technology Specifications, Services, Domains。分別把HL7整個很系統性的方式分類來理解,而不會被他複雜的內容嚇到,可幫助架構好對HL7的理解。

基礎內容Foundation :資料結構和名詞的定義

既然HL7是一套把醫院體系裡所有的行為資訊化,那麼就必須定義清楚每個醫療體系中參與的人事物,和溝通的資訊,並且有仔細的編碼,所以HL7有一套定義Data Type(Abstract, Vocabulary, Refinement, Constraint, Localization),在version 3還最重要的提供了Reference Information Model(RIM),有了這模型(完整使用軟體工程的概念),可以視覺化所有流程如下面這張圖。screenshot.png

系統架構Specification Infrastructure

整個HL7最精華的部分就是其底層架構的設計,因為HL7建立的目的就是為了讓醫療系統間和內部的資訊流通能有一個共通的基礎,所以裡頭所謂的Messaging(訊號傳遞)機制的設計是整個HL7的重點,包括訊息的標題和內容該如何規範等,訊息該怎樣設計傳輸到正確的角色上,訊息能寄給誰,誰有權限收到什麼樣的訊息,兩個物件間的溝通該怎麼開始、保持通訊和結束,過程中要用什麼來代表物件或是機構(OID)。

screenshot.png

所使用的現有技術Implementation Technology Specifications

實現整個資料傳輸的架構主要就是建立在XMLUML(Unified Modeling Language)上,XML主要用來紀錄資料本身(在Data Type那或是messaging裡便是用XML的格式),而UML則是紀錄或是表示整個系統中的物件、事件、角色等(Reference Information Model便是以UML的方式建構)。

screenshot.png

screenshot.png

HL7所提供的服務Services

Common Terminology Services是HL7提供用來轉換各種現有資訊平台的架構,讓他們能銜接到HL7 標準的服務,只要使用他們所提供的API,即使不使用HL7所用的資料結構,也能轉換成HL7(絕對有很大的限制性)

所採用的共通語彙Domains

HL7標準是為了轉換醫療體系中的人事物進入可紀錄和交換的資訊,定義每個細節便很重要,但是很多醫療體系中所使用的概念都已經很完整了,這時候HL7便會直接拿來用,而不用重新定義。