月份: 十一月 2016

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

install.packages("fortunes)

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

#隨機丟出一篇格言
fortunes()
#丟出第100號格言
fortunes(100)

 裡面有的對話蠻“幽默的”
 

近年來打者個人化醫療旗幟的各式各樣研究大量湧出，每個研究者切入的角度都很不同，這邊分享由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到次世代定去，還談了很多不同相關領域的應用。

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

使用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

ssh   accountname@ip -X

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

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

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

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

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

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

p是每個edge存在的機率

而不論是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

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

緊接者在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中，主要有三個函數帶有“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的位置，注意兩個都只會傳出“第一個”符合predict function的值！

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

以下摘錄自此篇翻譯：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醫學前沿》會得到您所需求的重要信息。

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

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()兩者的差異

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(.)))))