分類: 未分類

最近對於所謂的“Annotation”這個辭彙在生物資訊裡的意思越來越有感覺，而且深深覺得其對於整個生資分析中所佔有極重要的角色，經驗高低可能就決定在此，很多基本的分析流程都是可以很快速學會的，但是對於自己資料做好的annotation卻是沒有固定模式可學的。

那什麼是annotation的定義呢？(這邊是在探討比較廣義層面的)

下面是annotation在wiki裡面的定義，簡單而精準：

An annotation is a metadata attached to text, image, or other data. Often, annotations make reference to a specific part of the original data.

好的annotation其實就是善用全世界科學家花費時間和金錢做實驗得出的寶貴資料，將其應用在自己的分析專案中，怎麼善用全世界上千個資料庫的資訊，就必須要有廣泛的知識，且另一方面，進行annotation的過程，本身就有一些“目的”導向，通常你就要很清楚自己想要了解的方向，比如你拿到了一個Gene list，你好奇的是這群基因的序列關係、功能上的關係、藥物反應、多少比例是文獻中跟癌症幹細胞相關的標記、多少是跟生長因素相關的、分別牽涉到什麼生物路徑的，根據不同的好奇方向，其實就要去找不同的資料庫來佐證，且怎麼呈現又是另一門需要花費心思的事情。

從bioconductor上面的package總數，有一半以上都是跟annotation功能相關的

Bioconductor的計劃其中很重要的便是讓R變成在做生物資訊annotation的專業工具，所以開發出很多的annotation package來達成這個目的，且有一套良好的規範，分別讓想要將資料庫包成package，或是如何好的query資料，都能更加的統一性，其中AnnotationDbi就是其中一個package，將各種可以用來註釋的資料包成package，不用再使用web request的方式，這結構便是由AnnotationDbi packages所決定的，而這其中分成很多種類：

1. Organism level: e.g org.Mm.eg.db
2. Platform level: e.g. hgu133plus2.db
3. Homology level: e.g. hom.Dmi.inp.db
4. Transcriptom level:e.g.TxDb.Hsapiens.UCSC.hg19.knownGene,EnsDb.Hsapiens.v75
5. System-biology level: GO.db

這結構很好的把所有在bioconductor中的annotation資源分成五大類，並能根據這五大類來找尋自己所需要用的annotation 資源，而不會被搞糊塗。

現在主要生物醫學資料庫如ensembl, nubi都有開放API接口，可以直接用程式去提出一堆的請求，不用像傳統的網頁查找那樣，省掉很多麻煩，只要稍微瞭解一下http原理，各種程式語言如R的httr（中文教學）或 是Python的requests，使用起來都算方便，

這裡遇到一個小問題，應該說想要更好的解決它：如何加速查找速度，假如我有上千個想要查找的資料？
舉例來說，當你想要找一大串4000個基因變異、相關基因功能，雖然不需要在網頁上點選4000次（可能花好幾天吧），幾行程式碼就可以解決，但當要丟4000次請求，來取得資料時，假設1次請求，要花5秒種取得資料，那麼4000次，還是要花6個小時，那要怎麼改善這個問題呢？因為我們很貪心的想要花更少時間（假如網路頻寬很大的話）

一種最直覺的方式，便是把要查找的variant列表，切成多份，然後同時跑，最後再把資料整合起來。

另外，因抱持者想要探索看看新方法，便嘗試使用Python中的Queue和Thread，用多線程的方式來解決這個問題。

我基本思路是將我要查找的Rsid(Variant的資訊)，用iterator的方式來吐出每次要query的id，而Queue object可以把我處理好跟我所開啟thread間資訊同步的問題。

這邊是最後使用的代碼

#這邊是我想要查找的基因變異variant的ID
test = RSID.iteritems()

#我希望開的線程數量
concurrent = 10

#因為資料是以json檔傳回，最後我都要塞回這個test_list中
test_list = {}

#使用queue的重點，是要解決共享變數不能被thread同時讀取的問題，他可以有效解決thread間在取用共享變數的問題
que = queue.Queue()
for i in range(concurrent):
    que.put(next(test)[1])

