11 样本信息多轨道环形图（Nature Aging 2024 Fig 1a 复现）

#!/usr/bin/env Rscript
# ============================================================================
# 00_simulate_participants.R
#
# 按原文结果图的分布特征生成模拟数据 participants.csv。
# 原文 participants.csv 为付费软件 PlotCase 内置数据（未公开），此处仅复刻
# 其边际分布，保证汇总面板的数值与原图一致：
#   - Age / BMI：五数概括（范围 + 三分位数）与原图标签一致
#     Age  Q1=48.7  中位=55.7  Q3=61.0 （原图灰底标签）
#     BMI  Q1=25.2  中位=28.2  Q3=31.1 （原图灰底标签）
#   - 类别人数精确匹配原图百分比（n = 108）：
#     Gender   Female 56 (51.9%) / Male 52 (48.1%)
#     Ethnicity Caucasian 71 (65.7%) / Asian 22 (20.4%) / Black 8 (7.4%)
#                / Hispanics 4 (3.7%) / 缺失 3 (2.8% -> Missing)
#     IRIS     IR 36 (33.3%) / IS 34 (31.5%) / Unknown 38 (35.2%)
#   - ID 风格复刻：全 Z 开头 7 位字母数字
# 输出: participants.csv（与脚本同目录）
# ============================================================================

set.seed(20261006)

n <- 108

# 分段线性分位数映射：u ∈ (0,1) 升序 → 目标五数概括 [min, Q1, 中位, Q3, max]
# 用确定性格点 (k-0.5)/n 而非 runif，避免经验分位点随机漂移导致
# quantile() 结果偏离目标（type-7 分位与格点线性插值可精确对齐）
qmap <- function(n, fivenum, digits = 1) {
  u <- (seq_len(n) - 0.5) / n
  round(approx(c(0, .25, .5, .75, 1), fivenum, xout = u, ties = "ordered")$y,
        digits)
}

# 连续变量 --------------------------------------------------------------
age <- qmap(n, c(19, 48.7, 55.7, 61, 68), 1)      # Age 升序（后续按 Age 排序）
bmi <- sample(qmap(n, c(18.8, 25.2, 28.2, 31.1, 41), 1))  # BMI 打乱：与 Age 独立
bmi[sample(n, 2)] <- NA                            # 缺失 BMI：原图保留空白轨道

# 类别变量（人数精确匹配，整体打乱）--------------------------------------
gender    <- sample(rep(c("Female", "Male"), c(56, 52)))
ethnicity <- sample(rep(c("Caucasian", "Asian", "Black", "Hispanics", NA),
                        c(71, 22, 8, 4, 3)))       # 3 例缺失 -> Missing 2.8%
iris      <- sample(rep(c("IR", "IS", "Unknown"), c(36, 34, 38)))

# ID：Z 开头 7 位（有放回抽样，允许字符重复，贴近原图 ID 风格）-----------
make_id <- function() {
  paste0("Z", paste(sample(c(LETTERS, 0:9), 6, replace = TRUE), collapse = ""))
}
ids <- vapply(seq_len(n), function(i) make_id(), character(1))
while (anyDuplicated(ids)) {                       # 极小概率重复时重抽
  dup <- duplicated(ids)
  ids[dup] <- vapply(which(dup), function(i) make_id(), character(1))
}

participants <- data.frame(
  ID       = ids,
  Age      = age,
  BMI      = bmi,
  Gender   = gender,
  Ethnicity = ethnicity,
  IRIS     = iris,
  stringsAsFactors = FALSE
)

out <- "participants.csv"
script_arg <- grep("^--file=", commandArgs(), value = TRUE)
if (length(script_arg) > 0L) {
  out <- file.path(dirname(normalizePath(sub("^--file=", "", script_arg[1]),
                                          mustWork = FALSE)),
                   "participants.csv")
}

write.csv(participants, out, row.names = FALSE, na = "NA")

# ---- 分布自检（应与原图标签一致）---------------------------------------
cat("Age  Q1/中位/Q3:", paste(round(quantile(age, c(.25, .5, .75)), 1),
                              collapse = " / "), "\n")
cat("BMI  Q1/中位/Q3:", paste(round(quantile(bmi, c(.25, .5, .75), na.rm = TRUE), 1),
                              collapse = " / "), "\n")
