R處理小技巧:使用purrr來將data.frame轉成nested list

screenshot.png

其實這只是小小練習,想要將圖裡面原始的GO資料轉換成最後的list表單,這邊可以練習使用到split函數和map的subsetting功能,這兩個都是蠻有趣的工具,像是split本身其實是base函數,且新增了S4的method,另外其跟cut,tapply都被歸類在有“分組/分群”功能的函數。這邊還有使用到with,也是功能很強大的base函數,加在pipe之中,可以直接將data.frame切成vector,當然,with的功能其實是可以接受一個expression來。

第一步

all.up.GO.CLF_CLS.list <- all.up.GO.CLF_CLS %>% split(.$ONTOLOGY)

第二步

all.up.GO.BP.CLF_CLS   <- all.up.GO.CLF_CLS.list$BP %>% split(.$ENTREZID)

第三步

all.up.GO.BP.CLF_CLS   <- all.up.GO.BP.CLF_CLS %>% map(with, GO)

網絡視覺化分析工具Cytoscape簡介

screenshot.png    Cytoscape是一款開源免費的網絡視覺化和分析工具,專注在處理生物分子、蛋白質交互網絡,對於生物資料庫的整合頗佳,且其最有名的就是許多插件可以安裝,近年來在3.0版(2013)之後,Cytoscape軟件整個架構大改版,更模組化,且具有Restful API,對於一個從2001便開始的專案,真的是很努力在進步。

    在3.0版本之後,Cytoscape的核心開發團隊,並行了Cytoscape.js的專案,讓Cytoscape裡面視覺化的模塊可以有更多的影響力,而不侷限在生醫領域。

    這專案主要由著重的工作流有四個:Cytoscape Main(Core network analysks), Functional enrichment analysis, Network analysis of Proteomics data, Pathway analysis of gene expression data,且整個軟體的開發時程和想法非常公開透明,可以直接在網站頁面看到。核心的Cytoscape提供基本的資料整合、視覺化,更進階的功能他使用App的方式作為套件使用,而任何人都可以基於Cytoscape open API,用Java來撰寫這些App,許多生物資訊的論文便是將方法開發為Cytoscape App作為最後的應用。

      針對Cytoscape入門的書籍資料不多,最大宗的就是官方網站的資訊,另外,還有一本算蠻入門,但基本Cytoscape使用上的精髓都講的書可以用來參考:
Gang, Su.(2013)Instant Cytoscape Complex Network Analysis How-to. PackPub

Cover image for Instant Cytoscape Complex Network Analysis How-to

    整本書雖然截圖是以2.x版本的介面做展示,但每章節最後都會有可以連結到3.0相關資訊的網址,所以在瞭解上不會有太大的障礙,而書中主要是介紹Cytoscape核心功能為主。

    另外,關於學習Cytoscape使用方法,其實可以看一些他核心開發者Keiichiro Ono在LinkIn上的簡報,裡面有蠻多作為網絡視覺化分析工具開發者重視的點,和最近他蠻支持的做法,這邊其中三個我看完蠻受用的資料,分別是增進網絡分析重複性和更有效率的方法,以及基本建網絡的觀念。

下面是收集一些Cytoscape使用分享的網路文章,有興趣也可以看,但基本上覺得邊看Gang, Su.(2013)Instant Cytoscape Complex Network Analysis How-to. PackPub 上手比較快,尤其是已經稍微理解一些網絡資料格式,搭配Keiichiro Ono的簡報,這兩個看完功力就差不多惹!

由貼文時間排序:
2017
怎样用Cytoscape提升文章档次?
2016
一步一步,如魔鬼的步伐–教你製作簡單的network圖
cytoscape绘制网络图,导入数据的要求有哪些?
使用Cytoscape分析和可视化网络(1)
2010
R与Cytoscape配合在生物网络研究中的应用

論文閱讀:The Biology and Function of Fibroblasts in Cancer

Kalluri, R. (2016). The biology and function of fibroblasts in cancer. Nature Reviews Cancer, 16(9), 582–598. http://doi.org/10.1038/nrc.2016.73

screenshot.png

