2017年12月9日 星期六

是否有可能取代鴉片類止痛藥的 NFEPP

目前常用的止痛藥很多是鴉片類(opioids),像是海洛因、嗎啡等等,它們作用在 opioid receptors (ORs)。Opioid receptors 是一種 GPCRs,分佈在中樞神經系統(central nervous system, CNS)和周邊神經系統(peripheral nervious system, PNS)裡的痛覺神經網絡(nociceptive neural circuitry),opioid receptors 被啟動後,會抑制神經細胞,中止痛覺的傳送,所以有麻醉止痛的效果。鴉片類藥物雖然止痛(analgesic)效果很好,但其副作用也不小,包括暈眩、噁心、呼吸窘迫( respiratory depression)和上癮等症狀。近來北美最為氾濫的此類藥物當屬 fentanyl,其效用比嗎啡高上八十倍,因為可以自己合成,比海洛因便宜,且只要一點點就可以就可以達到效果,常被混在海洛因或夜店藥裡,因此容易讓人不小心過量致死,加拿大 BC 省光今年上半年(一到七月)因 fentanyl 過量致死的就有九百多人

ORs 分為四種:mu (μ), kappa (κ), and delta (δ) opioid receptors 和 opioid receptor like-1 (ORL1),最常用來止痛的鴉片類藥物多是作用在 μ opioid receptors (MORs),但因為 MORs 也參與了腦內回饋系統(reward system),所以也會產生上癮的現象。

這篇研究的作者們認為,鴉片類藥物之所以會有其他副作用,是因為 ORs 分佈在中樞神經和周邊神經系統,當藥物作用在中樞神經系統的 ORs,就會出現嘔吐或呼吸困窘等症狀,而痛覺接收和傳送是發生在周邊神經系統,如果能讓藥物只作用在疼痛部位的神經細胞,也許便可以避免掉副作用。在慢性疾病造成的發炎和疼痛情況下,會有細胞組織酸化的現象。跟正常的細胞環境相比,受傷的組織其酸鹼值比較低,因此他們比較了在正常情況下和發炎情況下,藥物和 MORs 的反應,希望能找到一個只有在酸性環境下才會啟動 MORs 的藥物,以避免在正常情況下也會啟動 MORs 而產生包括上癮等等的副作用。以 fentanyl 為例,它的 pKa (disscoiation constant) 大於 8,因此在正常生理情況下(pH7.4)和發炎組織(pH5-7)的情況下都會被質子化(protonated),進而啟動 MORs。作者者們認為在正常情況下啟動 MORs 會產生不必要的副作用,所以想找一個 pKa 小於七,只有在組織發炎、需要止痛的情況下才會啟動 MORs 的物質。

通常要讓一個化學物質的 pKa 降低,可以把它的 H (hydrogen) 改成 F (fluorine),少了一個可以做 H-bond 的連結,就會降低兩者(ligand & MOR)之間的 interaction。他們改造了 fentanyl 的幾個地方後,選中了這個 pKa = 6.8 的 NFEPP ((±)-N-(3-fluoro-1-phenethylpiperidin-4-yl)-N-phenyl-propionamide),因為他們認為當 pKa 在 6~7 之間的話,組織發炎時的酸鹼值就足以質子化 NFEPP。


Figure: V Spahn et al, Science (2017)

跟 fentanyl 相比,NFEPP 在正常的生理情況下不會啟動 MORs,靜脈注射(i.v.)只會使發炎的部位依劑量的不同出現不同程度的麻醉、止痛(analgesia)效果,並且高劑量(75g/kg)的注射也沒產生呼吸困窘的副作用。它不像 fentanyl 在低劑量(8-12ug/kg)時即可在健康部位和發炎部位皆出現明顯的止痛效果,並且在其三分之一的量(32ug/kg)下就會造成致死的呼吸窘迫症狀。除此之外,fentanyl 在皮下注射(s.c.)為 30-50ug/kg 的劑量時會產生的依賴藥物的現象,然而 NFEPP 在高劑量(150ug/kg)下並沒有讓老鼠出現上癮的現象。Fentanyl 在 20-60ug/kg 會依劑量高低產生不同程度的四肢麻痺和便秘情形,NFEPP 則在高劑量(150ug/kg)下也未出現這些症狀。

看來降低 pKa 是能解決 fentanyl 產生的問題,不過他們可能需要多加一個情況做比較。在這個研究裡,試驗的老鼠都是有發炎症狀的,只是一隻腳有一隻腳沒有,然後兩隻腳做比較。保險起見,也許要加一個完全健康沒任何發炎症狀的老鼠做比較,看看是否在健康的老鼠身上也不會出現上癮等等的副作用。



Article:

TM Villanueva, Analgesia: Designing out opioid side effects. Nature Review Drug Discov (2017)

PRF / A New Opioid Targets Active Sites of Inflammation to Relieve Pain


Paper:

V Spahn & G Del Vecchio al, A nontoxic pain killer designed by modeling of pathological receptor conformations. Science (2017)










2017年12月8日 星期五

解酒藥可以殺死癌細胞

有意思哩!六十年來被用來當戒酒藥的 disulfiram(商品名 Antabuse / Antabus / 戒酒硫)也可以抗癌。在 1971 年的時候有篇臨床報告說一位 38 歲的乳癌病患因為癌細胞擴散到脊髓無法醫治而開始酗酒,醫生為了讓她戒酒便給她吃戒酒藥。十年後病患過世了,沒想到解剖後發現,她的腫瘤消失了,脊髓裡只剩一點癌細胞而已。在那之後,有幾個臨床研究和老鼠實驗也顯示 disulfiram 可以阻止癌細胞生長,雖然如此,disulfiram 並沒有因此受到很大個關注。

Drug Bank: disulfiram (tetraethylthiuram disulfide, DSF)