cat("Gender:  ", paste(names(table(gender)), table(gender),
                       sep = "=", collapse = "  "), "\n")
cat("IRIS:    ", paste(names(table(iris)), table(iris),
                       sep = "=", collapse = "  "), "\n")
cat("written:", out, "\n")

#!/usr/bin/env Rscript
# ============================================================================
# 11_sample_info_multitrack_ring.R
#
# 样本信息多轨道环形图：108 位参与者按年龄排序后分布在同一圆周（270° 开口
# 朝左上），由内向外依次为胰岛素反应 / 族裔 / 性别三条类别条轨道，以及
# 年龄 / BMI 两条棒棒糖（lollipop）点轨道；环内左上开口嵌入样本汇总面板
# （Age/BMI 分布箱线图 + 三类别占比堆叠条），顶部为组合图例。
#
# 来源: 微信公众号「R语言数据分析指南」
#       《Nature Aging 图表复现：样本信息多轨道环形图》
#       https://mp.weixin.qq.com/s/_egWnpkkc56heFB4YGHCfA
# 论文: Shen et al. Nonlinear dynamics of multi-omics profiles during human
#       aging. Nature Aging (2024). Figure 1a
#
# 转录说明（原文代码以三段截图给出：右缘截断、p_summary/p_legend 构建缺失）:
#   [转录] = 截图原文逐字
#   [补全] = 截图右缘截断，按 R 语法与上下文补全（推断）
#   [重写] = 原文截图缺失，按结果图重写（推断）
# 输入: participants.csv（模拟数据，见 00_simulate_participants.R / README）
# 输出: output/sample_info_multitrack_ring.pdf
# 运行: Rscript 11_sample_info_multitrack_ring.R
# ============================================================================

suppressPackageStartupMessages({
  if (!requireNamespace("ggplot2", quietly = TRUE)) stop("缺少 ggplot2")
  if (!requireNamespace("dplyr",  quietly = TRUE)) stop("缺少 dplyr")
  if (!requireNamespace("tidyr",  quietly = TRUE)) stop("缺少 tidyr")
  if (!requireNamespace("patchwork", quietly = TRUE)) stop("缺少 patchwork")
  library(ggplot2); library(dplyr); library(tidyr); library(patchwork)
})

# ---- 路径：脚本所在目录为模块根 ------------------------------------------
args <- commandArgs(trailingOnly = TRUE)
script_arg <- grep("^--file=", commandArgs(), value = TRUE)
module_dir <- if (length(script_arg) > 0L) {
  dirname(normalizePath(sub("^--file=", "", script_arg[1]), mustWork = FALSE))
} else {
  getwd()
}

# ============================================================================
# 1. 数据读入与排序  [转录 + 补全]
# ============================================================================
participants <- read.csv(file.path(module_dir, "participants.csv"),
                         na.strings = c("", "NA")) %>%   # [补全] 原文无路径前缀
  arrange(Age) %>%
  mutate(participant_order = row_number(),
         Ethnicity = replace_na(Ethnicity, "Missing"))   # [补全] 行尾截断

# ============================================================================
# 2. 配色  [转录 + 补全]
# ============================================================================
category_colors <- c(
  IS = '#F01800', IR = '#3F9CAF', Unknown = '#BEBEBE',
  Female = '#F12D0C', Male = '#DFBF70', Asian = '#08739A',
  Black = '#A5CEDE',                 # [补全] 色值截断，按结果图取色
  Caucasian = '#EAD0B2', Hispanics = '#D29C4B', Missing = '#BEBEBE')
measurement_colors <- c(Age = '#702380', BMI = '#FF4B16')

# ============================================================================
# 3. 轨道布局（内 -> 外 5 条轨道的半径）  [转录 + 补全]
# ============================================================================
track_layout <- data.frame(
  # [补全] 截图仅可见前 3 项与全部半径；按原图核对：点轨 BMI 在内、Age 在外
  variable     = c('IRIS', 'Ethnicity', 'Gender', 'BMI', 'Age'),
  inner_radius = c(1.06, 1.66, 2.26, 2.85, 3.94),
  outer_radius = c(1.56, 2.16, 2.76, 3.44, 4.53))  # [补全] 后 3 项截断（轨道等厚 0.5）

