R| ggseg 绘制统计结果

最新推荐文章于 2025-05-04 19:41:26 发布

懒麻蛇

最新推荐文章于 2025-05-04 19:41:26 发布

阅读量3.3k

点赞数 12

文章标签： python java 大数据数据库 linux

本文链接：https://blog.csdn.net/lazysnake666/article/details/122404780

版权

ggseg是2018年出的R工具包，可以在R中绘制到皮层或者皮层下区域的统计结果。和brainconn一样，解决了需要从R导出数据用其他软件作图的问题。ggseg基于ggplot，因此能用ggplot的方式调整作图效果，比如facet_wrap和theme，同时还兼容了dplyr的管道命令(%>%)。

在最近1.6的版本中做出了很多主要的更新，引入了geom_brain/geom_sf，更符合ggplot的规范，让作图的逻辑更加清晰。之前使用的ggseg()的作图方式仍然会保留一段时间，但是将来会被geom的方式取代。

一直在用1.5.3不想升级的人👇

软件目前自带两个常用模板：

1. aparc (也就是Desikan-Killany atlas)

2. aseg (Automatic subcortical segmentation)

当然还可以加载一些作者制作的模板。

ggseg容易上手，需要做的就是告诉软件每个region/label的数据是什么。本文介绍基本功，基于目前最新的版本1.6.02，详情参考软件的文章和vignettes。

Mowinckel, A. M., & Vidal-Piñeiro, D. (2020). Visualization of Brain Statistics With R Packages ggseg and ggseg3d. Advances in Methods and Practices in Psychological Science, 3(4), 466-483.

安装

install.packages("remotes")
install_github("LCBC-UiO/ggseg", build_vignettes = TRUE)

⚠️github中有一个vignette有一些路径设置的问题，如果安装出错的话，build_vignettes选择FALSE。GitHub中的vignettes，但是不需要下载和knit，到软件主页查看作者织好的。

数据

比如需要绘制以下四个区域的统计值:

library(dplyr) 
someData = tibble(
  region = c("transverse temporal", "insula",
           "precentral","superior parietal"), 
  p = sample(seq(0,.5,.001), 4)
)

这里用region作为匹配，需要region和模板中的名称一致，可以看到它不是缩写，有空格，而且没有left/right的信息。另外还可以使用label作为匹配，同样需要和模板中一致，label的形式更接近于实际使用中的情况。

someData = tibble(
  label = c("lh_bankssts", "lh_insula",
             "lh_parsopercularis","lh_parsorbitalis"),
  p = sample(seq(0,.5,.001), 4)
)

下面的例子用region的数据，但实际使用中推荐label的方式。

作图

geom_brain() vs. geom_sf()

geom_brain和geom_sf用的是不同的风格，查看帮助文件会发现后者允许设置的内容更多，灵活度更高。

geom_brain

p1=ggplot(someData) +
  geom_brain(atlas = dk, 
             position = position_brain(hemi ~ side),
             aes(fill = p))p1

geom_sf

p1=someData %>% 
  brain_join(dk) %>% 
  reposition_brain(hemi ~ side) %>% 
  ggplot() + 
  geom_sf(aes(fill = p))p1

其中position_brain默认是水平，它设置是 rows ~ columns organisation。在geom_sf作图中需要使用reposition_brain的方式实现。

改theme

p1+theme_void()

使用theme_void()/theme_brain2()是快速去掉坐标轴的方法，也可以使用theme_custombrain()/theme()做更多设置。⚠️ggseg用了mono字体，如果要拼图的话需要注意一些字体的统一性。

theme_custombrain(
  plot.background = "white",
  text.colour = "darkgrey",
  text.size = 12,
  text.family = "mono"
)

改colorbar

scale_fill_系列的函数都适用。

p1+scale_fill_gradient(low = "grey", high = "brown")+
  theme_void()

p1+scale_fill_distiller(palette = "RdPu")+
  theme_void()

p1+scale_fill_viridis_c(option = "cividis", direction = -1)+
  theme_void()

基本功能介绍到此结束，一些细节的修改比如增加title，修改legend的标题和位置，设置colorbar的范围等和ggplot的设置类似。对于subcortical volume的作图方法类似，同样是添加相应region/label的数据。最后，推荐用label+geom_sf的方式进行作图。