這篇文獻回顧是關於纖維母細胞(Fibroblast)在癌症中角色,近年來對於癌症細胞所處的微環境在癌細胞癌化和轉移中扮演越來越重要的角色,其中纖維母細胞便是要角。
撰寫者為Raghu Kalluri教授,目前是MD Anderson Cancer Center擔任癌症生物研究中心(Cancer Biology Department)和癌症轉移研究中的負責人(Metastasis Research Center),本身是MD/PhD,其研究關注癌症組織環境對於癌症細胞的惡化和轉移的影響。這裡有Raghu Kalluri教授再冷泉港會議中接受訪談的影片,裡頭有其對於其研究室近年研究方向的想法。

微環境之於癌症細胞的重要性
關於身體對於癌組織的反應,對於最終腫瘤組成很重要。這些癌組織微環境還被稱為Desmoplastic reaction, tumour stroma, tumour microenvironment(TME)。像是微環境中的免疫細胞、angiogenesis、oxygen tension、interstitual pressure, ECM remodeling和代謝產物的改變都是近年來認為跟癌細胞發展和轉移重要的因素。其中,這些癌組織的周圍細胞中,纖維母細胞是最為重要的一個角色。這些跟癌症有關的纖維母細胞被稱作Cancer-associated fibroblast。

此篇文獻回顧內容大致可以分成這幾個要點:
1. 關於纖維母細胞研究的起源
2. 激化的纖維母細胞和未激化的纖維母細胞的差異
3. 纖維母細胞跟生體內相關生理和病理上的角色:傷口癒合,組織纖維化
4. 跟癌症的連結:所謂的癌症組織的纖維化現象(Cancer Fibrobsis)
5. 激化纖維母細胞的多樣性(來源、組織特異性、狀態特異性)
6. 表觀遺傳調控在纖維母細胞的表型影響
7. 癌症相關的纖維母細胞
     a.  對於癌化現象的角色
     b. 在代謝中扮演的角色
     c. 在免疫中扮演的角色
     d. 對於抗藥性的關係

關於纖維母細胞的活化
通常纖維細胞在組織內的狀態是“未激活”(quiescent, resting),在組織受傷的時候才會被激活,而這些狀態還可以分為reversible 或是irreversible兩種,可能由epigenetic modification造成的。所謂的活化狀態的纖維母細胞最一開始的定義是由ACTA2活化。這個基因剛好是調控體內平滑肌的。而這類型的活化纖維母細胞為疤和慢性傷口的主要參與纖維母細胞。而纖維母細胞對於環境壓力的適應能力,讓其被視為在癌組受損時,能參與癌症復發的角色。

The host response to evolving cancer cells results in the generation of tumour tissue that contains components of normal organs. Such robust host responses define complex heterogeneity interactions of cancer tissue with host cells and are known as desmoplastic reaction, tumour stroma or tumour microenvironment.

纖維母細胞研究起源:
1858年在顯微鏡下,對於那些在微血管旁邊,沒有跟basement membrane緊貼,不是免疫細胞,也不是上皮細胞,可能是結締組織的一部份稱作纖維母細胞. 在一開始認為纖維母細胞在一般情況下,對於組織內的代謝、各種表現不太有影響力,唯有在特定情況下才會激活,如受傷後急性發炎和慢性發炎以及組織纖維化的過程。

比較精準的定義,應該稱纖維母細胞為resting mesenchymal cell,且可以被激活成mesenchymal stem cell.這些細胞在TGFbeta, PDGF, IL-6刺激下,可以快速proliferate. 所以其實大多狀況下,我們所分離出來的纖維母細胞其實是就是所謂“活化”的,他們可以synthesizing ECM, 分泌cytokines, chemokines, 招集immune cell, exerting physical forces to modify tissue architecture。這樣的特性,在一些研究中被認為會加速癌細胞的轉移能力,其中有研究顯示BM-MSC會加速breast-cancer cell的轉移能力。

纖維母細胞在病理狀況下的角色

急性的傷口癒合

“any injuries to the functional parenchyma results in a host response. Injurious stimuli(mechanical trauma, radiation-induced or extreme temperature-induced damage, toxins, pathogenes or metabolic impairments.)”

標準的傷口復原反應會招集免疫、發炎細胞來促進血管新生和ECM deposition. ECM中的type I, type III, type IV, and type V collagenes, laminins, fibronectin都是由被活化的纖維母細胞所產生的。另一方面,纖維母細胞也是這些ECM-degrading protease的來源(如Matrix metalloproteinase, underscoring their crucila role in maintaing ECM homestasis by reguulation of ECM turnover).當傷口癒合後,活化的纖維母細胞會逐漸轉成不活化的狀態。