# ============================================================================
# 4. 数据整形（长表 + 半径映射 + 刻度 + 标签角度）  [转录 + 补全]
# ============================================================================
participant_categories <- participants %>%
  pivot_longer(c(IRIS, Ethnicity, Gender), names_to = 'variable',
               values_to = 'category') %>%               # [补全] values_to 截断
  left_join(track_layout, by = 'variable')

participant_measurements <- participants %>%
  pivot_longer(c(Age, BMI), names_to = 'variable', values_to = 'value') %>%
  left_join(track_layout, by = 'variable') %>%
  group_by(variable) %>%
  mutate(scale_min = .8 * min(value, na.rm = TRUE),
         scale_max = 1.1 * max(value, na.rm = TRUE),      # [补全] max 行尾截断
         radius = inner_radius + (value - scale_min) /
           (scale_max - scale_min) *
             (outer_radius - inner_radius)) %>%           # [补全] 行尾截断
  ungroup()

# 径向刻度：每变量取 [缩放下限, 值域中点, 缩放上限] 三点，映射到对应半径
radial_ticks <- participant_measurements %>%
  group_by(variable) %>%
  reframe(value = c(first(scale_min), mean(range(value, na.rm = TRUE)),
                    first(scale_max)),                    # [补全] 第 3 项截断
          radius = first(inner_radius) + (value - first(scale_min)) /
            (first(scale_max) - first(scale_min)) *       # [补全] 公式截断
            (first(outer_radius) - first(inner_radius)))

# 参与者 ID 的放射状标签角度：270° 均分，左上开口区文字翻转保持可读
participant_labels <- participants %>%
  mutate(angle = 90 - (participant_order - .5) / n() * 270,
         hjust = if_else(angle < -90, 1, 0),              # [补全] 截断
         angle = if_else(angle < -90, angle + 180, angle))

# ============================================================================
# 5. 环形主体  [转录 + 补全]
# ============================================================================
p_circle <- ggplot() +
  # 类别条轨道（IRIS / Ethnicity / Gender）
  geom_rect(data = participant_categories,
            aes(xmin = participant_order - .5, xmax = participant_order + .5,  # [补全]
                ymin = inner_radius, ymax = outer_radius, fill = category),
            color = '#333333', linewidth = .3) +          # [补全] 参数截断
  # 点轨道白底格（Age / BMI），缺失值留空
  geom_rect(data = participant_measurements,
            aes(xmin = participant_order - .5, xmax = participant_order + .5,  # [补全]
                ymin = inner_radius, ymax = outer_radius),
            fill = 'white', color = '#333333', linewidth = .3) +
  # 棒棒糖杆
  geom_segment(data = participant_measurements,
               aes(participant_order, inner_radius,       # [补全] y 截断
                   xend = participant_order, yend = radius, color = variable),
               linewidth = .45, na.rm = TRUE) +
  # 棒棒糖点
  geom_point(data = participant_measurements,
             aes(participant_order, radius, color = variable),  # [补全] 截断
             size = 1.6, na.rm = TRUE) +                  # [补全] size 截断
  # 参与者 ID（外圈放射状）
  geom_text(data = participant_labels,
            aes(participant_order, 5.02, label = ID,      # [补全] label 截断
                angle = angle, hjust = hjust), size = 3.1) +
  # 径向刻度数字（开口顶端，水平文字）—— 统一 2.6mm
  geom_text(data = radial_ticks,
            aes(x = .5, y = radius, label = round(value)),  # [补全] label 截断
            size = 2.6, hjust = 1.15) +                   # [补全] 参数截断
  scale_fill_manual(values = category_colors) +
  scale_color_manual(values = measurement_colors) +
  scale_x_continuous(limits = c(.5, nrow(participants) * 4 / 3 + .5),
                     expand = c(0, 0)) +                  # [补全] expand 截断
  scale_y_continuous(limits = c(0, 5.85), expand = c(0, 0)) +
  coord_polar(start = 0, clip = 'off') +
  theme_void() +
  theme(legend.position = 'none', plot.margin = margin(0, 0, 0, 0))

# ============================================================================
# 6. 样本汇总面板  [补全 / 重写]
#    原文三段截图不含 p_summary 构建代码（只到 summary_theme 定义），
#    以下按结果图重写：Age/BMI 分布面板 + 三类别占比面板。
# ============================================================================

# 四分位数（箱线图上的灰底数值标签）  [转录]
measurement_quartiles <- participant_measurements %>%
  group_by(variable) %>%
  reframe(value = quantile(value, c(.25, .5, .75), na.rm = TRUE))