Extra

原理

dk其实是一堆数据，如果直接画dk中data会发现

plot(dk$data)

软件作图的逻辑其实是在dk的数据中增加一个新的column，作图时指定需要fill的column是哪个一个。

brain_join可以根据label或者region添加新column到dk中。

someData %>%
  brain_join(dk)

当然也可以手动进行添加

dk %>% 
  as_tibble() %>% 
  left_join(someData)

添加p值的column后如果不指定需要画哪一个column，得到下图

someData %>% 
  brain_join(dk) %>% 
  plot()

使用facet_wrap

数据里增加分组的变量group，可以使用facet_wrap同时绘图。facet_wrap相比于先分组制图，后使用grid等工具拼接的方式，优点在于能用同一个colorbar。

someData <- tibble(
  region = rep(c("transverse temporal", "insula",
                 "precentral","superior parietal"), 2),
  p = sample(seq(0,.5,.001), 8),
  groups = c(rep("点赞组", 4), rep("不点赞组", 4))
)
someData %>%
  group_by(groups) %>%
  brain_join(dk) %>%
  reposition_brain(hemi ~ side) %>%
  ggplot(aes(fill = p)) +
  geom_sf(show.legend = FALSE) +
  facet_wrap( ~ groups) +
  theme_void()+
  theme(text = element_text(family='Kai'))

读取freesurfer文件作图

作者提供了读取freesurfer统计结果文件的方式，可以一次读取统计文件中ROI的各种指标。"aparc.stats$"是正则表达式，说明lh/rh都会读取。

dat <- read_atlas_files(subject_dir, "aparc.stats$")
dat
#> # A tibble: 68 x 11
#> subject label NumVert SurfArea GrayVol ThickAvg ThickStd MeanCurv GausCurv
#> <chr> <chr> <int> <int> <int> <dbl> <dbl> <dbl> <dbl>
#> 1 bert lh_bank… 1181 831 2297 2.77 0.428 0.116 0.024
#> 2 bert lh_caud… 843 572 1534 2.72 0.469 0.124 0.013
#> 3 bert lh_caud… 2758 1840 5772 2.80 0.526 0.114 0.021
#> 4 bert lh_cune… 2683 1654 3074 1.81 0.471 0.14 0.033
#> 5 bert lh_ento… 581 416 1840 3.33 0.667 0.116 0.032
#> 6 bert lh_fusi… 4113 2875 8519 2.72 0.599 0.131 0.025
#> 7 bert lh_infe… 4948 3466 10559 2.70 0.509 0.131 0.029
#> 8 bert lh_infe… 5056 3542 12358 2.98 0.625 0.138 0.033
#> 9 bert lh_isth… 1561 990 2350 2.09 0.745 0.113 0.023
#> 10 bert lh_late… 7961 5077 12743 2.30 0.588 0.139 0.033
#> # … with 58 more rows, and 2 more variables: FoldInd <int>, CurvInd <dbl>

画cortical thickness std

ggseg(dat, mapping = aes(fill = ThickStd))

同时画几个指标。这只是功能展示，因为几个指标不在一个scale上。

dat %>% 
  gather(stat, val, -subject, -label) %>% 
  group_by(stat) %>% 
  ggseg(mapping = aes(fill = val)) +
  facet_wrap(~stat)

实例

两个组s1和s2分别和正常组做了dk-68个区域皮层厚度的比较，结果如下：

制作effect size的统计图

#替换掉前缀后缀使得label和ggseg的一致
label=df$X
label=sub('combat_L','lh',label,ignore.case = TRUE)
label=sub('combat_R','rh',label,ignore.case = TRUE)
label=sub('_thickavg','',label,ignore.case = TRUE)

es=c(df$s1_es_vals,df$s2_es_vals)

#数据
someData = tibble(
  label = rep(label,2), 
  es = es,
  groups = c(rep("s1", 68), rep("s2", 68))
)
#作图
someData %>%
  group_by(groups) %>% 
  brain_join(dk) %>% 
  reposition_brain(hemi ~ side) %>% 
  ggplot(aes(fill = es)) + 
  geom_sf(show.legend = TRUE) + 
  facet_wrap( ~ groups) +
  scale_fill_viridis_c(option = "cividis")+
  theme_void()