組織纖維化
當對組織的傷害是持續性的,這類傷害可以是慢性物理性,毒素, 代謝改變或是自體免疫.這樣持續性的修復反應便會造成所謂的纖維化.這過程主要是由表觀遺傳學所調控,造成持續激活狀態的纖維母細胞,且活化 anti-apoptosis pathway and initiaitng proliferation. (51-54)

癌症纖維化(cancer fibrosis)
癌細胞就像是永遠不會癒合的傷口,所謂的cancer fibrosis又稱做desmoplastic reaction, tumour stroma and TME。目前已知道myofibroblast在癌症惡化和轉移有所關聯,但目前發現他的功能是很bimodal的。招集組織間纖維母細胞到癌症細胞周圍的,除了growth factor外,還有免疫細胞。在細胞損傷中,TGFbeta, PDGF, fibroblast growth factor 2(FGF2)都佔有很重要的角色。在很多癌症裡面,TGFbeta在招募活化的纖維母細胞都具有重要角色,且局部的癌症相關纖維母細胞(Cancer-associated Fibroblast)活化可能也是由在TME中的TGFbeta很有關係。

有研究顯示,在早期的癌化過程中,免疫細胞分泌的IL-1beta能活化纖維母細胞中的nNF-kB的訊息路徑,產生癌化前期的secretome。甚至有研究在具有癌症體質的成年公雞上,製造傷口,會使得傷口組織產生有侵襲性的癌症。在正常器官中,活化的纖維母細胞數量不多,但在癌症組織中,卻具有相當大量的活化纖維母細胞在組織之中,由組織裡上升的alphaSMA和FAP當作證據。在組織招募新的微血管和累積免疫細胞於癌症細胞的環境中,這些生長因子佔有非常重要的角色,如分泌VEGFA, PDGF, EGF,IL-6,IL-8。VEGF本身可以造成血管通透性增加,使得腫瘤富含ECM蛋白質,像是type I collagene, fibronectin等,可以刺激血管新生。隨者腫瘤長大,其會大量增加各種collagem, laminins, firbronectin, preoteoglycans, perostin and tenascin C。像是perostin/tenascin本來是不太會出現在正常現體組織,卻在乳癌中大量出現。換句話說,纖維母細胞催化癌化的過程,可能就是一種對於體內組織慢性傷口的反應。

纖維母細胞的功能和多樣性

現在已經知道有多個標誌可以用來觀察活化狀態的纖維母細胞,像是FSP1, S100A4, vimentin, alphaSMA, fap, pdgf RECEOTR-alpha, DDR-2。不同特異表現的纖維母細胞對於癌細胞的調控,也會有很大的差異,像是在乳癌研究中,缺少Caveolin 1表現的纖維母細胞,本身會對於癌症代謝上帶來調控(98,99),然後,高Caveolin 1表現的纖維母細胞也可以促進癌細胞的ECM modeling。作者也表示,這些活化纖維母細胞的特性可能也跟纖維母細胞原本組織來源相關,比如是從Marrow-derived precursors, endotheia cells, liver, pancreas stellate cells, resting tissue fibroblast。

纖維母細胞的表觀遺傳學調控

目前已知跟纖維母細胞調控癌細胞的相關訊息傳導路徑可能很多個(TCFbeta(18), 20),很多認為這些讓纖維母細胞活化,也可能跟表觀調控相關。

在Fibrobrosis-associated fibroblast(CAF)的研究中,發現使用5-azacytidine(一種demethylation藥物)可以使得FAFs的proliferate rate 下降,產生的collagen I 減少,alphaSMA減少。如活化JAK1,STAT3可已調控LIF,且同時抑制dnmt/jak活性,可以讓CAF轉換成具有invasive的能力。

癌症相關的纖維母細胞(CAFs)
在催生癌化上的角色
     本質上在癌組織中取得的纖維母細胞的特質就不太一樣,其展現的移動能力(migratory ability),自泌生長因子(autocrine growth factor),chemokine的能力都非常強。這些可能是表觀遺傳上的調控所造成的。下面分別是一些證據,在反應癌症相關纖維母細胞對於組織癌化所扮演的角色
