跟着PNAS学作图 | 提供全文数据和代码

2023-03-29 本文已影响0人小杜的生信筆記

论文

题目：Death rates at specific life stages mold the sex gap in life expectancy

网址： https://www.pnas.org/doi/full/10.1073/pnas.2010588118

代码网址

https://github.com/CPop-SDU/sex-gap-e0-pnas

该文章发表于2021年，论文中图形对我们一部分同学仍具参考价值。作者提供的全套的代码和数据，可以直接使用。此外，作者的数据和代码写的非常的规整。但是，需要看懂和运行代码，还是需要有一定的基础。

论文主图

论文主图仅有两张，如下图所示。

代码

Figure 1


# function to localize paths

devtools::source_gist("32e9aa2a971c6d2682ea8d6af5eb5cde")

# prepare session
source(lp("0-prepare-session.R"))


# theme -------------------------------------------------------------------
load("../dat/palettes.rda" %>% lp)

theme_custom <- theme_minimal(base_family = font_rc) +
    theme(
        legend.position = "bottom",
        strip.background = element_blank(),
        strip.text = element_blank(),
        panel.grid.minor =  element_blank(),
        panel.grid.major =  element_line(size = .25),
        panel.ontop = T
    )

作者将相关的代码编写在其他的R脚本中，使用时直接进行调用。
[图片上传失败...(image-122af8-1680158585814)]

# Fig 1 -- RELATIVE ----------------------------------
load("../dat/a6gap33cntrs.rda" %>% lp)

# relative
df6 %>% 
    filter(country %>% is_in(c("SWE", "USA", "JPN", "RUS"))) %>%
    mutate(
        name = name %>% 
            fct_recode(USA = "United States") %>% 
            fct_rev()
    ) %>%
    ggplot() +
    geom_col(
        aes(year, ctb_rel %>% multiply_by(100), fill = age_group),
        position = position_stack(reverse = TRUE),
        color = NA,
        width = 1
    ) +
    facet_grid(name ~ ., scales = "free_y", space = "free") +
    coord_cartesian(ylim = c(-10, 120), expand = FALSE)+
    scale_x_continuous(breaks = seq(1800, 2000, 50))+
    scale_y_continuous(breaks = seq(0, 100, 25), position = "right")+
    scale_fill_manual(
        values = pal_six, 
        guide  = guide_legend(ncol = 1, reverse = TRUE)
    ) +
    theme_minimal(base_family = font_rc, base_size = 20) +
    theme(
        legend.position = c(.6, .5),
        strip.background = element_blank(),
        strip.text = element_blank(),
        panel.grid.minor =  element_blank(),
        panel.grid.major =  element_line(size = .1),
        panel.spacing = unit(0, "lines"),
        panel.ontop = T
    )+
    labs(x = NULL,
         y = "Contribution, %",
         fill = "Age group")+
    # label countries
    geom_text(data = . %>% select(name, row, column) %>%  distinct(),
              aes(label = name, color = name), 
              x = 2015, y = 120, 
              hjust = 1, vjust = 1, size = 9, fontface = 2,
              family = font_rc)+
    scale_color_manual(values = pal_four %>% rev, 
                       guide = FALSE)

one_outer <- last_plot()
one_outer

# plot ratio
load("../dat/df4qx.rda" %>% lp)

df4qx %>%
    pivot_wider(names_from = sex, values_from = qx) %>% 
    ggplot(aes(age, y = m/f, color = country))+
    geom_hline(yintercept = 1, color = "gray25",  size = .5)+
    geom_smooth(se = F, size = 1, color = "#ffffff", span = .25)+
    geom_smooth(se = F, size = .5, span = .25)+
    scale_x_continuous(breaks = c(0, 15, 40, 60, 80))+
    scale_y_continuous(
        trans = "log", 
        breaks = c(.5, 1, 2, 3), 
        labels = c("", 1, 2, 3),
        limits = c(.75, 3.5)
    )+
    scale_color_manual(NULL, values = pal_four)+
    theme_minimal(base_family = font_rc, base_size = 16)+
    theme(
        legend.position = "none",
        panel.grid.minor = element_blank()
    )+
    labs(
        y = "Sex ratio, log scale",
        x = "Age"
    )+
    annotate(
        "text", x = 50, y = .9, 
        label = "Most recent year",
        size = 8.5, color = "grey50", alpha = .5,
        vjust = 1, family = font_rc, fontface = 2
    )

