月份: 三月 2017

藥物治療分類代碼系統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)研究筆記：五大基本架構

1 則迴響

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

實現整個資料傳輸的架構主要就是建立在XML和UML（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便會直接拿來用，而不用重新定義。

使用R dplyr串接sqlite資料庫

1 則迴響

最近十年R在很多方面都有進步，其中一個就是加強跟資料庫的串接方式，DBI就是其中一個奠基於Relational Database的對接方式，觀看R Special Interest Group on Databases名單：Hadley Wickham, Kirill Müller，都是頗為厲害的高手。

當然，當資料量還能塞在記憶體裡時，用資料庫的方式反而會浪費儲存空間，操作速度也不會過快，但當資料量大過於記憶體可以支撐，或是有一些同步資料等的需求時，利用R和資料庫的連接變很棒！

其中dplyr在這方面就默默做得超棒的，大多數人都知道dplyr對於資料的處理很方便，可以快速做select, filter, mutate, summarise等，但另一方面更棒的地方是他也能很好的跟資料庫對街，sqlite算是最簡單可以入手的，可以很快速的跟現有的工作流程整合在一起。

目前dplyr支持sqlite, mysql, postgresql這三大資料庫類型，基本上很夠用。他的目的是減少在R和SQL這兩類語言間換來換去的麻煩，另一方面，設計一個機制讓你可以安全地在R裡面用dplyr操作資料的verse轉換成SQL，所以相對來說可以轉換較順利，另一方面，他會有一個預算處理的機制（通常處理前幾行），所以可以簡單的預看自己的操作，等確定要撈資料的方式在執行這些query。

最主要提供的功能可以分成三大類：
1. 基礎database的串接(下面的代碼)
2. 使用dplyr的語法來進行sql query
3. 效調與資料庫的query:collect, compute, collapse

# 創建一個新的資料庫
my_db <- src_sqlite("Test.sqlite", create=T)
# 連接現有的資料庫
my_db <- src_sqlite("OldeDb.sqlite", create=F)

# 查看資料庫裡有哪些tables
src_tbls(my_db)

# 形成一種query status，沒有實際查看全部的資料
query <- tbl(my_db, "table_1")
# 實際查找資料庫，回傳資料到R裡面(利用collect)
query.result <- tbl(my_db, "table_1") %>% collect()
# 將R裡面的data.frame塞到資料庫裡成為table
copy_to(my_db, rdataframe, name="rtable", temporary=F)

#建立一個query物件
flights_sqlite <- tbl(nycflights13_sqlite(), "flights") #也可以直接使用sql的語法 tbl(my_db, sql("SELECT * FROM flights")) #奠基在dplyr語法來進行query select(flights_sqlite, year:day, dep_delay, arr_delay) filter(flights_sqlite, dep_delay > 240)
arrange(flights_sqlite, year, month, day)
mutate(flights_sqlite, speed = air_time / distance)
summarise(flights_sqlite, delay = mean(dep_time))

閱讀資料：
1. How to Use DBI: Connecting to Databases with R
2. DBI: R Database Interface
3. A Common Database InterfaceA Common Database Interface
4. Databases and dplyrDatabases and dplyr

使用R處理時間資料（DateTimeClasses）的格式（lubricate, POSIXlt,POSIXlc）

發表留言

「lubridate r package」的圖片搜尋結果

時間資料的處理，假如不借助現有的“時間類別”資料型態的話，會頗麻煩，比如有兩個時間要比看看誰先誰晚，想知道兩個時間的間隔，或是要進一步要篩選一堆時間資料根據某個特定的時間點，要自己用numeric來實現會麻煩一點，所以使用內建處理時間的函數像是as.POSIXlt, as.POSIXlc, as.Date會方便很多。他們本身就可以實現時間的加減等，另外，也可以使用lubricate軟件包，這是一個更為方便用來處理時間資料轉換的工具。

仔細研究會發現在電腦中的時間顯示也不是一件簡單的事，可以看此篇unix time，最重要的觀念是理解POSIX和epoch這兩種時間表達的方式，尤其是Epoch，其是所謂的reference time，它的目的是幫助簡化“時間值”彼此的加減運算，等於把時間轉換成距離某個“參考時間點”幾小時、幾分或是幾秒，那這樣就可以輕鬆地對兩組“時間值”彼此直接加減，不需要經過兩次換算，但這邊也容易造成使用者混淆，因為每個平台或是程式語言可能有彼此不一樣的epoch參照點，像是R的參照時間就是1970/01/01，所有的時間假如直接轉成數字的話，便是距離此時間點的天數（或是小時），端看時間裡所使用的最小單位。