def doWork(*args):
    tmp = args[0]
    while tmp.qsize() > 0:
        rsid = tmp.get()
        test_list[rsid]=rsid
        #這邊是ensembl Rest API的接口
        server = "http://grch37.rest.ensembl.org"
        url_1 = "/variation/human/"
        url_2 = "?pops=1"
        ext = url_1+rsid+url_2
        r = requests.get(server+ext, headers={ "Content-Type" : "application/json"})

        if not r.ok:
            r.raise_for_status()
            sys.exit()
        decoded = r.json()
        test_list[rsid] = decoded

print("\t[Info] Queue size={0}".format(que.qsize()))

# 這邊是線程開啟和運行的地方
for i in range(concurrent):
    thread_name='Thd'+i
    t = threading.Thread(target=doWork, name=thread_name, args=(que,))
    t.daemon = True
    t.start() 

td = datetime.datetime.now() - st
print("\t[Info] Spending time={0}!".format(td))

P.S:後記，這種方式其實頗危險的，假如所使用的服務管理者，會特別管理高頻要資源的IP那麼就要注意，可能從此ip就被封了，所以最好可以掛載proxy在進行query，另一方面，時間允許下，就用單線慢慢抓吧！甚至必須去測試對方有沒有限制每秒能query的數量，就要去調適最佳化的參數！

閱讀參考：

Python Queue 模塊
what is the fastest way to send 100,000 HTTP requests in Python?
using iterators and generators in multi-threaded applications
iterables/iterators/generators
How many threads is too many?
Python模塊學習 – threading 多線程控制和處理

這本書The Biostar Handbook最近總是在Biostar論壇裡看到他的廣告，其實作者stván Albert很早就開始預告這本書的寫作計畫了，整個書的landing page很早就出來了，但最近才開始賣這本書，身為資訊焦慮患者，就敗了一本來研究一下。

Stvan Albert本身為Pennsylvania State University 生物資訊教授， Biostars: Bioinformatics Questions and Answers（就是Biostar論壇）就是他在維護和管理的，另外，Galaxy platform是他其中一個參與的專案，所以算是教學經驗豐富，本身做了許多有趣的生資項目。

整本書算是很入門的測序資料分析分享，可以看到部分章節不太完整，開頭談了一點生物資訊的入門起手觀念，目前三大主流測序機器pacficbio/illumina/miniON，帶一點unix指令教學，分別介紹基本的alignment 資料結構、vcf的資料結構、常用的分析方式和網站(GO analysis/BLASK)，基本上一個主題提三個例子，大多點到為主，Variant calling和RNAseq為兩大分析主題簡介，為起手的介紹。

整本書我個人喜歡的地方是他解說文章寫得很親民，不太會文鄒鄒的，且通常會附許多很經典的網路文章，來幫助他把一些“網友”的觀點原汁的加到他的書裡面（蠻好奇引用部落格文章的趨勢，會不會慢慢成主流），可以點選超連結，再進一步閱讀（買電子書的好處），所以買這本電子書，可以想成他幫你整理網路上很多經典文章的連結於其中，節省你上網大海撈針的時間，直接看這些經典文章，從如何做好專案資料整理到統計的p value該如何處理等，相對於很硬派的教科書，這本書應該是本“小書”，蠻適合作為定序資料分析的入門書。

問題其實很簡單，在使用panther database來幫手上的基因找尋他們的功能，通常一個基因都會對應到5-6個功能，而最後希望能用一張表格來輕鬆瞭解自己找到的基因功能，處理前的資料如下：

希望處理後，檔案可以看起來如下面這幫

實際上，這個問題使用aggregate就可以解決，輕鬆一行代碼的事情，裡面總共有三種寫法，供參考使用

#method 01
up.result.protein.function %>% aggregate(CLASS_TERM ~ ENTREZ, data=.,c)