# 汇总面板主题  [转录 + 补全]
summary_theme <- theme_bw(base_size = 10) +
  theme(panel.grid = element_blank(), axis.title = element_blank(),
        axis.text = element_blank(),                      # [补全] 截断
        axis.ticks = element_blank(), strip.background = element_blank(),
        # 统一字号体系（与图内 geom 文字 2.6mm 对齐）：strip = 7.4pt ≈ 2.6mm
        strip.text = element_text(size = 7.4),            # [重写] 面板标题样式
        plot.margin = margin(4, 2, 2, 2), legend.position = 'none')

# ---- 6a. Age / BMI 分布面板：竖箱线 + 抖动点 + 四分位灰底标签  [重写] ----
labeller_measures <- labeller(variable = c(Age = 'Age\n(years)', BMI = 'BMI'))

p_dist <- ggplot(participant_measurements, aes(variable, value)) +
  geom_jitter(aes(color = variable), width = .35, height = 0,
              size = .9, na.rm = TRUE) +
  # 箱线画在点之上、填充透明：框线清晰且 box 内可见点（对齐原图）
  geom_boxplot(color = 'black', fill = NA, linewidth = .35,
               width = .7, outlier.shape = NA, na.rm = TRUE) +
  geom_label(data = measurement_quartiles,
             aes(variable, value, label = sprintf('%.1f', value)),
             inherit.aes = FALSE, fill = 'grey62', color = 'white',
             fontface = 'bold', size = 2.6, linewidth = 0) +
  scale_color_manual(values = measurement_colors) +
  # free：x 离散 scale 不跨 panel 共享，单 level 才落在 panel 中心
  # （默认 fixed 时 Age=位置1/BMI=位置2 落在共享范围 [0.4,2.6] 内 → 各偏半格）
  facet_wrap(~ variable, scales = 'free', labeller = labeller_measures) +
  summary_theme

# ---- 6b. 类别占比面板：单列堆叠竖条 + 百分比标签  [重写] ------------------
category_summary <- participant_categories %>%
  mutate(variable = factor(variable, levels = c('Gender', 'Ethnicity', 'IRIS'))) %>%
  count(variable, category, name = 'n') %>%
  group_by(variable) %>%
  # 标签区间与 geom_col 的堆叠同序：堆叠按 category 字母序反向（底->顶），
  # cumsum 前先按字母序反排，保证百分比标签落在自己的色块上
  arrange(desc(category), .by_group = TRUE) %>%
  mutate(pct = 100 * n / sum(n),
         ymax = cumsum(pct), ymin = ymax - pct) %>%
  ungroup() %>%
  mutate(
    # 统一字号分级（2.6mm 主档 / 窄段降档防跨出色块：7.4% 段仅 ~26px 高，用 1.6mm）
    label_size = ifelse(pct < 5, 1.1, ifelse(pct < 15, 1.6, 2.6)),
    # <5% 的极窄段（3.7/2.8 相邻）水平错开，避免竖排标签互叠
    nudge = ifelse(pct < 5, ifelse(pct < 3.2, 0.18, -0.18), 0))

labeller_cats <- labeller(variable = c(
  Gender = 'Gender\n(%)', Ethnicity = 'Ethnicity\n(%)',
  IRIS = 'Insulin\nresponse\n(%)'))

# 按色块亮度自动选择标签字色（浅块深字、深块白字）
cat_luma <- vapply(category_summary$category,
                   function(cl) sum(grDevices::col2rgb(category_colors[cl])) / 3,
                   numeric(1))
category_summary$label_color <- ifelse(cat_luma > 170, '#333333', 'white')

p_cat <- ggplot(category_summary, aes(variable)) +
  geom_col(aes(y = pct, fill = category), width = .7) +
  geom_text(aes(y = (ymin + ymax) / 2, label = sprintf('%.1f', pct),
                color = label_color),
            # size 作图层参数（向量逐行生效），避免被 scale_size 二次映射；
            # nudge 水平错开极窄段标签（见 category_summary$nudge）
            size = category_summary$label_size,
            position = position_nudge(x = category_summary$nudge),
            angle = 90, fontface = 'bold', show.legend = FALSE) +
  scale_fill_manual(values = category_colors) +
  scale_color_identity() +
  # free_x：每 panel 只保留自己的 level → 条形居中（同 p_dist 的共享 scale 问题）
  facet_wrap(~ variable, scales = 'free_x', labeller = labeller_cats) +
  summary_theme

