R fortunes package:R大師們的智慧語錄(冷笑話+kuso)

fortunes package是R裡面一個很有趣的包,裡面有379則各式各樣與R相關的格言和句子,形式多樣,有的是大師們的對話,有的就是很傳統的那種格言,可以將其設定在.Rprofile中,那樣每次打開R session就會看到名言,讓硬梆梆的分析腦,加入一點柔性的語料。

install.packages("fortunes)

使用上可以指定要那一篇,或是要隨機丟出來

#隨機丟出一篇格言
fortunes()
#丟出第100號格言
fortunes(100)
 裡面有的對話蠻“幽默的”
 screenshot.png

結合GWAS和RNAseq來實現個人化醫療:WHOLE approach(一)

近年來打者個人化醫療旗幟的各式各樣研究大量湧出,每個研究者切入的角度都很不同,這邊分享由Greg Gibson教授所努力的其中一種方式:WHOLE(Wellness and Health Omics linked to the Environment) approach。

Greg Gibson教授是Georgia Institute of Technology中Predictive Health Institute and School of Biology的負責人,本身一開始從果蠅的GWAS研究起家,近年來開始轉做人類醫學的研究,他有寫一本關於基因體學的教科書A Primer of Genome Science,輕鬆易入手,娓娓地談關於基因體研究歷年來的發展,從microarray到次世代定去,還談了很多不同相關領域的應用。
screenshot.png

Greg Gibson所提出的The WHOLE approach,其實就是將臨床上的危險因子評估與一套Genomic Risk Score相結合,這Genomic Risk Score則是整合GWAS研究找到的各式各樣SNP來進行計算,希望能提供病人更多的資訊。

 

screenshot.png

分別整理一下這位教授的研究,將其分類來做深入探討,尤其是他如何計算Genomic Risk Score的,這套方式是一種很貼近臨床的approach!
Transcriptomics for Translational Medicine

        Gibson, G., Marigorta, U. M., Ojagbeghru, E. R., & Park, S. (2015). PART of the WHOLE: A case study in wellness-oriented personalized medicine. Yale Journal of Biology and Medicine, 88(4), 397–406.


        Gibson, B. G. (2015). GTEx detects genetic effects. Science, 348(6235), 640–641. http://doi.org/10.1126/science.aab3002

        Gibson, G. (2014). Directions for the drivers. Nature, 512(7512), 31–3. http://doi.org/10.1038/512032a

         J., K., N., G., D.J., E., N.C., C., J.D., S., A.A., Q., … Gibson, G. (2014). Gene expression profiles associated with acute myocardial infarction and risk of cardiovascular death. Genome Medicine, 6(5), 1–13. http://doi.org/10.1186/gm560 (v)

        Hemani, G., Shakhbazov, K., Westra, H. H.-J., Esko, T., Henders, A. K., McRae, A. F., … Powell, J. E. (2014). Detection and replication of epistasis influencing transcription in humans. Nature, 508(7495), 249–253. http://doi.org/10.1038/nature13005

        Zhao, J., Akinsanmi, I., Arafat, D., Cradick, T. J., Lee, C. M., Banskota, S., … Gibson, G. (2016). A Burden of Rare Variants Associated with Extremes of Gene Expression in Human Peripheral Blood. American Journal of Human Genetics, 98(2), 299–309. http://doi.org/10.1016/j.ajhg.2015.12.023 (v)

        PRASHER, B. H. A. V. A. N. A., GIBSON, G. R. E. G., & MUKERJI, M. I. T. A. L. I. (2016). Genomic insights into ayurvedic and western approaches to personalized medicine. Journal of Genetics, 95(1), 209–228. http://doi.org/10.1007/s12041-015-0607-9

         Gibson, G., Powell, J. E., & Marigorta, U. M. (2015). Expression quantitative trait locus analysis for translational medicine. Genome Medicine, 7(1), 60. http://doi.org/10.1186/s13073-015-0186-7(v)

Omic Personality and the Genomics of Wellness
Tabassum, R., Sivadas, A., Agrawal, V., Tian, H., Arafat, D., & Gibson, G. (2015). Omic personality: implications of stable transcript and methylation profiles for personalized medicine. Genome Medicine, 7(1), 88. http://doi.org/10.1186/s13073-015-0209-4(v)
Wellness and Health Omics Linked to the Environment: The WHOLE Approach to Personalized Medicine(v)
Visscher, P. M., & Gibson, G. (2013). What if we had whole-genome sequence data for millions of individuals? Genome Medicine, 5(9), 80. http://doi.org/10.1186/gm484
Gibson, G., & Visscher, P. M. (2013). From personalized to public health genomics. Genome Medicine, 5, 60. http://doi.org/10.1186/gm464
Patel, C. J., Sivadas, A., Tabassum, R., Preeprem, T., Zhao, J., Arafat, D., … Gibson, G. (2013). Whole genome sequencing in support of wellness and health maintenance. Genome Medicine, 5(June), 58. http://doi.org/10.1186/gm462
(v)

使用OSX遠端ssh到server IGV發生錯誤

使用Broad  institute開發的Integrative Genomics Viewer(IGV)來看splicing graph或是bam檔算是很基本的技能,但在本機開啟IGV來讀取bam檔通常是痛苦的事情,所以使用server IGV遠端連線開啟就變得很必要,不然好幾GB的bam檔,傳來傳去或是光打開都是很吃電腦效能的。

這邊分享一個在mac上使用遠端ssh來登入執行server IGV會發生的問題:下面是error message

INFO [2016-11-28 16:10:38,623]  [DirectoryManager.java:163] [main]  IGV Directory: /home/weitinglin66/igv
INFO [2016-11-28 16:10:38,623] [DirectoryManager.java:163]  IGV Directory: /home/weitinglin66/igv
INFO [2016-11-28 16:10:38,632]  [Main.java:99] [main]  Startup  IGV Version 2.3.69 (98)03/08/2016 04:20 PM
INFO [2016-11-28 16:10:38,632]  [Main.java:100] [main]  Java 1.7.0_111
INFO [2016-11-28 16:10:38,633]  [DirectoryManager.java:72] [main]  Fetching user directory…
INFO [2016-11-28 16:10:38,764]  [Main.java:101] [main]  Default User Directory: /home/weitinglin66
INFO [2016-11-28 16:10:38,764]  [Main.java:102] [main]  OS: Linux
GConf 發生錯誤:無法連繫設定伺服器;可能是您需要啟用 ORBit 的 TCP/IP 網路連線,或由於上次系統當機導致有殘留的 NFS 鎖定檔。 請參考 http://projects.gnome.org/gconf/ 以取得更多資訊。(詳情 –  1: 不在使用的作業階段中執行)
INFO [2016-11-28 16:10:42,238]  [Main.java:150] [main]  A later version of IGV is available (2.3.84)
ERROR [2016-11-28 16:10:42,244]  [DefaultExceptionHandler.java:49] [main]  Unhandled exception
java.awt.HeadlessException:
No X11 DISPLAY variable was set, but this program performed an operation which requires it.
     at java.awt.GraphicsEnvironment.checkHeadless(GraphicsEnvironment.java:207)
     at java.awt.Window.<init>(Window.java:535)
     at java.awt.Frame.<init>(Frame.java:420)
     at java.awt.Frame.<init>(Frame.java:385)
     at javax.swing.JFrame.<init>(JFrame.java:174)
     at org.broad.igv.ui.Main.main(Main.java:81)
INFO [2016-11-28 16:10:42,248]  [ShutdownThread.java:51] [Thread-2]  Shutting down
這邊可能會有好多個“坑”!!
首先,要先講一些背景知識,IGV本身是一個java程式,我們將它在linux server中運行時,他其實是要開啟“視窗”,所以他其實需要使用linux 下的x window支援(又稱X11),然後這些顯示器的訊息其實是要串接到我們本機的MAC螢幕,所以我們osx系統也要能解讀那邊拋過來的訊息。所以可以先閱讀鳥哥這邊關於x window的介紹,另外,在osx系統有開發針對x window相容的軟件計劃XQuartz,用來將X11系統移植到OSX中。
那要如何解決這個error呢?
所以我們在ssh進入server時,要先開啟-X指令模式如下
ssh   accountname@ip -X
這樣的意思是讓server那邊要是在執行依賴X11相關的程式(如IGV, FastQC),可以將訊號從ssh丟過來對接,可以看下面從ssh manual截屏下來的
screenshot

瞭解生物內複雜網絡過程中Network Model的演進:Regular Network, Random Network, Small World Network, Scale-Free Network

“Models are all wrong, but some are useful"-George Edward Pelham Box
生物體內是由一個很複雜的系統所組成,藉由基因、轉錄組、蛋白質等各式各樣的調控而維持者生命。為了瞭解這樣複雜的系統,最好的方法就是去嘗試建立一個模型。所謂的模型本質上是用來解釋真實世界的資料而存在的,本質上是簡化一個複雜系統,且用數理等可量化的方式來描述,並且可以用這個模型來預測這個系統!

最一開始,關於網絡模型(network model)的探索,可以說是從基礎的數理分析上開始,從建立最基本的數理基礎,而Regular Networks是最簡單的一個網絡模型:

E = \frac{1}{2}nr

E代表的是Edge,n代表網絡中的節點數,r代表的是節點的degree。簡單說這樣的網絡每個節點的連結都是一模一樣的。可以參考下圖:

screenshot.png

這圖每個regular network都是由n=8, r=4所組成,其中1-lattice是這類regular network中的特殊形態,就是每個node間的連結模式都是一模一樣的!

很明顯的,這樣的regular network是建立在完美數學模型上所形成的,跟生物體內的網絡狀態很不一樣,應該不能用來解釋生物內的網絡。

random network的數學模型則是另一個極端的模型,在1950年代末期由Erdos和Renyi所提出,使用下面的方程式來描述這樣的網絡:

random \quad network \quad = \quad G(n,p) 
p是每個edge存在的機率

screenshot.png

而不論是random network或是regular network似乎都不太能完整的描述真實生物體網絡之現況,兩個洽好像是處在網絡模型的兩個極端,一個代表極端“完美”,一個代表極端“隨機”,在2000年前後,由Barabasi提出的scale-free networks和Strogatz提出的small world network可以更貼近地去形容生物內網絡模型的架構,剛好處在regular network和random network的中間。

Barabasi提出的scale-free network,其描述的網絡具有不同edge數的節點個數是成“幕次”的關係,簡單說,連接數多的節點只佔有非常少的數量,而大多數的節點都是跟這類hub的節點相接,這樣造成任兩個在網絡中的node最短距離其實是小的。

而Strogatz提出的small world network則是描述此網絡中任兩點的距離,是正比於總節點數的N數,簡單說,就是此網路節點數的增加,遠高過任兩點距離的增加,這類網絡的特性便是path length短,而clustering coefficient高。

就會長得如下:

第一個是small world network, 第二個是random network,第三個是scale-free network

screenshot.png

可以理解成少數的節點具有非常大量的連接數,而大多數的節點都不是,它們兩者間的數量級是呈指數關係的,好比現在富人和窮人在社會上的關係鏈結。

緊接者在2000-2010年中間,許多針對細胞內複雜網絡的分析文章開始如春筍出現,各種針對局部網絡motif和特定細胞內特定網絡的關係開始被探索!

參考閱讀:
random network經典論文
“On random graphs", Publ. Math. Debrecen, 1959, and “On the evolution of random graphs", Publ. Math. Inst. Hung. Acad. Sci, 1960.

Scale-free network經典論文
Barabási, Albert-László; Albert, Réka. (October 15, 1999). “Emergence of scaling in random networks". Science. 286 (5439): 509–512. arXiv:cond-mat/9910332

Small-world network經典論文
Watts, Duncan J.; Strogatz, Steven H. (June 1998). “Collective dynamics of ‘small-world’ networks". Nature. 393 (6684): 440–442. Bibcode:1998Natur.393..440W. doi:10.1038/30918. PMID 9623998. Papercore Summary http://www.papercore.org/Watts1998

Raval, A., & Ray, A. (2013). INTRODUCTION TO BIOLOGICAL NETWORKS.

使用purrr package學functional programming的觀念(四):search

在purrr package中,主要有三個函數帶有“search”的意涵,接受atomic vector和list兩種資料結構,再針對裡面的element作處理。這三個函數分別為contains(), detect(), detect_index()

#注意,丟進去的是list還是vector
#X為list
x <- list(1:10, 5, 9.9)
x %>% contains(1)
[1] FALSE
X %>% contains(1:10)
[1] TRUE
#Y為atomic vector
Y <- c(1,2,3)
Y %>% contains(1)
[1] TRUE
Y %>% contains(1:3)
[1] FALSE

#detect()和detect_index()兩者本質類似,都是使用predict function,只是一個回傳elements的值,一個回傳elements的位置

detect(20:40, function(x){
  x > 22 && x %% 2 == 0
})
[1] 24

detect_index(20:40, function(x){
  x > 22 && x %% 2 == 0
})
[1] 5

NEJM新英格蘭雜誌發佈新的中文平台:《NEJM医学前沿》

以下摘錄自此篇翻譯:Making An Impact on Clinical Practice and Research in China

有點震驚這個新聞,可以感受到華人醫療領域在中國經濟力量崛起後,會越來越重要!

兩週前,美國 NEJM集團與中國上海的嘉會醫學研究與教育集團JMRE 攜手啟動了一項戰略合作項目——《NEJM醫學前沿》。這一線上醫學資源由《NEJM醫學前沿》官網(www.nejmqianyan.cn)和移動客戶端組成。該項目匯集了發表於《新英格蘭醫學雜誌》(NEJM)和《新英格蘭醫學雜誌期刊薈萃》(NEJM Journal Watch)的 一批重量級論文的中文翻譯,今後每週將與《新英格蘭醫學雜誌》幾乎同步發布重要文章的譯文;同時,部分譯文附有由中國專家原創的相關述評。

近幾十年來,中國在醫學研究和患者醫療方面取得了長足的進步。但隨著經濟發展,中國人口的快速老齡化及生活方式的改變,導致了心血管疾病、腦卒中、糖尿病和癌症等慢病負擔日益增長。作為權威性科學期刊,NEJM 囊括了對中國醫療從業者有廣泛價值的資料。而這些內容只有便於讀者獲取,才能發揮最大的作用,廣泛服務於大眾;因此,採用其母語(中文)才能實現準確快速的傳播。 《NEJM醫學前沿》將翻譯來自於《新英格蘭醫學雜誌》和《新英格蘭醫學雜誌期刊薈萃》的高質量醫學和科學內容,並把高品質譯文帶給中國的醫務工作者、科研人員和教育工作者。我們致力於改進中國的臨床實踐,使中國融入國際醫學的主流。

《NEJM醫學前沿》創刊號中收錄了近五年來發表的眾多精品論文的譯文,其中150餘篇選自NEJM,100餘篇選自《NEJM期刊薈萃》。論文題材偏重於癌症、心腦血管疾病和糖尿病,其中近20篇文章附有由中國醫生科學家撰寫的述評。

《NEJM醫學前沿》計劃每週翻譯和出版發表於《新英格蘭醫學雜誌》和《新英格蘭醫學雜誌期刊薈萃》的新內容。我們的選題著重於上述三大疾病領域,今後會逐漸擴展到更為廣泛的醫學範疇。 《NEJM醫學前沿》還全文翻譯NEJM每週的目錄,並提供目錄的NEJM官網鏈接。讀者完成方便而安全的註冊後,即可免費獲取《NEJM醫學前沿》的內容,還可以通過NEJM官網獲取《NEJM醫學前沿》上的所有譯文鏈接。

我們的宗旨是通過提供可靠的醫學信息來促進全人類健康,相信您在《NEJM醫學前沿》會得到您所需求的重要信息。

目前看起來翻譯的品質還需要再提高!且似乎也是走付費訂閱路線!

使用AmiGO和QuickGO來獲取Gene Ontology資料

獲取由Gene Ontology Consortium提供的Gene ontology資料有很多種方式,主要可以分為:
1. 網頁介面
2. 直接下載
3. 使用API

這邊分享兩個主要的網頁介面可以用來獲取Gene Ontology的資料:AmiGO、QuickGO

AmiGO
為Gene Ontology Consortium官方支持的網頁工具,由NHGRI-NIH所贊助,所以會根據NHGRI-NIH所補助的方向來發展,在網頁裡可以直接query, browse, visualizing GO的資料和從MOD, UniProtKB及其他地方搜集而來的annotation,詳細的來源可以參考,具有wizard interface,最近已大幅度提升運行速度,增加許多其他的資料型態可供使用如display of annotation extensions, protein form(splice variants and proteins with post translational modification)。

QuickGO
由EMBL-EBI資助,為The Gene Ontology Annotation Project中的一部份項目在,他比較獨特的是提供filter capabilities for GO annotation, powerful integrated subset/slim interface, 和historical view of the term.

使用purrr package學functional programming的觀念(三):reduce

reduce()和reduce_right()可以將list或是vector中的element依序由左而右,或是由右到左的處理,提供一些基本的數理處置。

可以參考下面的範例來使用reduce()

1:3 %>% reduce(`+`)
#[1]6
1:10 %>% reduce(`*`)
#[1]3628800
5 %>% replicate(sample(10, 5), simplify = FALSE) %>% reduce(intersect)
#[1]7 #每次都不同,因為有sample這函數在裡面
x <- list(c(0, 1), c(2, 3), c(4, 5))
x %>% reduce(c)
#[1] 0 1 2 3 4 5
x %>% reduce_right(c)
#[1] 4 5 2 3 0 1
#可以用上面的例子來理解reduce()和reduce_right()兩者的差異

使用purrr package學functional programming的觀念(二):map family

map family是purrr package中很重要的一群函數,本質上用來取代R base中的apply, lapply , vapply, sapply等,基本功能就是一個可以使用另一個function作為參數的函數,而使用這參數所代表的函數來針對特定資料結構(主要是vector和list)中的element作處理,然後輸出成vector或是list的資料格式。

主要有:map_lgl(), map_chr(), map_dbl(),map_int(), map_if(), map_at(), map2_chr(), pmap()

map_lgl(), map_chr(), map_dbl(),map_int()這四組函數的差別,主要是固定輸出的資料格式,分別就輸出成logical、 numeric、 character 、 integer.

map_lgl(c(1,2,3), is_even)
#[1] FALSE TRUE FALSE
map_chr(c(1,2,3), is_even)
#[1] "FALSE" "TRUE" "FALSE"
map_dbl(c(1,2,3), sqrt)
#[1] 1.000 1.414 1.732
map_int(c(1,2,3), is_even)

map_at則是可以指定針對哪幾個element來使用參數中的函數來modify,輸出成一個list,而map_if則是先對每一個element使用一個邏輯判斷的函數(prediction function),再來看是否使用最後參數的函數來modify。

map_at(c(1,2,3,4), c(2), square)
# [[1]]
# [1] 1
# 
# [[2]]
# [1] 4
# 
# [[3]]
# [1] 3
# 
# [[4]]
# [1] 4
 map_if(c(1,2,3,4),function(x){x>2},square)
# [[1]]
# [1] 1
# 
# [[2]]
# [1] 2
# 
# [[3]]
# [1] 9
# 
# [[4]]
# [1] 16
[/code


map2則針對兩個以上不同的資料格式,一起使用參數中的function來處理。

map2(c(1,2,3),c("a","b","c"),paste)
# [[1]]
# [1] "1 a"
# 
# [[2]]
# [1] "2 b"
# 
# [[3]]
# [1] "3 c"

進階一點的使用方法是搭配dplyr的處理邏輯來使用,可以參考一下的寫法:

1:10 %>%
map(rnorm, n = 10) %>%
map_dbl(mean)
# Or use an anonymous function
1:10 %>%
map(function(x) rnorm(10, x))
# Or a formula
1:10 %>%
map(~ rnorm(10, .x))
# A more realistic example: split a data frame into pieces, fit a
# model to each piece, summarise and extract R^2
mtcars %>%
split(.$cyl) %>%
map(~ lm(mpg ~ wt, data = .x)) %>%
map(summary) %>%
map_dbl("r.squared")
# Use map_lgl(), map_dbl(), etc to reduce to a vector.
# * list
mtcars %>% map(sum)
# * vector
mtcars %>% map_dbl(sum)
# If each element of the output is a data frame, use
# map_df to row-bind them together:
mtcars %>%
split(.$cyl) %>%
map(~ lm(mpg ~ wt, data = .x)) %>%
map_df(~ as.data.frame(t(as.matrix(coef(.)))))

做GSEA分析的重要基因組(gene set)參考資料庫: The Molecular Signatures Databases (MSigDB)

在當下這個高通量測序的時代,找到幾百個有表現差異的基因後,通常就要進行一系列更複雜的分析,這些分析便是建立在用目前已知的生物醫學知識建立相關的調控路徑等,做更進一步的路徑分析或是網絡分析來看這些有表現差異的基因是用什麼方式在什麼pathway上對生物體造成影響所以找出一組(set)跟目標生物功能或路徑相關的基因,便是很重要的,而The Molecular Signatures Databases是Broad institute維護用來提供gene set作為Gene Set Enrichment Analysis(GSEA)的工具,這裡頭便是放了很多大家用來分析的基因set。

MSigDB資料庫跟常見的gene set資料庫(廣義來看)如GO, BioCarta, GeneMAPP不同的地方是裡面儲存的資料便是為了導入GSEA分析的格式,另外,相對於GO, BioCarta, GeneMAPP在註釋上較為嚴謹和封閉,MSigDB背後維護的團隊很樂意大家把自己研究領域的gene set寄過去跟他們分享,另外,MSigDB便沒有像上述資料庫有所謂的pathway diagram,畢竟他curated的方式主要是用crowdsourcing,且目的不同。

MSigDB將裡頭的gene set分成幾類:c1,c2,c3,c4,c5,c6,c7

分類     介紹   
  H: Hallmark Gene Set       為MSigDB團隊整合整個資料庫裡來自各個使用者所貢獻的Gene Set,使用一整套計算分析的流程,去除掉重複性,梳理比較具有代表性的founder sets     
C1: Positionl Gene Set       裡頭包含326組gene sets, 這邊的Gene Sets是以染色體上的位置來分類的,可以用來找出跟染色體上序列缺失、或是片段放大,表觀遺傳學上或是區域相關的效應
C2: Curated Gene Set     裡頭包含4729組gene sets,整合來自其他的資料庫如pubmet、reactome 、pathway database、BioCarta pathway database、KEGG gene sets等
C3: Motif Gene Set      裡頭包含836組gene sets,為基因cis-regulatory motif相關的基因資訊,包含跟promoter、3-UTR相關的資訊、transcription factor targets(主要由TRANSFAC資料庫而來)、microRNA targets  
C4: Computational Gene Set       包含858組gene sets,主要computational的方法從microarray的資料得出跟癌症相關的資訊,其中分成Cancer Gene Neighborhoods、Cancer Module兩類。 
C5: GO Gene Set       包含6166組gene sets,跟Gene Ontology內的分類一樣,分成BP、CC、MF 
C6: Oncogenic Signatures      包含189組gene sets,主要為跟癌症dis-regulated相關的細胞內路徑,由NCBI GEO和一些未公佈的實驗數據,以及文獻中已有的。
C7: Immunologic Signatures      包含4872組gene sets,主要在看不同細胞狀態下免疫系統相關的改變,由老鼠或是人類的實驗而來。    

可以先從Hallmark Gene Set裡頭來看,因為其是由所有其他分類整理而來,去除掉很多重複,所以使用Hallmark Gene Set當作起手,而假如針對特定很細的路徑,則再去細看那個collection中已有的各種gene sets。