1. CAFs可以催化non-invasion的癌細胞表現invasion的能力
2. CAFs藉由分泌SDF-1(Stromal cell-derived factor 1/CXCL12)來催生血管新生
3. CAFs上升的HSF1(heat shock factor 1)和癌細胞的HSF-1 drivedn tumorigenic programme相關
4. Yes-associated protein 1(YAP-1)會造成ECM stiffening
5. CAFs產生如MMPs(ECM degrading protease)讓癌細胞較高的motility/invasion ability
6. CAFs缺少分泌TIMPs,加強分析MMPs和DISINTEGRIN,adam10,這些跟刺激癌細胞的幹細胞特性相關
7. 在肺癌中,癌細胞分泌THBS2可以活化Fibroblast,如此造成轉移的能力

ECM remodeling by fibroblast may participate in the generation and maintenance of the cancer stem cell niche. Perostin(POSTN)可能就是其中的關鍵,對於Wnt signalling-mediated cancer stem cell niche maintenance. Elevated WNT signaling is found in colon cancer stem cells that are proximal to CAFs

p.s
這篇文獻回顧跟癌細胞會剝離原本的ECM造成轉移這主題相關

Buchheit, C. L., Weigel, K. J., & Schafer, Z. T. (2014). Cancer cell survival during detachment from the ECM : multiple barriers to tumour progression. Nature Publishing Group, 14(9), 632–641. http://doi.org/10.1038/nrc3789

Zachary T. Schafer 研究這個主題的教授(ECM, anoikis-resistance)

在癌症細胞相關代謝
在CAFs內的代謝情況就跟highly proliferating的細胞一樣,在能量產生上很大程度依賴glycolysis,導致這些代謝調控的改變,可能是由於如TGFbeta, PDGF, hypoxia, hypoxia-inducible factor 1 alpha, ROS-mediated suppression of CAV1所造成的。這樣的現象(Warburg effect)同時伴隨者跟增加的catabolic activity/autophagy有關(TGF-beta signaling,IDH3 down-regulation)。CAF可能扮演重要的anabolic/catabolic balance in the regionally diverse milieu(環境)。Tryptophan/arginine starvation對於T cell活化和淋巴球的功能具有重要角色,CAFs藉由調控代謝環境,藉此對於免疫細胞的反應做出調控。

這篇論文在談論為何這些高度proliferate cell換依賴glycolysis呢?

Heiden, M. G. Vander, Cantley, L. C., Thompson, C. B., Mammalian, P., Exhibit, C., & Metabolism, A. (2009). Understanding the Warburg Effect : Cell Proliferation. Science, 324(May), 1029–1034. http://doi.org/10.1126/science.1160809

這篇論文在談論癌症細胞在能量代謝平衡的問題

Fiaschi, T., Marini, A., Giannoni, E., Taddei, M. L., Gandellini, P., De Donatis, A., … Chiarugi, P. (2012). Reciprocal metabolic reprogramming through lactate shuttle coordinately influences tumor-stroma interplay. Cancer Research, 72(19), 5130–5140. http://doi.org/10.1158/0008-5472.CAN-12-1949

這篇文章有講到CAF/CSC裡 MCT1 increased, GLUT1 decreased

Oncogenes induce the cancer-associated fibroblast phenotype: Metabolic symbiosis and “fibroblast addiction” are new therapeutic targets for drug discovery
Catabolic fibroblasts donate the necessary fuels (such as L-lactate, ketones, glutamine, other amino acids, and fatty acids) to anabolic cancer cells, to metabolize via their TCA cycle and oxidative phosphorylation (OXPHOS).Thus, oncogene activation (RAS, NFkB, TGF-β) and/or tumor suppressor loss (BRCA1) have similar functional effects on adjacent stromal fibroblasts, initiating “metabolic symbiosis” and the cancer-associated fibroblast phenotype.Targeting “fibroblast addiction” in primary and metastatic tumor cells may expose a critical Achilles’ heel, leading to disease regression in both sporadic and familial cancers.”

在癌症組織中纖維母細胞與免疫的關聯

CAF對於免疫功能的調控(pleiotropic immunomodulatory functions)主要是由他的secretome所主導的,分泌ECM protein, ECM-remodeling enzymes, production of a plethora of cytokines, chemokines. 目前的研究主要是展現了fibroblast secretome的heterogeneity,its differential impact on tumour biologt.未來可能需要更多fibroblast-specific deletion cytokine/chemokines in preclinical tumour models的實驗,可能才可以提供比較精準的模型。