one_a <- last_plot()
one_a
# Death risk Ratio, Sweden, years 1750, 1800, 1850, 1900, 1960, 2019
# plot qx
load("../dat/qxdiff.rda" %>% lp)

qxdiff %>% 
    filter(country == "SWE", 
           year %>% is_in(c(1800, 1900, 1960, 2019 ))) %>% 
    ggplot(aes(age, y = ratio, color = year %>% factor))+
    geom_hline(yintercept = 1, color = "gray25",  size = .5)+
    geom_smooth(se = F, size = .75, span = .4)+
    scale_x_continuous(breaks = c(0, 15, 40, 60, 80))+
    scale_y_continuous(
        trans = "log", 
        breaks = c(.5, 1, 2, 3), 
        labels = c("", 1, 2, 3),
        limits = c(.75, 3.5)
    )+
    scale_color_viridis_d(end = .97)+
    theme_minimal(base_family = font_rc, base_size = 16)+
    theme(
        legend.position = c(.85, .75),
        legend.spacing.x = unit(.1, "line"),
        legend.key.height = unit(1, "line"),
        panel.grid.minor = element_blank()
    )+
    labs(
        color = "Year",
        y = "Sex ratio, log scale",
        x = "Age"
    )+
    annotate(
        "text", x = 50, y = .9, 
        label = "Sweden",
        size = 8.5, color = "#009C9C", 
        vjust = 1, family = font_rc, fontface = 2
    )

one_b <- last_plot()
one_b

# plot difference
df4qx %>%
    pivot_wider(names_from = sex, values_from = qx) %>% 
    ggplot(aes(x = age, y = m - f, color = country, group = country)) +
    geom_path(size = .5)+
    scale_color_manual(NULL, values = pal_four)+
    scale_x_continuous(breaks = c(0, 15, 40, 60, 80))+
    scale_y_continuous(
        trans = "log",
        breaks = c(.0001, .001, .01, .05),
        labels = c(.0001, .001, .01, .05) %>% paste %>% str_replace("0.", "."),
        limits = c(9e-6, .1)
    )+
    theme_minimal(base_family = font_rc, base_size = 16)+
    theme(legend.position = c(.77, .25),
          legend.spacing.x = unit(.1, "line"),
          legend.key.height = unit(1, "line"),
          legend.text = element_text(size = 16),
          panel.grid.minor = element_blank())+
    labs(
        y = "Sex gap, log scale",
        x = "Age"
    )

one_c <- last_plot()
one_c

# arrange and save
blank <- ggplot(tibble(x = 1, y = 1), aes(x, y))+
    geom_rect(xmin = -Inf, xmax = Inf,
              ymin = -Inf, ymax = Inf,
              fill = "#ffffff",
              color = NA)+
    theme_void()

library(cowplot)
one <- ggdraw() +
    draw_plot(one_outer) +
    # white space for plots
    draw_plot(blank, x = 0, y = .75, width = 0.7, height = 0.25)+
    draw_plot(blank, x = 0, y = .55, width = 0.33, height = 0.42)+
    draw_plot(blank, x = 0, y = .33, width = 0.33, height = 0.67)+
    # inset plots
    draw_plot(one_a, x = 0, y = .66, width = .33, height = .33)+
    draw_plot(one_c, x = .34, y = .66, width = .33, height = .33)+
    draw_plot(one_b, x = 0, y = 0.35, width = .33, height = .33)+
    # annotate plot letters
    draw_text(
        LETTERS[c(1,3,2,4)],  
        x = c(.01, .35, .01, .01),
        y = c(.99, .99, .66, .3), 
        hjust = 0,  vjust = 1, size = 20, 
        family = font_rc, fontface = 2
    )

ggsave(
    filename = "out/main-one.png" %>% lp, 
    plot = one, width = 10, height = 10, 
    type = "cairo-png"
)