Disulfiram (DSF) 作用在 aldehyde dehydrogenase (ALDH),ALDH 是酒精代謝過程中需要的酵素(如圖),負責把 acetaldehyde 轉換成 acetate,disulfiram 抑制了 ALDH 後會造成 acetaldehyde 堆積在體內,會讓人感到極度的不適。



在這期的 Nature 中有篇由美國、捷克和丹麥合作的研究,檢視了在 2000 年到 2013 年間 240,000 位癌症患者和他們的用藥,發現在三千位有用 Antabuse 的患者中,有持續用藥的患者死亡率比後來停用的患者低惹 34%。這研究中檢視的不只是乳癌,還有攝護腺癌、直腸癌和其他種癌症,disulfiram 對它們的效果都是一樣的。此研究也顯示 disulfiram 可以減緩老鼠體內乳癌腫瘤的生長,並且和銅(copper, Cu)一起使用的話效果更好,喂有 disulfiram 的老鼠其腫瘤比沒有吃的小 57%,而和銅(DSF/Cu)一起吃的老鼠其腫瘤比沒吃的小 77%。他們也比較了這兩種飲食在老鼠體內代謝情形,disulfiram 在體內的主要代謝物 ditiocarb (diethyldithiocarbamate, DTC)會和銅結合成 CuET。吃了 DSF/Cu 的老鼠,CuET 在肝臟和血清(plasma)的含量比只吃 DSF 的高(感覺是廢話 XD),但有趣的是吃 CuET 的老鼠其腫瘤細胞裡的 CuET 量比其他器官高(只吃 DSF 的老鼠也有同樣的情形,但沒有高那麼多),而在細胞株上的檢測也顯示 CuET 的毒性(toxicity)比代謝物 DTC 本身還高。

他們發現 CuET 會和 NPL4 結合,NPL4 是 ATPase p97 (VCP) 的 adaptor protein,它通常會和另一個 cofactor, Ufd1, 一起和 p97 結合變成一個 complex (Npl4-Ufd1-p97)。p97 有很多不同的 cofactors,和不同的 cofactor 結合會執行不同的功能,其中主要的是把被 ubiquitin (Ub) 標記的蛋白質(Ub-proteins)從細胞內的細胞膜(i.e., ER membrane)中或是 protein complex 中抓出來,讓其能夠被降解(degradation),它參與的蛋白質降解管道有 proteasomal degradation 和 ERAD (ER-associated degradation),而和 NPL4 合作後的功用是負責 Ub-proteasome-mediated degradation,但是 CuET 和 NPL4 結合的話便會抑制它的功能。p97 和 NPL4 都是 cytosolic proteins,但當把 CuET 加到細胞上後,NPL4 無法執行功能,Ub-proteins 無法被降解而堆積在細胞中,引起後續的死亡反應。他們發現 CuET 會引發細胞的 UPR (unfolded protein response) 和 HSR (heat-shock response)。

所以目前可能的情況是 disulfiram 的代謝物 DTC 和銅結合後會抑制 p97-NPL4-Ufd1 參與的蛋白質降解功能,大量無法被降解的 Ub-proteins 在細胞裡堆積,引發癌細胞的死亡反應。因為 CuET 只會堆積在癌細胞腫瘤,所以這也可能是為什麼正常人吃 disulfiram 除了不適外沒有其他問題,只是 ..... 為什麼 CuET 只會堆積在癌細胞腫瘤中,原因還不清楚。

另外一個我覺得有趣的點是根據 Drug Bank 的資料,disulfiram 還有一個功能是:inhibits the peripheral benzodiazepine receptor.



Article:

Science / An old drug for alcoholism finds new life as cancer treatment


Paper:

Z Skrott et al, Alcohol-abuse drug disulfiram targets cancer via p97 segregase adaptor NPL4. Nature (2017)









2017年12月6日 星期三

關於偏頭痛(migraine)

(舊文轉過來)

你有偏頭痛嗎?這兩天 AAAS 在聖地牙哥有個年會,今天便在旗下的《Science》發了篇文章,告訴大家他們和 University of Miami Health System 的神經學家 Teshamae Monteith 討論了什麼關於偏頭痛的最新進展。偏頭痛不是普通的頭痛,它是一陣一陣難以忍受的痛,外加通常會噁心反胃,而且對聲音、光和味道很敏感。(嗯,這樣我好像沒偏頭痛過耶。)全世界大概有十分之一的人有偏頭痛,雖然目前無藥可治,但科學家們仍努力解開偏頭痛之謎,看能不能找出有效的治療方法。

目前我們對偏頭痛的了解是什麼?

這是個複雜的疾病,之前研究學者們都以為是跟血管相關的疾病,因為有的患者會在偏頭痛發作時感覺到血管規律的脈動。(阿,這我有過。)現在呢,科學家則認為應該是跟感官有關的疾病,因為很多的感官,包括光線、聲音、味覺、聽覺等等,都有受到影響。當偏頭痛發作時,患者常常無法專心、胃口和心情都會改變,也很難入眠。最引起 Teshamae Monteith 注意的是偏頭痛患者常被各種不同的徵狀困擾,像是在每波的陣痛之間對光的敏感度會增加,表示他們的神經或傳導系統各有不同和改變。大約三分之二的急性偏頭痛患者有異常性頭痛(allodynia),會對某種刺激特別敏感,例如沖澡時產生的水蒸氣都會令他們疼痛難當。不同的偏頭痛徵狀可能是因為每個人感官對刺激的敏感度不同。

目前偏頭痛的最新療法是什麼?

現在最主要的治療是一種叫 triptans 的藥物,它作用在血清素接收器(serotonin receptors),血清素被認為是參與偏頭痛的神經傳遞物質(neurotransmitter),血清素(5-HT)在平常的時候是低的,但偏頭痛發作的時候會升高。憂鬱症和偏頭痛也有強烈的關係,憂鬱症被指出和異常的血清素量有關,而憂鬱症患者較有可能有偏頭痛,偏頭痛患者也比較容易憂鬱。目前還不知道 triptans 是如何作用的,但它的確對某些偏頭痛患者有用,能夠止痛。它是個好藥,但可惜不是對每個人都有用,所以目前在研發各種新藥,希望是不像 triptans 這種會造成血管收縮的藥。