一般來說,認為CAFs可以提供immunosuppressive TME.但整體還是要看在哪類的TME中。目前認為跟CAFs相關的secretome裡有IL-6, IL-4, IL-8, IL-10, TNF, TGFbeta, C-C motif chemokine ligand 2(CCL2), CCL5(RANTES),CXCL9, CXCL10,SDF1, PGE2, NO, HGF, Human Leukocyte Antigen G(HLAG)。

IL-6 signaling被認為可能會抑制DCs的成熟,讓T cell無法活化,造成T cell anergy,另外,IL-6 可以使得monocyte differentiate into macrophage lineage和活化mast cells.但還不清楚是單純IL-6就可以,還是要其他的cytokine參與,其中TGFbeta也被認為在其中扮演重要的角色,可能和其調控Thelp 17 cell相關。

CXCL14從CAF而來,也被認為可以造成macrophage recrute到腫瘤。
CCL2 in breast cancer 則認為跟這件事有關,且其同時會被macrophage,DCs所分泌。
T cell recruitment所需的cyokine在CAF的secretome也可以找到(CXCL9, CXCL10, SDF-1),同時像是negative co-regulatory immune signal(PDL-1/2)跟抑制T cell activation有關。

targeting FAP+ CAF顯示出anitumour的效果,其藉由intramoural recruitmnet of CD+8 T cell SC04 cd+8 cell-mediated cancer cell killing

癌症相關纖維母細胞和抗藥性的關係

“mechanisms of resistance involving the stroma include the modulation of pathways involving the stroma include the modulation of pathways involving cacner cell-ECM/ CAF-ECM adhesion, cytokine, chemokine…。"

下面這三篇都在談論TEM對於藥物耐受性的關係

Correia, A. L., & Bissell, M. J. (2012). The tumor microenvironment is a dominant force in multidrug resistance. Drug Resistance Updates, 15(1–2), 39–49. http://doi.org/10.1016/j.drup.2012.01.006

 

Dittmer, J., & Leyh, B. (2015). The impact of tumor stroma on drug response in breast cancer. Seminars in Cancer Biology, 31, 3–15. http://doi.org/10.1016/j.semcancer.2014.05.006

 

Ohlund, D., Elyada, E., & Tuveson, D. (2014). Fibroblast heterogeneity in the cancer wound. Journal of Experimental Medicine, 211(8), 1503–1523. http://doi.org/10.1084/jem.20140692

CAF-mediated immune modulation, pro-angiogenic actions and metabolic reprogramming of the TME might aid in cancer cell survival and facilitate escape from therapy-induced cancer control

R處理資料小技巧:使用Non-standard evaluation來加速資料整理(函數quote,parse,assign,eval,expression)

在使用R處理資料的時候,很多時候會遇到文字型態的資料,他有表達某種運算或是簡化的觀念,這邊舉一個例子,比如一比標註藥物基因體資料內容如下:“CYP2C191-10”,對我們來說,可以理解成CYP2C191到CYP2C1910,總共十筆資料,這時候該如何用R處理呢?此時便可以使用Non-standard evaluation的觀念來處理,下面是範例,這邊會使用到quote,parse,eval等在處理此類問題的函數,其實Hadley Wickham正在開發一個用來處理Non-standard evaluation的工具包Lazyeval

基本上,quote, parse, eval這三個最基本的函數便可以演示一下所謂的non-standard evalution在做什麼,簡單說就是當syntax非正常console可接受語法的時候,該如何讓其繼續發揮功能!

那什麼是expression呢?
我們每個在R command line下所寫的算式或是函數的使用,其實就是一種expression,尤其是在像R這類非常interactive的程式語言,當我們輸入進去R console時候便是讓R去將expression執行出來,也就是去evaluation expression。正常在console下打1+1,按下輸入後,便會跑出2,但假如我們的輸入是“1+1”的字串呢?那麼這就必須使用non-standard evaluation來處理這類型的問題。

比如quote和expression這兩個函數可以用來包著這樣的表示

quote(1+1)
expression(1+1)

單這樣的話,代表用來處理尚未被evaluation的expression,這是後便用eval函數來執行這些expression:

eval(quote(1+1))
eval(expression(1+1))

screenshot.png

你也可以使用parse來將字串轉化為expression,那麼便可以執行,比如
screenshot.png

