不積跬步，無以至千里
本期我們嘗試復(fù)現(xiàn)2023年7月7日發(fā)表在Nature Communications上的Molecular profiling of aromatase inhibitor sensitive and resistant ER+HER2- postmenopausal breast cancers文章中的Fig3a塑径。
以下是原圖：
數(shù)據(jù)可以自行下載拾酝，也可評論區(qū)留言我私發(fā)給你卓嫂。
代碼及分析
之前分享代碼都沒有講代碼在干什么蹋辅，所以這次做一次嘗試進(jìn)行代碼分析。
-
加載包
library(dplyr)
library(ggplot2)
library(magrittr)
library(tidyr)
library(vroom)
library(forcats)
-
數(shù)據(jù)讀取
plotdata <- vroom(file = 'fig3a.csv', col_names = TRUE)
這里用到了vroom包的vroom()函數(shù)對數(shù)據(jù)進(jìn)行讀取缔俄，與傳統(tǒng)的read.table()函數(shù)相比疗绣，其主要有兩個優(yōu)勢：（1）讀取速度更快线召；（2）自動判斷文件分隔符。
查看一下數(shù)據(jù)：
head(plotdata)
# A tibble: 6 × 13
# Description `GeneRatio PRs vs GRs` `NES PRs vs GRs` `p.adjust PRs vs GRs`
# <chr> <dbl> <dbl> <dbl>
#1 HALLMARK_A… 70 2.78 8.33e-10
#2 HALLMARK_I… 59 2.67 8.33e-10
#3 HALLMARK_E… 45 -2.49 8.33e-10
#4 HALLMARK_E… 50 2.16 8.33e-10
#5 HALLMARK_T… 48 2.20 8.33e-10
#6 HALLMARK_I… 53 2.18 8.33e-10
-
數(shù)據(jù)清洗
這里我只保留了文章中的PATHWAY多矮，并且按照文章的順序進(jìn)行了排列缓淹，所以采用了原始粗暴的方式：
levels <- c(
'E2F TARGETS', 'G2M CHECKPOINT', 'MITOTIC SPINDLE', 'MYC TARGETS V1', 'MYC TARGETS V2', 'P53 PATHWAY',
'ESTROGEN RESPONSE EARLY', 'ESTROGEN RESPONSE LATE', 'IL2 STAT5 SIGNALING', 'KRAS SIGNALING DN', 'KRAS SIGNALING UP', 'MTORC1 SIGNALING', 'TNFA SIGNALING VIA NFKB',
'ALLOGRAFT REJECTION', 'COMPLEMENT', 'IL6 JAK STAT3 SIGNALING', 'INFLAMMATORY RESPONSE', 'INTERFERON ALPHA RESPONSE', 'INTERFERON GAMMA RESPONSE',
'HYPOXIA',
'EPITHELIAL MESENCHYMAL TRANSITION',
'GLYCOLYSIS',
'APICAL JUNCTION'
)
colors <- c(
rep('#279D77', 6),
rep('#CF6611', 7),
rep('#7974A1', 6),
'#E2348E',
'#E5B63D',
'#91793C',
'#686868'
)
這里的顏色是用FastStone軟件吸取的。
接下來工窍，我只對原文件plotdata進(jìn)行了過濾割卖，只保留了這些PATHWAY.
plotdata <- plotdata %>%
pivot_longer(cols = starts_with('NES'), names_to = 'Compare', values_to = 'NES') %>%
pivot_longer(cols = starts_with('GeneRatio'), names_to = 'GeneRatio_Group', values_to = 'Ratio') %>%
pivot_longer(cols = starts_with('p.adjust'), names_to = 'p.adjust_Group', values_to = 'p.adjust') %>%
mutate(Description = gsub(pattern = 'HALLMARK_', replacement = '', x = Description)) %>%
mutate(Description = gsub(pattern = '_', replacement = ' ', x = Description)) %>%
filter(Description %in% levels) %>%
mutate(Description = fct_relevel(Description, rev(levels))) %>%
mutate(Compare = gsub(pattern = 'NES ', replacement = '', x = Compare)) %>%
mutate(Compare = fct_relevel(Compare, c('PRs vs GRs', 'PRs ESR1 HIGH vs GRs', 'PRs ESR1 LOW vs GRs', 'PRs ESR1 LOW vs PRs ESR1 HIGH')))
這里有幾個知識點(diǎn)；
i. pivot_longer()用來改變數(shù)據(jù)框結(jié)構(gòu)患雏，對應(yīng)的還有pivot_wider()鹏溯，都是tidyr包里的函數(shù)，可以自主學(xué)習(xí)一下淹仑；
ii. 用gsub()函數(shù)對特定的字符進(jìn)行替換丙挽；
iii. 用fct_relevel()函數(shù)來對數(shù)據(jù)框中的列進(jìn)行因子化，來自于forcats包匀借，可以自主學(xué)習(xí)一下颜阐。
現(xiàn)在這個數(shù)據(jù)已經(jīng)變成這樣了：
head(plotdata)
# A tibble: 6 × 7
# Description Compare NES GeneRatio_Group Ratio p.adjust_Group p.adjust
# <fct> <fct> <dbl> <chr> <dbl> <chr> <dbl>
#1 ALLOGRAFT RE… PRs vs… 2.78 GeneRatio PRs … 70 p.adjust PRs … 8.33e-10
#2 ALLOGRAFT RE… PRs vs… 2.78 GeneRatio PRs … 70 p.adjust PRs … 1.67e- 9
#3 ALLOGRAFT RE… PRs vs… 2.78 GeneRatio PRs … 70 p.adjust PRs … 1 e- 9
#4 ALLOGRAFT RE… PRs vs… 2.78 GeneRatio PRs … 70 p.adjust PRs … 1 e- 9
#5 ALLOGRAFT RE… PRs vs… 2.78 GeneRatio PRs … 60 p.adjust PRs … 8.33e-10
#6 ALLOGRAFT RE… PRs vs… 2.78 GeneRatio PRs … 60 p.adjust PRs … 1.67e- 9
-
可視化
分析這個圖，有幾個關(guān)鍵點(diǎn)：
i. 首先這里肯定有分面吓肋，所以一定會用到facet_*函數(shù)凳怨；
ii. GeneRatio映射到了點(diǎn)的大小上；
iii. p value映射到了點(diǎn)的顏色上是鬼；
iiii. 點(diǎn)和x = 0之間有連線肤舞，這也就是我們常說的棒棒糖圖。
所以先簡單可視化一下：
p1 <- plotdata %>%
ggplot(aes(x = NES, y = Description)) +
geom_vline(xintercept = 0, color = 'grey', linewidth = 1) +
geom_segment(aes(x = 0, xend = NES, y = Description, yend = Description), color = 'grey', linewidth = 1) +
geom_point(aes(color = -log10(p.adjust), size = Ratio)) +
scale_color_gradient(low = 'blue', high = 'red', breaks = seq(3, 9, 2), limits = c(1, 9), name = 'Significance\n(-log10 FDR)') +
scale_size_continuous(name = 'GeneRatio(%)', range = c(2, 5), limits = c(30, 60)) +
scale_x_continuous(limits = c(-3, 3), breaks = seq(-3, 3, 1), expand = c(0, 0)) +
labs(x = 'Normalized Enrichment Score (NES)', y = NULL) +
theme_bw() +
theme(
panel.grid = element_blank(),
panel.border = element_rect(linewidth = 1.5, color = 'black'),
panel.spacing.x = unit(0.15, units = 'in'),
strip.background = element_rect(linewidth = 1.5, color = 'black'),
axis.ticks = element_line(color = 'black'),
axis.text.x = element_text(family = 'sans', colour = 'black', size = 11),
axis.text.y = element_text(family = 'sans', size = 11, face = 'bold', color = rev(colors)),
axis.title.x = element_text(family = 'sans', color = 'black', size = 14),
legend.title = element_text(family = 'sans', color = 'black', face = 'bold', size = 12),
legend.text = element_text(family = 'sans', color = 'black', face = 'bold', size = 11)
) +
facet_wrap(~Compare, ncol = 4)
到這一步：
解釋一下部分參數(shù)：
（1）geom_vline()添加豎線均蜜，對應(yīng)的geom_hline()添加橫線;
（2）geom_segment()添加棒棒糖的“棒棒”李剖；
（3）scale_color_gradient()與scale_size_continuous()更改映射詳細(xì)內(nèi)容，如顏色囤耳，大小范圍篙顺；
（4）scale_x_continuous()對x軸坐標(biāo)進(jìn)行個性化修改偶芍，包括范圍，斷點(diǎn)德玫；
（5）theme()中一系列的參數(shù)：panel改變圖形背景匪蟀，其中grid去掉了背景中的線條，border改變背景邊框宰僧，spacing.x改變圖形左右兩側(cè)空白的大小萄窜，防止分面圖形相隔太近；strip改變分面屬性撒桨；axis與legend改變坐標(biāo)、圖例文字等信息键兜。
截至目前，我們還有一個關(guān)鍵的問題沒有解決：原圖中有不同的方塊兒背景色，而我們還沒有坊秸，所以現(xiàn)在需要完成這個事情柑船，這也是這篇帖子最重要也是最難的一部分。
為了完成這個任務(wù)现诀，我選擇geom_rect()函數(shù)夷磕，但是這里有一個問題，上面的圖中縱坐標(biāo)為一個離散值仔沿，所以不可避免的坐桩，兩個相鄰矩形之間一定會存在空隙，所以這里我們要想辦法把離散值轉(zhuǎn)換為連續(xù)值封锉，這樣的話就能實(shí)現(xiàn)兩個矩形之間的“無縫銜接”绵跷。
這里我們用一個新的數(shù)據(jù)框作為geom_rect()函數(shù)的輸入，這個數(shù)據(jù)框是這樣的：
rect.data <- data.frame(
ymin = c(0, 1.5, 2.5, 3.5, 4.5, 10.5, 17.5),
ymax = c(1.5, 2.5, 3.5, 4.5, 10.5, 17.5, 23.5),
colors = letters[1:7]
)
head(rect.data)
# ymin ymax colors
#1 0.0 1.5 a
#2 1.5 2.5 b
#3 2.5 3.5 c
#4 3.5 4.5 d
#5 4.5 10.5 e
#6 10.5 17.5 f
現(xiàn)在可以添加矩形框了：
p1 +
geom_rect(data = rect.data, aes(xmin = -Inf, xmax = Inf, ymin = ymin, ymax = ymax, fill = colors), inherit.aes = FALSE, alpha = 0.2, show.legend = FALSE) +
scale_fill_manual(values = c('#686868', '#91793C', '#E5B63D', '#E2348E', '#7974A1', '#CF6611', '#279D77')) +
guides(color = guide_colorbar(order = 1), size = guide_legend(order = 2))
出圖：
這里有兩個點(diǎn)需要強(qiáng)調(diào)一下：
（1）geom_rect()一定要加inherit.aes = FALSE成福，如果不加就會報錯碾局，因?yàn)樗^承的映射在這個數(shù)據(jù)框中沒有相應(yīng)的元素；
（2）使用guides()函數(shù)改變圖例順序奴艾，否則是和原文不一樣的净当。
寫在最后
其它與原圖不符的我會選擇直接在Adobe Illustrator進(jìn)行更改，更加方便蕴潦，在代碼方面就不贅述了像啼。