# 左侧面板 = 分布面板 | 占比面板（类别格更宽以容纳 7.4pt strip 标题全称）
p_summary <- (p_dist | p_cat) + plot_layout(widths = c(1, 1.15))

# ============================================================================
# 7. 顶部组合图例  [重写]
#    原文截图缺失 p_legend 构建；图例 18 项 3 行 6 列照抄结果图
#    （含 10 个组学数据类型 —— 论文原图图例项，本复现 5 轨道未使用，
#    结果图亦如此保留）。
# ============================================================================
legend_items <- data.frame(
  label = c('Insulin sensitivity', 'Insulin resistance', 'Asian', 'Black',
            'Caucasian', 'Hispanics',
            'Female', 'Male', 'Transcriptomics', 'Proteomics',
            'Metabolomics', 'Cytokine',
            'Clinical lab test', 'Lipidomics', 'Gut', 'Skin', 'Oral', 'Nasal'),
  color = c('#F01800', '#3F9CAF', '#08739A', '#A5CEDE', '#EAD0B2', '#D29C4B',
            '#F12D0C', '#DFBF70', '#7BC8A4', '#F5E642', '#9B93C8', '#F87171',
            '#5B9BD5', '#F5924B', '#8CC63F', '#F8B3D9', '#9B59B6', '#B7E4A0'),
  stringsAsFactors = FALSE)   # 后 10 项组学配色为结果图近似取色 [重写]

legend_items$label <- factor(legend_items$label, levels = legend_items$label)

p_legend <- ggplot(legend_items, aes(fill = label)) +
  geom_col(aes(1, 1), width = 1) +
  scale_fill_manual(values = setNames(legend_items$color, levels(legend_items$label)),
                    name = NULL) +
  guides(fill = guide_legend(nrow = 3, byrow = TRUE)) +
  theme_void() +
  theme(legend.position = 'bottom', legend.margin = margin(0, 0, 0, 0),
        legend.key.size = unit(7, 'pt'),
        # 统一字号体系：图例文字 7.4pt ≈ 2.6mm（与 strip / geom 主文字同档）
        legend.text = element_text(size = 7.4),
        legend.spacing.y = unit(2, 'pt'))

# ============================================================================
# 8. 组装输出  [转录 + 补全]（inset 位置参数按本版汇总面板比例微调）
# ============================================================================
legend_table <- ggplotGrob(p_legend)
# ggplot2 >= 3.5 时 5 个 guide-box-* 占位并存（右侧为空 zeroGrob），
# 需取非空的那个（原文截图条件 == 'guide-box' 截断，此处补全其等价逻辑）
guide_idx <- which(grepl('^guide-box', legend_table$layout$name))
guide_idx <- guide_idx[!vapply(legend_table$grobs[guide_idx], function(g) {
  inherits(g, 'zeroGrob')
}, logical(1))]
p_legend <- wrap_elements(full = legend_table$grobs[[guide_idx[1]]])

p_body <- p_circle +
  # 面板精确卡入环形左上 1/4 象限（以实测像素为准：框右贴 θ0 射线、
  # 框底贴 270° 射线=环心线、框顶贴环顶 r=4.53 弧；inset 与框之间差
  # strip/plot.margin 固定偏移，已折算进参数）
  # 面板与图形外缘的 ID 文字排衔接（不留左上空段）：
  #   左缘 = 最左 ID 文字列的右端（文字排尽头线）
  #   顶缘 = 最上 ID 文字列的底端（同理，strip 紧贴其下）
  #   底缘 = 270° 射线（环心水平线）  右缘 = θ0 射线（环心垂直线）
  inset_element(p_summary, left = .124, bottom = .505, right = .502, top = .835)

p <- (p_legend / p_body) +
  plot_layout(heights = c(.13, 1))

out_dir <- file.path(module_dir, 'output')
if (!dir.exists(out_dir)) dir.create(out_dir, recursive = TRUE)
out_pdf <- file.path(out_dir, 'sample_info_multitrack_ring.pdf')
ggsave(out_pdf, p, width = 10, height = 10.3, device = 'pdf')
cat('written:', out_pdf, '\n')