#method 02</pre>
up.result.protein.function %>% aggregate(CLASS_TERM ~ ENTREZ, data=.,toString

#method 03
up.result.protein.function %>% aggregate(CLASS_TERM ~ ENTREZ, data=.,paste,collapse=",")
 

這問題也可以使用dplyr中的summarise來解決！

閱讀參考：
Aggregating by unique identifier and concatenating related values into a string

有時我們會需要讀入多個不同編碼或是代號儲存的，最麻煩的是如何產生這些檔案路徑的各種組合，此時便可以使用expand.grid函數，這函數可以產生各種組合的data.frame，然後使用do.call或是reduce，來將這些組合合併成路徑，另外，cross_n函數等，跟傳統expand.grid的功能是差不多，有興趣可以查一下他的函數介紹。

#method 1: with all dplyr tool
expand.grid(c("DNA_","RNA_"),c("sample-1","sample-2","sample-3"),c(".csv") %>% as.list %>% reduce(paste) %>% map(read_csv) %>% reduce(rbind)

#method 2: with do.call
expand.grid(c("DNA_","RNA_"),c("sample-1","sample-2","sample-3"),c(".csv")  %>% do.call(paste) %>% map(read_csv) %>% reduce(rbind)

其實這問題的解決方法有蠻多的，下面的討論串分別是用不同的概念和想法去實現的，而使用reduce和map這兩個的方式算是比較簡單，且可以輕易的串在分析流程的方式。
 

相關閱讀參考：
R – concatenate row-wise across specific columns of dataframe
concatenating rows of a data frame
How do I concatenate a vector of strings/character in R?
Finding all possible combinations of a three strings

rvest是由Hadley Wickham所開發的R的package，其結合httr和xml2，提供更簡潔方便的來抓取網頁的html, xml資料。單純使用httr其實就可以很方便地使用http的方式來抓取網站資源，但還需要再另外使用其他處理回傳資料的軟件包，所以rvest就解決這個問題，可以直接指定css或是xpath來抓取特定的xml/html nodes。

關於如何選定css或是xpath，可以使用Hadley Wickham推薦的selectorgadget，可以很方面的在網頁上選取想要抓取的部分。

這邊想要分享更進階的應用，使用rvest來爬取需要先登入帳號密碼的網站。

問題
主要是想用rvest來爬取藥物基因體學PharmGKB的資料庫，但此資料庫需要先登入

(此資料庫是可以申請License授權資料的，這邊只是想做一個練習！)

那要怎麼解決這個問題呢？

可以利用rvest裡面的html_session, html_form, submit_form, jump_to這幾個函數來解決！

# create a session
# 首先，第一個使用的url為有填寫登入帳號密碼的頁面
url.address       <-"https://www.pharmgkb.org/sign-in-start.do"  

## 使用html_session建立一個session(session其實是用package所建立的，可以保持你在此網站互動的足跡，想更深入可以稍微瞭解一下http協定和這篇文章)
pgsession <-html_session(url.address)
## 抓取填寫帳密欄位的表格（這邊的2是要根據每個人所登入的網頁不同調整，這邊帳密表格剛好在第2個）
pgform    <-html_form(pgsession)[[2]]       ## pull form from session

## 填寫帳密
filled_form <- set_values(pgform,
                          `userId` = "＊＊＊＊@gmail.com",
                          `password` = "＊＊＊＊＊")

## 送出帳密到這個session中，取得進一步爬取資料的授權
submit_form(pgsession,filled_form)

## 接者要跳轉入自己真正有興趣的網頁
rs.url.address <- paste0("https://www.pharmgkb.org/view/variant-annotations.do?variantId=",rsID)

# 使用jump_to函數在此session中送出url request（模擬個人瀏覽網頁的狀況）
web.page <- jump_to(pgsession, rs.url.address) %>%
                    read_html

參考閱讀：
關於http協定和session
1. how programming language implement “session"
stackover討論串
1. scrape password-protected website in R
2. automating the login to the uk data service website in R with Rcurl or httr
3. scraping password protected forum in r
4. Using rvest or httr to log in the non-standard forms on a webpage
rvest資料
1.rvest package on Github
2.rvest documentation on CRAN
3.rvest tutorial

學習LaTeX語言，且用其來寫書、文章、論文等結構性文體，一直都是一件很吸引我的挑戰，或許有人覺得用WYSIWYG工具（ex. Office）來寫便很方面，且現在Word支援數學公式的寫入。但仔細去學習後，會發現練習使用LaTeX其實能精進自己寫好“結構性”文章的能力，同時多體會一個語言背後的思維，哈不過也是一種自我約束的練習。

LaTeX是Donald Knuth所發明的排版程式語言，有完整的語法結構，Donald Knuth本身對於小細節是非常執著的電腦科學大師，而他所寫的TeX語言，基本這三十幾年改變不到，也就是一開始就寫得非常仔細，而學習一門新的程式語言最棒的地方在於去理解創造者背後賦予這語言的“哲學”，而LaTeX就是希望能讓書寫文字這件表達人類智慧的事情能更加“刻意”，你會因為要使用LaTeX，所以你會去注意關於文字呈現的眉眉角角。雖然需要花點時間學習，但算是為了讓自己更“organized”的刻意練習吧！另一方面，LaTeX中的很多語法和函數的用法，可以學到很多“古典”英語的詞彙，頗富趣味。

使用LaTeX寫好後的代碼，要經過轉譯才能形成最後的輸出，如pdf，因為我不想再多安裝一個專屬LaTeX的IDE，想盡量限縮我所使用的工具，後來發現Atom是一款不錯的工具，且搭配套件的話，基本上，可以完成編寫和轉譯的工作，這邊是目前整理起來的整個工作環境建置：

安裝LaTex：
-mac的用戶，安裝MacTex

安裝Atom：
點選連結

安裝Atom上的套件：

1. latex-plus
2. language-latex
3. autocomplete-bibtex
   4.linter-chktex
   5.spell-check
   6.latexer
   7.pdfview

這邊用LaTeX來排個有名的小詩，超級…自我感覺良好的排版，可以當作是一種自我癖好。當然，這個技能可以幫助你創建一些很有趣的履歷模塊，也比較能看懂別人的模塊，並且自己微調。

相關閱讀，有比較詳盡的介紹：Atom for LaTex Editing, Create Your First PDF With Atom and LaTeX

PharmGKB資料庫主要是儲存基因變異與藥物反應的資料，此資料庫在2000年發佈，由NIH NIGMS所資助的研究計畫，為NIH Pharmacogenomics Research Network的一部份，資料庫目前是由史丹佛大學的研究團隊所維護。

PharmGKB主要藉由下面八類工作來建立此藥物基因體資料褲

1. 建立基因變異和藥物與疾病的關係，藉由文獻回顧的方式
2. 整理在藥理上重要的基因，並且針對其基因上的變異、相關藥物和影響的路徑做彙整
3. 針對FDA藥物標誌上跟基因相關的資訊做整理
4. 提供作為檢驗藥物基因體學領域重要問題的資料
5. 根據所彙整的藥物基因體學資訊，來製作藥物劑量的臨床指引
6. 藉由建立合作的方式參與將藥物基因體學應用於臨床的研究項目
7. 發表PharmGKB的分析和註釋結果，主要是基因變異註釋、基於藥物基因體學的藥物劑量指引、
8. 將資訊公開於網站上，且提供下載

經由上面八類工作後，PharmGKB會將結果以下面六類形式在網站上提供瀏覽或是資料直接下載

1. Variant Annotations(Research-level annotations of individual publications describing the relationship between genetic variants and drugs; these are created on a paper-by-paper basis)
2. Drug-Centered Pathway
3. Very Important Pharmacogene Summaries
4. Clinical Annotations (Genotype-based pharmacogenomic relationships summarizing all variant annotations regarding the same genetic variant-drug association)
5. Pharmacogenomics-Based Drug-Dosing Guidelines
6. Drug Labels with Pharmacogenomic Information

PharmGKB的資料是來自於公開的研究文獻，使用自然語言處理和人工註釋這些文獻後，整理成知識資料庫（Konwledgebase），這些資料可以用來當作藥物與基因變異（或多型性）的支持證據，而為了有效將這些資料整理成對於臨床有影響力的指引，PharmGKB成立一個協會叫Clinical Pharmacogenetic Implementation Consortium(CPIP)，這協會主要負責擔任將PharmGKB中的資料整理成臨床指引，目的是幫助臨床可以有效地利用基因檢測的資料於藥物治療決策上。

那PharmaGKB裡面到底有哪些資料呢？主要分成六大類：VA(Variant Annotation), PW(Drug-Centered Pathway), VIP(Very Important Pharmacogene Summarries), CA(Clinical Annotations), DG(Pharmacogenomics-Based Drug-Dosing Guidelines), DL(Drug Labels with Pharmacogenomic Information)

這邊收集蠻多相關的資料，作為之後閱讀的素材：
1. CPIC的網站，裡頭資料豐富，甚至提供ppt素材，可以下載
2. PharmGKB網站，除了有原始六大類檔案外，還有一些guidelines可以下載
3. PharmGKB資料下載處