什麼是目前最有可能的新藥?

有個新藥是針對 CGRP (calcitonin gene-related peptide) ,它是一個被認為在急性偏頭痛發作時釋放出來的其中一種 peptide。在一個隨機、雙盲的預防性治療研究中,給患者不同劑量的 CGRP 藥物 telcagepant,雖然數據顯示效果不錯,但因為在一小部分有服用藥物的患者中,發現他們的某些肝臟酵素出現異常,以至於此藥物無法通過安全測試。不過現在有種抗體藥已進入第二階段(Phase II)的測試,因此被認為是目前最熱門的新東西,雖然第一階段的測試數量太小,無法確定它的效用有多大,但仍然是令人興奮的,因為它是第一個顯示安全的 CGRP 藥物。

偏頭痛遺傳的機率有多大?

關於致病基因的部分還需要更多的研究,因為偏頭痛這個病症太複雜,但目前大多數的患者都有家族史。Monteith 在 2008 年當住院醫師的時候,發現三個基因跟偏癱性偏頭痛(hemiplegic migraine)有關,這個病會使身體有一邊無力。而現在,有一堆基因被發現和偏頭痛有關,其中有些和 glutamate 有關。glutamate 是主要興奮性(excitatory)神經傳導物質,其功能是激化神經細胞,但因為腦部系統很多是用 glutamate 來傳遞訊息,所以要抑制它的功能有困難度。

偏頭痛的致病因子是什麼?

偏頭痛通常在青春期或青春期前發作,女性通常在生理期前後發作,男性的頭痛症狀也和他們的賀爾蒙變化有關。較低的社經地位、肥胖和不好的睡眠品質都會增加偏頭痛發作的頻率,腦部組織損傷(infarctions)是偏頭痛的併發症,其他像是抽菸和口服避孕藥則會增加偏頭痛的機率。

研究偏頭痛最困難的地方在哪?

我們(Teshamae Monteith)目前遇到的瓶頸是動物實驗的部分,因為牠們沒辦法和人腦完全一樣,也無法有極類似的偏頭痛症狀。人類的腦部影像技術有比較大的發展空間,在偏頭痛研究的應用上也比較讓人興奮。透過腦部影像,我們能夠更進一步了解跟偏頭痛有關的腦部構造、神經連結和化學物質。另外,研究曾經有過偏頭痛的已逝患者大腦,也可能對了解和偏頭痛有關的腦部損傷有所幫助,不過這種研究有其難解釋的部分,因為造成死亡的原因,例如中風、癌症等等,都會改變大腦狀態。(編按:如果大腦狀態因為真正的死因而改變,就無法知道哪些改變是真的跟偏頭痛相關,或是因為偏頭痛而造成的。)



原文:E Underwood. The Science of Migraines. Science News (Feb 2015)









2017年12月4日 星期一

R | ggplot: summary & error bar

這篇要講的應該是很常會用到的,就是加 error bar。

課堂講義請看這:R Graphics with Ggplot2 - Day 3

這邊我們要用到檔案資料 diamonds。

library(ggplot2)
library(dplyr)
data(diamonds)

我們先用 xtabs() 來看檔案中 color 的分佈情形。(xtabs() 的功能請看這篇。)

xtabs(~ color, diamonds)

Table 1
color
    D     E     F     G     H     I     J
 6775  9797  9542 11292  8304  5422  2808 

下面用 bar graph 看分佈,然後和上面做比對。

ggplot(diamonds, aes(color)) + geom_bar()

也可以分開設:(這個部分可以參考這篇會比較清楚)

d <- ggplot(diamonds, aes(color))
d + geom_bar()

#1 (color only)


接著,來看看每種 color 的價格,這裡可以用 by() function 來看每種鑽石顏色的價格總整理(summary)。

語法:by(response var, predictor var, summary)

這裡的 predictor variable (又稱 explanatory or independent variable)通常是個 factor (不是數字),會放在圖的 X-axis,response (or dependent) variable 指的是對應 factor 出現的反應,因此通常是放在圖的 Y-axis。

by(diamonds$price, diamonds$color, summary)

Table 2
diamonds$color: D
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    357     911    1838    3170    4214   18690
----------------------------------------------------------
diamonds$color: E
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    326     882    1739    3077    4003   18730
----------------------------------------------------------
diamonds$color: F
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    342     982    2344    3725    4868   18790
----------------------------------------------------------
diamonds$color: G
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    354     931    2242    3999    6048   18820
----------------------------------------------------------
diamonds$color: H
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    337     984    3460    4487    5980   18800
----------------------------------------------------------
diamonds$color: I
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    334    1120    3730    5092    7202   18820
----------------------------------------------------------
diamonds$color: J
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    335    1860    4234    5324    7695   18710 

也可以用 tapply() 功能來看。

tapply(diamonds$price, diamonds$color, summary)

$D
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    357     911    1838    3170    4214   18690

$E
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    326     882    1739    3077    4003   18730

$F
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    342     982    2344    3725    4868   18790

$G
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    354     931    2242    3999    6048   18820

$H
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    337     984    3460    4487    5980   18800

$I
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    334    1120    3730    5092    7202   18820

$J
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    335    1860    4234    5324    7695   18710 

如果只想看平均值,就把 summary 改成 mean。

tapply(diamonds$price, diamonds$color, mean)

       D        E        F        G        H        I        J
3169.954 3076.752 3724.886 3999.136 4486.669 5091.875 5323.818

在圖裡,可以用 summary_bin 這個功能來看,語法是:stat = "summary_bin"

或是:stat = "summary" (這兩個是一樣的)

summary_bin 裡面的功能有:fun.data, fun.y, fun.ymax (設最大值), fun.ymin (設最小值)

ggplot(diamonds, aes(color, price)) +
       geom_bar(stat = "summary_bin", fun.y = "mean")