下面是當我們想要進一步瞭解後面細節遇到的名詞：
1. POSIX(Portable Operating System Interface)
這其實是IEEE的標準中為了讓各個操作系統能保持一致的，裡頭包含了許多系統底層的規範包括Process, Signal, Memory violation, Pipe, C library, 當然timer也是其中的規範之一
2. Epoch(reference time)
Epoch是傳統年代表中用來參照的時間點，比如西元前和西元後中的epoch就是耶穌誕生日，而在計算機領域中不同的操作系統會有不同的epoch，像MATLAB就是以Jan 1, BC、GO則是用Jan 1, AC，Microsoft則是1990-1-0（真的是Jan 0），大多數的系統和語言像是JAVA, PHP, Ruby, Javascript, Perl，unix系統等則是使用1970,1,1。

R系統中DataTimeClasses主要有兩個物件：POSIXct, POSIXlt
1. POSIXct代表的是總秒數從1970/1/1開始算起(R裡面的epoch)
2. POSIXlt則是一個列表，裡面有跟時間相關的數值（秒、分、時、日、月、年等）

基本R時間數字的處理和轉換
可以將R內時間資料處理分成兩個大觀念來理解：資料格式和常用操作

在資料格式上，R裡面有兩大類：Date和POSIXct/POSIXlt
POSIXct/POSIXlt為承襲unix系統中處理時間的方式，自由度較高，也較複雜，使用的函數為as.POSIXct/as.POSIXlt
Date為簡單版的時間資料格式，就是由1970/01/01作為epoch的，使用的函數為as.Date

在處理時間的資料上，觀念上其實只有兩類操作：
1. 時間“資料格式”轉換顯示方式（通常被忽視）
2. 非時間“資料格時”（為字串）轉換為時間”資料格式“
但通常沒有區分開這兩類時，很容易會感到挫折，因為…會發現函數的表現跟自己預期不太一樣，或是搞不清楚文檔在說什麼。當可以區分這兩類操作後，再來則是理解一點關於“電腦世界”制定時間的一些方法，才不會被一堆名詞搞矇，接者可以開始找一些不錯的函數包lubricate在處理一些時間資料的操作。

#第一類：本為時間資料格式，拿來轉換顯示格式
Sys.time()
#[1] "2017-03-19 23:14:28 CST"
class(Sys.time())
#[1] "POSIXct" "POSIXt"
time <- Sys.time()
format(Sys.time(), "%Y-%m-%d")