parse(text = "1+1")
eval(parse(text = "1+1"))
#處理基因變異資料的轉換,使用到non-standard evalution的觀念
"*3-*8, *11-*16, *19-*21" %>% str_split(",") %>%
                              unlist() %>%
                              str_replace_all("\\*","") %>%
                              str_replace_all("-",":") %>%
                              map(~parse(text=.)) %>%
                              lapply(eval) %>%
                              unlist %>%
                              as.character() %>%
                              str_replace_all("^","CYP2D6\\*")

另外,在處理這種將“text”轉變成為變數這件事上,assignget這兩個函數非常好用,可以放在expression處理這部上層,來解決一些複雜變數處理的問題。assign和get兩個函數的功能剛好相反,可以參考下面的範例

assign("A",c(1,2,3)
#> A
#> [1] 1 2 3
assign("A",data.frame(a=c(1,2),b=c(1,23)))
#> A
#> a b
#  1 1
#  2 23
get("A")
#> a b
#  1 1
#  2 23

相關閱讀參考:
Advance R:Non-standard evaluation

Non-standard-evaluation and standard evaluation in dplyr

Struggling with Non Standard Evaluation

Fun Standardizing Non Standard Evaluation

I hate non-standarad evalution

用R做hash table來對照rsid和allele names

Hash table是程式語言中很重要的一種資料結構,可以方便我們用相對較省時和符合直覺的方式用來查找某些資料,在不同程式語言的實現方式都不太一樣,在R中我們其實有多種做法,都可以達到Hash table的效果,這邊舉一個例子,目的是為了設計用來對應基因變異與其變異的命名轉換,通常這些資料間的轉換出現的場景非常普遍,所以這謝的技術蠻實用的。

在R中可以分別使用vector, list和environment三種資料型態來達成這個目的。

# Study the R hash function

# 測試用來查找的資料
test <- allele_to_rs_table$`Allele Name`[sample(198,size = 10)]

# use Vector --------------------------------------------------------------

# _創建自己的hash table
vector.hash <- allele_to_rs_table$`Rs#`
names(vector.hash) <- allele_to_rs_table$`Allele Name`

# _直接把test放進去vector.hash來完成
vector.hash[test]

# use List ----------------------------------------------------------------
list.hash <- list()

#_方法一建立list.hash
lapply(allele_to_rs_table$`Allele Name`,function(x){list.hash[[x]] <- allele_to_rs_table$`Rs#`[ allele_to_rs_table$`Allele Name` == x ] })

#_方法二建立list.hash,使用傳統的方式
for (i in allele_to_rs_table$`Allele Name` ){
    list.hash[[i]] <- allele_to_rs_table$`Rs#`[ allele_to_rs_table$`Allele Name` == i ]
}

#_直接查找
list.hash[["CYP2C9*13"]]
lapply(test,function(x){list.hash[[x]]})
# use Environment ---------------------------------------------------------

#_創建新子環境
e <- new.env()

#_將環境中的object做處理
#_方法一
with(e, "CYP2C19*2" <- "rs4244285" )
#_方法二
assign("CYP2C19*2", "rs4244285", envir = e)

#_查找
base::get("CYP2C19*2", envir=e, inherits = TRUE)

當然這些是直接使用base R的函數來實現,當要對應的key多的時候可以考慮進階使用如hash/digest這兩個package來幫助提高速度。

閱讀參考:
1.Hash Table Performance in R: Part I
2.Hash Table Performance in R: Part II
3.Hash Table Performance in R: Part III

使用關聯性分析來探索臨床或基因表現資料

關聯性分析(association rules analysis)在商業上早已很普遍的被使用,尤其是零售業者用來發覺顧客對購買不同商品間有沒有特殊的習慣,其中購物籃分析(market basket analysis)是最廣為人知且基礎的一種關聯性分析方法。這種分析方式有一個很重要的關鍵,便是需要大量的樣本,否則便沒有太大意義,另一方面,這不是一種能產生所謂“統計顯著”的分析方式,更多的是提供一個探索性的工具,能不能從結果產生有意義的東西,還是需要依據“領域知識”。

在過去使用這種“購物籃分析“方式於醫療或是生醫資料會是非常”可笑的“,因為每個實驗或是臨床試驗都是精準設計樣本數,想問的問題都是非常明確,然後才去設計實驗或是收集資料,希望從結果得到一個”顯著的意義“,但隨者越來愈多基因資料或是臨床資料的出現,很多時候其實也是在做一種”大海撈針“,或許這種購物籃分析也能使用來”探索“這些生醫資料,來產生一些有趣的”關係“,再進一步設計實驗來回答這種關係!

那什麼是購物籃分析呢?

購物籃分析是頻繁模式探勘(Frequent Pattern Mining)中的一種實現方法,目的是為了找出在資料中某些頻繁出現的連結,而使用association rule是這方法的關鍵,所以有時候又有人會稱為關聯性分析(association rule analysis)。這方法出現的背景是當我們想要探討某些東西/事物/購買物品一起出現或是某些東西/事物/購買物品不會一起出現的關係。在商場上購物的例子中,就是探討顧客所購買的商品之間,是否會有一些特出的關係存在!

這方法的邏輯脈絡如下,我們會先嘗試找出每次購物中各種物件/物品間一起出現或是反之的關係,然後可以從這些關係組合出各種可能,然後可以用下面的方式表現,即所謂association rule,每種關係都代表一種可能的關聯規則:

computer=>antivirus\_software\[support=60%,confidence=60%\]

接下來,就是要有能評估一個關聯規則是不是所謂的frequent pattern,是不是我們所要找的,而這邊支持度(support)和置信度(confidence)這兩個指標邊能幫助我們理解個別項關聯規則對我們是否有意義,這兩個指標簡單來說分別代表這項關聯規則是否有用(usefulness)和可以確信(certainty)。所以通常我們會設下最下支持度闕值(minimum support threshold)和最小置信度闕值(minimum confidence threshold),當某項關聯規則能超過我們所設的闕值,才會成為我們所相信的關聯規則,而這闕值要怎麼設,則由領域知識來決定,並沒有所謂的標準,下面是者兩個指標的定義:

支持度(support)
support(A=>B)=P( A \cup B)
資料中item set A和iterm set B在每一筆記錄中同時存在的機率

置信度(confidence)
confidence(A=>B)=P(B|A)
在每一筆item set A存在的資料中,item set A和item set B同時出現的機率

有了上面的基礎架構,接者便是介紹將其實現的演算法,Apriori演算法算是最基本且常用的算法,R裡面的arules package便是使用這演算法

這邊放arules的使用代碼

library(arules)
#建立模擬的資料,裡面A和B這兩個item常常一起出現和消失
A <- c(rep(1,200),sample(c(1,0), size=800, replace = TRUE, prob = c(0.1,0.9)))
B <- c(rep(1,200),sample(c(1,0), size=800, replace = TRUE, prob = c(0.1,0.9)))
C <- c(sample(c(1,0), size=1000, replace = TRUE, prob = c(0.1,0.9)))
D <- c(sample(c(1,0), size=1000, replace = TRUE, prob = c(0.1,0.9)))
E <- c(sample(c(1,0), size=1000, replace = TRUE, prob = c(0.1,0.9)))
mimic <- data.frame(A,B,C,D,E)

#轉換成apriri函數可以接受的格式
trans <- as(data.matrix(mimic),"transactions")
rule.iter <- apriori(trans, parameter = list(support=0.2,
                                             minlen = 1,
                                             maxlen = 3,
                                             ext = FALSE,
                                             confidence=0.6,
                                             smax = 1,
                                             arem ="none",
                                             minval = 0.1,
                                             originalSupport = TRUE,
                                             maxtime = 5))
#看運算後的結果
inspect(rule.iter)

關聯性分析的閱讀參考:
1. 白話大數據與機器學習, 第十一章
2. 數據挖掘:你必須知道的32個經典案例,第五章/第七章

在臨床或是生醫領域的應用論文:
1. Distance-Enhanced Association Rules for Gene Expression, Proceeding BIOKDD’03 Proceedings of the 3rd International Conference on Data Mining in Bioinformatics, 2003
2. Mining gene expression databases for association rules, Bioinformatics, 2003
3. Integrated analysis of gene expression by association rules discovery, BMC Bioinformatics 2006
4.Efficient mining of association rules for the early diagnosis of Alzheimer’s disease, Physics in Medicine and Biology ,2011
5. icuARM-An ICU Clinical Decision Support System Using Association Rule Mining, IEEE Journal of Translational Engineering in Health and Medicine,2013
6. Association rule mining based study for identification of clinical parameters akin to occurrence of brain tumor, Bioinformation, 2013
7. Association rule mining to detect factors which contribute to heart disease in males and females, Expert Systems with Applications, 2013
8. Validation of an Association Rule Mining-Based Method to Infer Associations Between Medications and Problems, Applied Clinical Informatics, 2013
9. Mining Association Rules in Dengue Gene Sequence with Latent Periodicity, Computational Biology Journal, 2015
10. Dynamic association rules for gene expression data analysis, BMC Genomics, 2015
11. Prediction of Metabolic Pathway Involvement in Prokaryotic UniProtKB Data by Association Rule Mining, PloS One. 2016

將有興趣的基因列表使用org.Hs.eg.db查找最新的對應GO terms,再使用topGO做GO analysis

topGO package在bioconductor裡已經超過10年了,算是最早R裡面可以用來做GO term enrichment analysis的工具,但也因此裡面處理輸入資料的方式相對的來說頗不方便,語法比較接近R base,所以對於用習慣tidyverse體系的package的人來說會有點頭痛,這邊則來分享一下,如何使用topGO來做GO analysis

使用topGO做GO analysis主要有五個步驟:
1. 將基因列表轉換為Entrez ID
2. 使用org.Hs.eg.db資料庫分別找出每個Enrez ID所對應到的GO terms
3. 將資料結構處理成topGO中所使用的topGOdata object
4. 使用runTest進行GO analysis
5. 使用GenTable來獲取結果成data.frame

最麻煩的地方是理解topGOdata object建立所需要的常數,和每個輸入的結構!

這邊補充一下,各個下面代碼的重要資料結構長得樣子

screenshot.png

library(annotationDbi)
library(org.Hs.eg.db)
library(purrr)
library(dplyr)
library(topGO)

all.up.GO.Sphere_CLS.entrezID  <-  up.tukey.aov.Sphere_CLS %>% with(EntrezID)

#使用annotationDbi中提供的select方法,query的結果為data.frame的形式
all.up.GO.Sphere_CLS <- AnnotationDbi::select(org.Hs.eg.db,                                            keys = all.up.GO.Sphere_CLS.entrezID [!is.na(all.up.GO.Sphere_CLS.entrezID )] %>% as.character,
                                           columns = "GO",
                                           keytype = "ENTREZID")

#將entrezID對應GO term的data.frame格式轉換成list-like結構

all.up.GO.Sphere_CLS.list <- all.up.GO.Sphere_CLS %>%
                        filter(EVIDENCE != "IEA") %>% split(.$ONTOLOGY, .$ENTREZID)
all.up.GO.BP.Sphere_CLS   <- all.up.GO.Sphere_CLS.list$BP
all.up.GO.BP.Sphere_CLS   <- all.up.GO.BP.Sphere_CLS %>% split(.$ENTREZID) %>% map(with, GO)
all.up.GO.BP.Sphere_CLS.GO2geneID <- inverseList(all.up.GO.BP.Sphere_CLS)

#建立topGO所需要輸入參數allGenes, gene2GO兩個輸入變數geneList,all.up.GO.BP.Sphere_CLS
geneNames  <- names(all.up.GO.BP.Sphere_CLS)
geneList <- factor(as.integer(geneNames %in% (up.Sphere_CLS$EntrezID[!is.na(up.Sphere_CLS$EntrezID)] %>% as.character)))
names(geneList) <- geneNames

# construct topGOdata
GOdata <- new("topGOdata", ontology = "BP", allGenes = geneList,
              annot = annFUN.gene2GO, gene2GO = all.up.GO.BP.Sphere_CLS)

# calculate the result
resultFisher <- runTest(GOdata , algorithm = "classic", statistic = "fisher")
resultKS <- runTest(GOdata, algorithm = "classic", statistic = "ks")
resultKS.elim <- runTest(GOdata, algorithm = "elim", statistic = "ks")

# query the result
allRes  <- GenTable(GOdata, classicFisher= resultFisher,
                    classicKS = resultKS, elimKS = resultKS.elim,
                    orderBy = "classicFisher", ranksOf = "classicFisher", topNodes = 200)
allRes  <- allRes %>% transmute(GO.ID=GO.ID,
                                Term=Term,
                                Annotated=Annotated,
                                Significant=Significant,
                                Expected=Expected,
                                classicFisher=as.numeric(classicFisher),
                                classicKS=as.numeric(classicKS),
                                elimKS=as.numeric(elimKS))

allRes.up.Sphere_CLS <- allRes %>% filter(classicFisher < 0.05)