在上面的語法裡,Y-axis 雖然已經設好是 price,但你需要在 geom_bar() 裡面的設定 Y-axis 是 price 的平均值(mean)或是中間值(median)等等,設定的語法是:fun.y = 

也可以用上面的已設定好的 d <- ggplot(diamonds, aes(color)) 把上面的語法改成下面這樣:(上面的語法是把 Y-axis 設在 ggplot() 裡面,但是 d 裡面只設了 X-axis,所以在下面的語法中,需要把 Y-axis 設在 geom_bar() 裡面。)

d + geom_bar(aes(y = price),
             stat = "summary_bin", fun.y = "mean")

#2 (color vs. price mean)


可以把圖和上面的 summary 表格做比較,看畫出來的是不是 mean。(和上面不同的是表一和圖一的指的是每種顏色的鑽石有幾個,所以圖 #1 的 Y-axis 是 count,而表二和圖二不只是每種顏色有幾個,而是每種鑽石顏色的價格,所以它的 Y-axis 指的不只是價格,還要指定是每種鑽石顏色的平均價格,還是每鑽石顏色其價格的中間值等等。)

下面來畫看看 median 的圖。

d + geom_bar(aes(y = price),
             stat = "summary_bin", fun.y = "median")

#3 (color vs. price median)


除了用 geom_bar(),我們也可以直接用 stat_summary_bin()stat_summary() 的功能。

stat_summary: summarise y values at distinct x values.

ggplot(diamonds, aes(color, price)) +
       stat_summary_bin(fun.y = "mean", geom = "bar")

或是用上面設好的 d。stat_summary_bin() 和 stat_summary() 是一樣的東西,上面和下面的這兩個語法畫出來的圖和上面用 geom_bar() 畫出來的 mean 那張圖是一樣的(圖 #2)。

d + stat_summary(aes(y = price),
                 fun.y = "mean", geom = "bar")

我們也可以畫點圖,用 geom_point() 就是下面那樣:

ggplot(diamonds, aes(color, price)) +
       geom_point(stat = "summary_bin", fun.y = "mean")

#4


上圖的 Y-axis 不是從 0 開始,而是從三千,我們可以把它設成從 0 開始。下面用上面設好的 d 和 stat_summary來畫圖,但基本上跟上面是一樣的,只是多了 ylim() 的語法。

ylim(min, max): Y-axis 的最小值和最大值

在下面的語法中,Y-axis 的最小值我們設為 0,最大值沒設,用 NA。

d + stat_summary(aes(y = price),
                 fun.y = mean, geom = "point") +
    ylim(0, NA)

#5


可以看到上圖的 Y-axis 是從 0 開始,而不是像圖 #4 那樣從三千開始。

接下來,試試在圖中加個 error bar,有四種方法。

注意,下面四種的預設都是 standard error (SE),除非你有指定另外的算法在裡面。

1. geom_pointrange(): vertical line with centre (中間點和其 error bar / standard error)

ggplot(diamonds, aes(color, price)) +
       geom_pointrange(stat = "summary_bin")



我們可以改變點和其 error bar 的大小和顏色。

ggplot(diamonds, aes(color, price)) +
       geom_pointrange(stat = "summary_bin",
                       size = 0.2, color = 'blue')

## No summary function supplied, defaulting to `mean_se()



summary_bin 裡面的 error bar 預設是 standard error (SE),你也可以把他設成別的,例如把它設成最大值 fun.ymax = max 和最小值 fun.ymin = min

ggplot(diamonds, aes(color, price)) +
       geom_pointrange(stat = "summary",
                       fun.y = mean,
                       fun.ymin = min,
                       fun.ymax = max)



可以把 pointrange() 和其他圖合在一起,這裡要注意一下,放在後面的會覆蓋住前面的,所以需要把 geom_point() 放在 geom_pointrange() 後面,如果相反的話會在圖裡看不出 geom_point() 裡面的設定。下面把兩個用不同顏色和不同大小的點區分。

ggplot(diamonds, aes(color, price)) +
       geom_pointrange(stat = "summary_bin",
                       size = 0.5, color = 'red') +
       geom_point(stat = "summary_bin", fun.y = "mean",
                  size = 1.5, color = 'blue')



Error bar 除了以上兩種外,也可以設其他的,例如讓它的範圍為中間 50%,最小值為第 25%,最大值為第 75%。

下面先用 group_by() 的功能把 color 分類,然後再用 summarize() 的功能在每種顏色裡算出其的價格。先把每種顏色的平均價格設為 y = mean(price),把誤差值(SE)的最小值指定為 ymin,設其為每種鑽石顏色的第 25%;把最大值指定為 ymax,設其為每種顏色的第 75%。

ps. quantile 分為:0 (0%), 0.25 (25%), 0.75 (75%), 1 (100%)

在下面先把算出來的指定為一個 data frame:df_color_price

df_color_price <- group_by(diamonds, color) %>%
  summarize(ymin = quantile(price, 0.25),
            ymax = quantile(price, 0.75),
            y = mean(price)),
            n = n())

# A tibble: 7 x 5
  color   ymin    ymax        y     n
  <ord>  <dbl>   <dbl>    <dbl> <int>
1     D  911.0 4213.50 3169.954  6775
2     E  882.0 4003.00 3076.752  9797
3     F  982.0 4868.25 3724.886  9542
4     G  931.0 6048.00 3999.136 11292
5     H  984.0 5980.25 4486.669  8304
6     I 1120.5 7201.75 5091.875  5422
7     J 1860.5 7695.00 5323.818  2808

上面語法中的 n() 是指在那個分類裡有幾個,在 R 裡面的解釋是:"The number of observations in the current group."(也就是說,上面得出的 n 的數字應該是要表 #1 裡的是一樣的。)

從上面的表可以看出 summarize() 算出每種鑽石顏色的平均價格(y)和前後 25% 的價格(ymin, ymax),還有每種顏色的鑽石有多少個(n),接著我們用這個 data frame 畫圖。

ggplot(df_color_price, aes(color, y,
                       ymin = ymin, ymax = ymax)) +
       geom_pointrange() +
       geom_point(aes(size = n))

在前面 df_color_price 的語法裡加入 n = n() 是要讓點圖裡的點會依據每種顏色的數量呈現出不同的大小,例如在前面 pointrange() 的例子裡加入 size =  的話可以指定點的大小。之前都是設數字(例如 size = 2),但是在這裡設 size = n 的話,因為在前面已經設了 n = n(),表示 n 是每種顏色的數量的話,那點的大小就會因數量而異。

也可以把上面兩個語法用 %>% 合在一起。

group_by(diamonds, color) %>%
         summarize(ymin = quantile(price, 0.25),
                   ymax = quantile(price, 0.75),
                   y = mean(price),
                   n = n()) %>%
ggplot(aes(color, y, ymin = ymin, ymax = ymax)) +
       geom_pointrange() +
       geom_point(aes(size = n))



也可以和 geom_bar() 合用。

ggplot(diamonds, aes(color, price)) +
       geom_bar(stat = "summary", fun.y = "mean") +
       geom_pointrange(stat = "summary",
                       size = 0.1, color = 'blue')



在上面的圖中,因為每種鑽石顏色的 standard error (SE) 很小,所以看不出來。

可以用上面已經設好的 df_color_price 來畫圖:

ggplot(df_color_price, aes(color, y,
                       ymin = ymin, ymax = ymax)) +
        geom_bar(stat = "summary") +
        geom_pointrange(size = 0.2, color = 'blue')

因為在 geom_pointrange() 裡面並沒有設 stat = "summary",所以圖呈現出來的會是設在第一層 ggplot() 裡的 error bar (ymin, ymax)。

也可以用 %>% 把兩個語法合在一起,變成下面這樣:(這裡沒用 n = n() 是因為並非畫點圖,所以不需要。)

group_by(diamonds, color) %>%
  summarize(ymin = quantile(price, 0.25),
            ymax = quantile(price, 0.75),
            y = mean(price)) %>%
ggplot(aes(color, y, ymin = ymin, ymax = ymax)) +
       geom_bar(stat = "summary") +
       geom_pointrange(size = 0.2, color = 'blue')



下面把 geom_pointrange() 和 geom_line() 合在一起用。

ggplot(diamonds, aes(color, price)) +
     geom_pointrange(stat = "summary_bin",
                     size = 0.2, color = 'steelblue') +
     geom_line(stat = "summary_bin", fun.y = "mean",
               group = 1, alpha = 1/2)



geom_line() 裡面的語法可以參考這篇:ggplot | Regression Line & Grouping

2. geom_linerange(): vertical line (只有一條直線 error bar / standard error)

如果單用的話就會像下圖那樣只有一條線,所以通常會和其他圖和用。

ggplot(diamonds, aes(color, price)) +
       geom_linerange(stat = "summary_bin")



可以和點圖和 bar graph 合用。和點圖一起用的話,其實畫出來的圖會跟只用 pointrange 差不多。你也可以幫 geom_linerange() 指定顏色,就是那條代表 SE 的直線,可以設成和點不一樣的顏色,如果直接用 pointrange() 的話就無法用不同的顏色。

ggplot(diamonds, aes(color, price)) +
       geom_point(stat = "summary_bin", fun.y = "mean",
                  size = 2, color = 'skyblue') +
       geom_linerange(stat = "summary_bin", color = 'blue')



同樣可以和 bar graph 合用。

ggplot(diamonds, aes(color, price)) +
       geom_bar(stat = "summary_bin", fun.y = "mean") +
       geom_linerange(stat = "summary_bin", color = 'red')



預設的 SE 跟之前 pointrange 的圖一樣太小看不出來,下面再用之前設的 df_color_price 再畫一次。

ggplot(df_color_price, aes(color, y,
                       ymin = ymin, ymax = ymax)) +
       geom_bar(stat = "summary") +
       geom_linerange(size = 0.5, color = 'steelblue')



geom_linerange() 裡面的 size =  是條線條的粗細。

3. geom_errorbar(): error bars.

如果只有 error bar 的話是這樣:

ggplot(diamonds, aes(color, price)) +
       geom_errorbar(stat = "summary_bin")



通常都是和 bar graph 合在一起用。

ggplot(diamonds, aes(color, price)) +
       geom_bar(stat = "summary", fun.y = "mean") +
       geom_errorbar(stat = "summary",
                     width = 0.1, color = 'grey')



geom_errorbar() 裡面的 width =  是用來設 error bar 的寬度。

下面用上面設好的 df_color_price 再畫一次。

ggplot(df_color_price, aes(color, y,
                           ymin = ymin, ymax = ymax)) +
       geom_bar(stat = "summary") +
       geom_errorbar(width = 0.1, color = 'steelblue')



4. geom_crossbar(): vertical bar with center.

終於來到最後一個。只用 geom_crossbar() 的話是這樣:

ggplot(diamonds, aes(color, price)) +
       geom_crossbar(stat = "summary_bin")



用上面的 df_color_price 畫圖:

ggplot(df_color_price, aes(color, y,
                           ymin = ymin, ymax = ymax)) +
       geom_crossbar(width = 0.2, color = 'steelblue')



接下來我想用檔案資料 mpg,用裡面的 class (X-axis) 和 cty (Y-axis) 畫圖。下面跟前面一樣先做總整理。

by(mpg$cty, mpg$class, summary)

mpg$class: 2seater
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
   15.0    15.0    15.0    15.4    16.0    16.0
----------------------------------------------------------
mpg$class: compact
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  15.00   18.00   20.00   20.13   21.00   33.00
----------------------------------------------------------
mpg$class: midsize
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  15.00   18.00   18.00   18.76   21.00   23.00
----------------------------------------------------------
mpg$class: minivan
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  11.00   15.50   16.00   15.82   17.00   18.00
---------------------------------------------------------
mpg$class: pickup
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
      9      11      13      13      14      17
----------------------------------------------------------
mpg$class: subcompact
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  14.00   17.00   19.00   20.37   23.50   35.00
----------------------------------------------------------
mpg$class: suv
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
   9.00   12.00   13.00   13.50   14.75   20.00 

把點圖和 crossbar 合用:

ggplot(mpg, aes(class, cty)) +
       geom_point(alpha = 0.5) +
       geom_crossbar(stat = "summary",
                     width = 0.2, color = 'red')



也可以像上面一樣,把 error bar 設為中間的 50%。

df_mpg <- group_by(mpg, class) %>%
  summarize(cty_ymin = quantile(cty, 0.25),
            cty_ymax = quantile(cty, 0.75),
            y = mean(cty))

ggplot(df_mpg, aes(class, y,
                   ymin = cty_ymin,
                   ymax = cty_ymax)) +
       geom_crossbar(width = 0.2, color = 'steelblue') +
       ylim(10, 35)



我把上面的 Y-axis 設了最大值和最小值 ylim(10,35),讓它和上面點和 crossbar 合圖的 scale 相同,比較好做比較。上面總整理的表格有 1st Qu. 和 3rd Qu.,Qu. 即是 quantile;1st Qu. 就是第一個 quantile,即第 25% ;3rd Qu. 是第三個 quantile,即第 75%。也就是說 1st Qu. 和 3rd Qu. 即是我們在 df_mpg 裡設的 cty_ymin 和 cty_ymax。可以把 crossbar 的最大值和最小值和表格裡的 1st. Qu 和 3rd Qu. 做比對,看看是否一樣。

下面同樣是把點圖和 crossbar 合用,但是用 stat_summary() 來畫。在 error bar 的部分另外設成 sd,用的是 package "Hmisc" 裡面的 mean_sdl

mean_sdl: smean.sdl computes the mean plus or minus a constant times the standard deviation

library(Hmisc)

ggplot(mpg, aes(class, cty)) +
       geom_point(alpha = 0.5) +
       stat_summary(aes(y = cty), fun.y = mean,
                    fun.data = mean_sdl,
                    width = 0.2, color = 'pink',
                    geom = "crossbar")



mean_cl_normal: computes 3 summary variables: the sample mean and lower and upper Gaussian confidence limits based on the t-distribution.

ggplot(mpg, aes(class, cty)) +
       geom_point(alpha = 0.5) +
       stat_summary(aes(y = cty), fun.y = mean,
                    fun.data = mean_cl_normal,
                    width = 0.2, color = 'pink',
                    geom = "crossbar")





最後,來練習一下。

比較 cut 和 depth 的關係,用點畫出中間值(median)和中間 50% 的 error bar。

group_by(diamonds, cut) %>%
  summarize(ymin = quantile(depth, 0.25),
            ymax = quantile(depth, 0.75),
            y = median(depth),
            n = n()) %>%
  ggplot(aes(cut, y, ymin = ymin, ymax = ymax)) +
  geom_pointrange(color = 'steelblue') +
  geom_point(aes(size = n), color = 'skyblue')




【 其他補充 】

有的地方可以看到 stat = "identity",這和 stat = "summary_bin" 有什麼不同呢?

stat = "identity" 通常是用在每個項目只有一個,例如在 diamonds 的檔案裡面,每個 color 的觀察只有一個,也就是說當你畫 price vs. color 的圖的時候,每個 color 裡不需要算價格的平均值或找它的中間值,因為就只有一個數(一個觀察)可以畫在 Y-axis。但是當資料變多的時候,也就是每種 color 裡面有好幾個,每個都是不同的價格的時候,那你畫 price vs. color 時就需要算出每種 color 的平均值,才能畫出 Y-axis,這時候就需要用 stat = "summary_bin" 來算出平均值或找出中間值。

另外,當你沒設 stat = " "  的時候,geom_bar() 裡的預設是 stat = "bin"

ggplot2: geom_bar() 網頁裡的解釋是這樣:

The heights of the bars commonly represent one of two things: either a count of cases in each group, or the values in a column of the data frame. By default, geom_bar uses stat="bin". This makes the height of each bar equal to the number of cases in each group, and it is incompatible with mapping values to the y aesthetic. If you want the heights of the bars to represent values in the data, use stat="identity" and map a value to the y aesthetic.



其他網頁參考:

My ggplot2 cheat sheet: Search by task

R document / ggplot2: geom_bar() function

ggplot 2 / Bar charts

ggplot2 / Vertical intervals: lines, crossbars & errorbars

ggplot2 / Summarise y values at unique/binned x










2017年11月27日 星期一

R | ggplot: 其他圖型變化 - 密度和重疊圖

這篇要介紹其他的圖型變化,例如用點的大小、深淺(透明度),或是兩種圖重疊來表現。

首先,安裝下面這些要用到的 library。

library(ggplot2)
library(dplyr)
library(hexbin)

用點的大小或深淺來表現是什麼意思呢?我們同樣用檔案資料 mpg 來舉例,假設我們要看 class 和 cty 的關係,先用 xtabs() 功能來看看。

xtabs() 是 cross tabulation,是用來看一個 variable 裡面的分佈,或是用來看兩個 variables 之間的關係。

如果想看某個 variable 裡面的分佈,語法是:xtabs(~ var1, data)

我們可以先來看看 class 裡面的分佈。

xtabs(~ class, mpg)

  class
2seater compact midsize minivan pickup subcompact  suv
      5      47      41      11     33         35   62

會顯示出 2seater 有五個,compact 有 47 個,minivan 有十一個等等。

如果想看兩個 variables 之間的關係的話,語法是這樣:xtabs(~ var1 + var2, data)

例如想看 cty 和 class 之間的關係,也就是 class 裡面各種車型裡 cty 的分佈情形。

xtabs(~ cty + class, mpg)



從上面的表格可以看到 2seater 的五個裡面,cty = 15 的有三個,cty = 16 的有兩個。

接下來,我們可以用點圖來表現。

ggplot(mpg, aes(class, cty)) + geom_point()



如果想把點變小,可以這樣:(不過在這個例子裡把點變小沒什麼特別的意義)

ggplot(mpg, aes(class, cty)) + geom_point(shape = '.')



可以把上面的表和圖做對照,例如在 2seater 裡面有兩種 cty, 15 的有三個,16 的有兩個。不過在上圖看不出來 15 的比 16 的多,這時我們可以讓點依數量而深淺(或透明度)不同,用的語法是在 geom_point() 裡面加入:alpah = 1/n(分母越大越透明)

ggplot(mpg, aes(class, cty)) + geom_point(alpha = 1/10)



除了用點的深淺或透明度呈現每個點的大小外,也可以用 geom_jitter()。當每個點有一個以上的量時,jitter 會用把它隨機分散。例如上面 2seater 裡的 15 有三個,在點圖上只呈現一個點,但在 jitter plot 上面會呈現出三個點。

ggplot(mpg, aes(class, cty)) + geom_jitter()

geom_jitter() 原本是下面這樣,後來簡化成上面那樣。

ggplot(mpg, aes(class, cty)) +
            geom_point(position = 'jitter')



可以看到 2seater 的 15 和 16 那附近有兩、三個點,但是因為預設值太分散了,所以分不清楚,這時我們可以設定其散落的範圍,例如把垂直移動(散落)的範圍設為 0,如此就能區分出 cty 的每個點,像是可以把 15 和 16 分開。

ggplot(mpg, aes(class, cty)) +
            geom_jitter(width = 0.2, height = 0)



也可以和 box plot 重疊。Box plot 的中間那一條是 medium (中間值,不是平均值),就是一個數列中,中間的那一個。上面那條線代表數列中的 25th percentiles (百分位),也就是中間值到最大值間的中間值,叫 Q3;下面那條線則是數列中的 75th percentiles,就是中間值到最小值間的中間值 ,叫 Q1。Q1 和 Q3 之間的距離,也就是 Q3-Q1,叫做 IQR (interquartile range)。(還是不懂的可以看這個影片

Box plot 的上下會有兩條鬚鬚,叫 whisker。它的長度在 R 裡的預設是 1.5,也就是:
max: Q3 + 1.5 * IQR
min: Q1 - 1.5 * IQR



我們先來看看 whisker 的預設是 1.5 是長怎樣(用預設的意思就是沒特別去設定),同時也把 outlier 設了顏色和形式,看看圖會變成怎樣。

ggplot(mpg, aes(class, cty)) +
            geom_boxplot(outlier.colour = "red",
                         outlier.shape = 1) +
            geom_jitter(width = 0.2, height = 0)



outlier.shape =  的語法裡: 0 = square, 1 = circle, 2 = triangle

設定 whisker 長度的語法是:coef = 

如果沒設的話就是 1.5 (如上圖)。如果設 coef = 0 的話,就是沒有那兩條鬚鬚。如果設為 coef = NULL,鬚鬚就會延伸到最大值和最小值。

ggplot(mpg, aes(class, cty)) +
            geom_boxplot(coef = NULL) +
            geom_jitter(width = 0.2, height = 0)



或是和 violin plot 重疊,之前因為覺得這種圖用不上,所以沒介紹,不過在這邊可以和 jitter plot 重疊看看。

ggplot(mpg, aes(class, cty)) +
            geom_violin() +
            geom_jitter(width = 0.2, height = 0)



再來,我們換特別的圖:geom_bin2b()

ggplot(mpg, aes(class, cty)) + geom_bin2d()



這個跟點圖的透明度類似,用顏色的深淺區分那個點的 count 有多少。

現在視覺話很夯,我們也可以用其他特別的圖來做視覺話,例如六角形的點。這個圖需要用到 library(hexbin)

在這個圖裡,我們把 class 裡面的 2seater 的資料拿掉,把拿掉後的資料設為 mpg_df。如果你還記得的話, !=  是指不等於,在這裡我們用 filter() 挑出 class 不是 2seater 的資料。

mpg_d1 <- filter(mpg, class != '2seater')

再用 mpg_df 來畫圖。

ggplot(mpg_d1, aes(class, cty)) + geom_hex()

也可以用 %>% 把上面兩個合在一起。(這邊看不懂的話,請看這篇。)

ggplot(mpg %>%
       filter(class != '2seater'), aes(class, cty)) +
       geom_hex()



最後,來做個練習。

Visualize the relationship between weight and price of the diamond in the “diamonds” dataset. Try the strategies we’ve just discussed to deal with overplotting.

ggplot(diamonds, aes(carat, price)) +
                 geom_point(alpha = 1/15)



也可以這樣畫。

ggplot(diamonds, aes(carat, price)) + geom_hex()












2017年11月19日 星期日

R | ggplot: Faceting

這篇要介紹的是 faceting,就是把檔案資料畫成一組圖,而不是只有一張圖。要看懂這篇需要先看前面三篇哦。

第一篇:Ggplot2 | Point plot & Box plot

第三篇:Ggplot2 | Regression Line & Grouping

上課講義請參考:R Graphics with Ggplot2 - Day 2

這邊同樣要用到下面這些:

library(ggplot2)

分圖要用的指令是:facet_grid(row ~ col)

或是:facet_wrap( ~ variable)

裡面放的是要依照哪個 variable 來分隔圖,前面是指用 row (列)來分圖,後面是指用 column (行)來分圖,可以只放其中一個。

我們先看下面的例子,一樣是用檔案資料 mpg。

ggplot(mpg, aes(displ, hwy, shape = drv)) +
            geom_point(aes(colour = factor(cyl)))



我們可以用 cyl = 5 的那個點來做 marker,5 只有兩個,是綠色三角形(drv, f),所以當你用 drv 來分圖的時候,只有 f 那個圖會有綠色的那兩點。

現在來用用看 facet_grid(),我們用 drv 並以 column 來分圖,所以會是這樣:

ggplot(mpg, aes(displ, hwy)) +
            geom_point(aes(colour = factor(cyl))) +
            facet_grid(. ~ drv)

前面沒放 variable 的話,就用點點 "." 。

也可以把 colour = factor(cyl) 這個放在前面的 ggplot() 裡面,變成這樣:

ggplot(mpg, aes(displ, hwy, colour = factor(cyl))) +
            geom_point() +
            facet_grid(. ~ drv)



注意到了嗎?綠色那兩點指出現在中間 f 那個圖裡。

如果是依 row 來分圖,就是這樣:

ggplot(mpg, aes(displ, hwy, colour = factor(cyl))) +
            geom_point() +
            facet_grid(drv ~ .)



接下來,我們可以依兩種 variables 分圖,行用 year,列用 drv。同樣用 cyl = 5 那兩點的做 marker,那兩點的年是 2008,所以它會出現在 f 和 2008 那一格圖中。

ggplot(mpg, aes(displ, hwy, colour = factor(cyl))) +
            geom_point() +
            facet_grid(year ~ drv)



上面用的 variable 中的元素都不會太多,頂多分個三、四個圖,但如果像是 manufacturer,就會變得很擠,像下圖那樣。

ggplot(mpg, aes(displ, hwy, colour = factor(cyl))) +
            geom_point() +
            facet_grid(manufacturer ~ .)



這時候可以用 facet_wrap(),會把圖由左往右,由上往下排出來。

ggplot(mpg, aes(displ, hwy, colour = factor(cyl))) +
            geom_point() +
            facet_wrap(~ manufacturer)



你也可以指定要有幾個 columns,可以用  ncol =

ggplot(mpg, aes(displ, hwy, colour = factor(cyl))) +
            geom_point() +
            facet_wrap(~ manufacturer, ncol = 5)



接下來,來點複雜的。我們可以不用整個 dataset,如果只想要其中幾種資料,或是要過濾掉幾個元素的話,可以用 filter() 的功能挑出你要用的。

先設一個 dataframe, d0。用 filter() 功能從 mpg 裡面挑出 cyl 不是 5,且 drv 裡面符合 "4" 和 "f" 的。

var %in% c(" ", " ") 是用來從你想要的 variables 那個 column 裡面找出符合條件的。%in% 意指 match。如果還記得的話, !=  是指不等於,可以參考基本指令這篇

d0 <- filter(mpg, cyl != 5, drv %in% c("4", "f"))



可以看見,d0 裡面的 cyl 沒有 5,drv 裡面也只有 f 和 4。

接下來,用這個畫圖。先設第一層 p0。(關於 %>% 的解釋和用法請看這篇。)

p0 <- d0 %>%
         ggplot(aes(cty, hwy))

再用 p0 設第二層,用法和解釋請看上篇

p4 <- p0 + geom_point()



上圖是用 d0 的資料畫出來的點圖 p4。

也可以依 variable 用顏色區分。

p5 <- p0 + geom_point(aes(colour = drv))



接下來,我們把圖 p5 用分成小圖。如果列的部分是用 drv 分圖,因為 drv 我們擷取出 4 和 f 的部分,所以只會有兩列,一列是 4,一列是 f。因為是也用 drv,所以可以明顯看出兩個顏色各一邊。

p5 + facet_grid(. ~ drv)



如果是用另一個 variable,例如分用 cyl 分圖,就會每個小圖都有兩種顏色。因為我們過濾掉 5,只剩下 4, 6, 8,所以只會有三列(columns)。

p5 + facet_grid(. ~ cyl)



因為已經把 p5 設好了,所以只在在後面直接加 facet_grid() 就好了。

p5 + facet_grid(drv ~ cyl)



也可以直接加在 p0 後面,像這樣:

p0 + geom_point(aes(colour = drv)) +
     facet_grid(drv ~ cyl)



如果是用 facet_grid() 裡面的另一個 variable,也就是 cyl,因為它不是 factor 或 character,出來的就會是色階,但也可以看到顏色上的區分。

p0 + geom_point(aes(colour = cyl)) +
     facet_grid(drv ~ cyl)



如果不是用已經用在 facet_grid() 裡面的 variable,而是用另外的 variable(例如 class),就會變成這樣:

p0 + geom_point(aes(colour = class)) +
     facet_grid(drv ~ cyl)



也可以同時用 drv 和 cyl 分列,變成這樣:(這裡有在點的部份指定顏色)

p0 + geom_point(color = 'navy') +
     facet_grid( ~ drv + cyl)



這樣很擠,可以用 facet_wrap()。

p0 + geom_point(aes(colour = cyl)) +
     facet_wrap( ~ drv + cyl)



如果把 cyl 設為 factor,就會變成這樣:

p0 + geom_point(aes(colour = factor(cyl))) +
     facet_wrap( ~ drv + cyl)



如果你有注意到的話,會發現分成小圖後,每個小圖的 scale 都是一樣的,這是預設好的,如果要讓每個圖有自己的 scale,可以在 facet 裡面加:scale = " "

讓 X-axis 有自己的 scale,用的是:scales = "free_x"

p5 + facet_grid(~cyl, scales = "free_x")



注意上圖中,每個小圖的 X-axis 刻度都不一樣,但是 Y-axis 的刻度是一致的。下圖則是設成讓 Y-axis 有自己的刻度。

p5 + facet_grid(~ cyl, scales = "free_y")



如果要讓兩個軸都有自己的 scale,就用:scales = "free"

p5 + facet_wrap(~ cyl, scales = "free")





最後,來練習一下吧。

Exercise. Use a histogram to explore the distribution of engine size for different vehicle classes.

ggplot(mpg, aes(displ)) +
            geom_histogram() +
            facet_wrap(~ class)





(天阿,這篇花了我好多時間,但是不快點打完怕我會忘記。)