#第二類：本來是“字串”，要轉換成時間資料格式
as.POSIXct("2014-05-01", format="%Y-%m-%d")
[/code
或也可以使用lubricate函數來處理文字的轉換

ymd("20110604")
mdy("06-04-2011")
dmy("04/06/2011")

轉換成時間資料格式後，很多處理就變成很直覺性了，你可以直接加上秒數，或是兩個時間相減或是相加

Sys.time()
#[1] "2017-03-20 21:23:12 CST"
Sys.time() + 1000
[1] "2017-03-20 21:38:16 CST"
as.Date(Sys.time()) + 10
#[1] "2017-03-30"
as.Date(Sys.time()) + 15
[1] "2017-04-04"

閱讀參考：
1. Do more with dates and times in R with lubridate
2. Converting time zones in R: tips, tricks and pitfalls
3. Date data class
4. Date-Time class

bookdown: 用R寫書的好工具

發表留言

近幾年來在R語言界越來越講究所謂的reproducible of analysis，當生物醫學的資料分析越來越複雜，牽涉的團隊越來越多的時候，讓分析的過程和文件能容易分享開始引起大家注意，所以緊接者如pandoc, knitr, rmarkdown等等讓人直接將分析代碼和註解文字寫在一起，甚至後面更進階的shiny，都有想要達成讓分析流程可以更容易分享，且讓不懂代碼的人也能接近，反過來，寫書算是一個最大型的書寫成果，要是裡頭含有代碼和運算（這年頭生物資訊的書不可能純生物，幾乎夾雜者數學公式、運算和軟體包）那這本書就會有點難處理，所以knitr的作者Yihui Xie，開發了這個bookdown，以Rmarkdown為撰寫主體，並且整合htmlwidge, shiny app, citation reference等，可以滿足之前單用latex無法做到的書寫體驗。最後使用bookdown寫完後，可以輸出成html, pdf, epub等格式，還可以上傳到bookdown.org網站上面，最棒的是跟RStudio能整合在一起，基本上創造了新的書寫模式。

有在用Rmardown的人在使用bookdown的話，上手會非常快，而且作者也有寫一本書bookdown: Authoring Books and Technical Documents with R Markdown，把如何使用bookdown用的方式很清楚的解釋和教學。

想要快速地開始，可以直接下載bookdown的template來修改：
1. 中文版型
2. 英文版型

screenshot.png

進階R語言爬蟲：在rvest或是httr中使用proxy和timeout 連線問題的處理

發表留言

在撰寫爬蟲程式的時候，有一個很重要的問題要注意便是當要爬取特定網站的大量資料時，通常網站管理員會對你的爬蟲做“處理”，針對你的ip做禁止，而每個網站處理的方式不同，可以看這些資料：
How do websites block web crawlers?
Web Crawlers: Love the Good, but Kill the Bad and the Ugly

看完就會知道原來…自己的爬蟲會遭遇什麼可怕的事情，這時候處理注意一些連線數和頻率的問題外，可以使用proxy來保護自己的ip，至少被禁止後還可以換個proxy來爬，另一方面，使用proxy最大的問題就是連線不穩定，所以就會需要寫一些error handler的機制，這邊主要分享一下使用httr和rvest時，如何使用proxy，另外，還有分享使用tryCatch來處理connect problem。

在httr包中，有use_proxy這個函數，可以用來設定httr::GET裡面的網路設定，裡面還可以自己定義很多http裡頭的設定，像是add_headers, authenticate, config, set_cookies, timeout, user_agent, verbose都是。

set_proxy <- use_proxy("64.251.21.73", 8080, username = NULL, password = NULL,
  auth = "basic")
GET("http://had.co.nz", use_proxy("64.251.21.73", 8080))
GET("http://had.co.nz", set_proxy)
#url 放proxy的ip，為字串
#port 放埠口，為數字（容易弄錯）

而基本上，rvest是使用httr為基底的，所以也可以使用，這邊是read_html的範例代碼，可以看出其關係

# From a url:
google <- read_html("http://google.com", encoding = "ISO-8859-1")

# From an httr request
google2 <- read_html(httr::GET("http://google.com"))
# 所以可直接將use_proxy用在rvest裡面，雖然rvest的doc文件沒有提得很仔細，但由其原理可以推出來
google <- read_html("http://google.com",use_proxy("64.251.21.73", 8080))
#同樣的，html_session可以如法泡製
google <- html_session("http://google.com", use_proxy("64.251.21.73", 8080))

那加入proxy後，最常見的問題便時連線時間變很長，容易不穩定，此時有幾個解決方法。
使用error handler並且調整time out的設定。
1. 調整time out設定
在httr中可以設定timeout函數，設定等待回傳的時間值
2. error handler
可以使用purrr::safely或是base::tryCatch來解決拋出問題的時候，該如何處理，可以設定一個轉換proxy ip的設定，但某個proxy不穩定的時候，置換另一個，或是直接不使用proxy本地來爬取

提升VIM使用技巧：移動+搜尋/取代+複製貼上

發表留言

雖然在一些流行的程式語言，通常都有各自專門的IDE(整合開發環境)，但避免不了的便是在各種unix-like環境中闖蕩，而VIM是這些環境中幾乎都有的編輯器，所以多少都會有一些使用經驗，比如insert模式、normal模式等等，但其實只要再加上多一點點的快捷鍵和使用習慣，可以再讓自己的開發能力上升。這邊推薦在盡量練習關於“移動”的快捷鍵，如何搜尋和替代，以及直接不用滑鼠的複製和貼上、剪下。這三個類型的任務都有頗不錯的參考可以閱讀。(基本上，每個鍵盤的按鈕在vim中都有意義)

移動

搜尋和替代

複製、貼上、剪下

下面的書和資源蠻不錯的：

screenshot.pngscreenshot.png

R資料處理小技巧：使用dplyr的filter加上stringr的str_detect，可以用來過濾複雜文字資料

發表留言

有時候在處理一些資料，其特定欄位的描述是很文字的時候，此時可以將dplyr中的filter搭配上stringr包裡頭的str_detect來處理，可以發揮很好的效果，這邊舉例來說明。

下面這筆資料，假如我想要知道有幾個人他的姓名中含有“秀”的

screenshot.png

那麼這類問題就可以用這種搭配來解決，可看下面的代碼

library(dplyr)
library(stringr)
data %>% filter(str_detect(`<code>人名`,"秀"</code>)) %>% summarise(number=n())