MicrobiomeStatPlot | 线图教程Line plot tutorial

学术   2024-11-09 00:18   美国  

线图简介


线图最常用于绘制连续的数据。线图中点与点连成线,暗示了数据间的联系以及数据的趋势。线图可以用来反映某一变量随着另一变量变化的趋势或者说速度,一般用线段的升降来表示数值的大小变化。线图有多种多样的展现形式。

标签:#微生物组数据分析  #MicrobiomeStatPlot  #线图  #R语言可视化 #Line plot

作者:First draft(初稿):Defeng Bai(白德凤);Proofreading(校对):Ma Chuang(马闯) and Jiani Xun(荀佳妮);Text tutorial(文字教程):Defeng Bai(白德凤)

源代码及测试数据链接:

https://github.com/YongxinLiu/MicrobiomeStatPlot/项目中目录 3.Visualization_and_interpretation/LinePlot

或公众号后台回复“MicrobiomeStatPlot”领取


线图应用案例1

这是Yonatan H. Grad课题组2022年发表于Plos Biology上的文章,第一作者为Daphne S. Sun,题目为:Analysis of multiple bacterial species and antibiotic classes reveals large variation in the association between seasonal antibiotic use and resistance https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3001579

图 1 |  各类抗生素使用的季节性模式。

(A) 2011 年至 2015 年马萨诸塞州波士顿市按日历月划分的每万人日均抗生素使用量。线条表示LOESS平滑曲线,阴影区域表示95%CI。(B)每月处方率的正弦模型拟合。点表示历月每万人日均抗生素处方量的月平均季节偏差,误差条表示平均值的标准误差。线表示正弦模型振幅和相位的点估计值。阴影区域表示振幅的95%CI。星号表示季节性振幅具有统计学意义(FDR < 0.05)。FDR,错误发现率。

结果

本研究中包括的 5 类抗生素在统计上都显示出显著的季节性使用模式(图 1A)。青霉素类和大环内酯类是最常用的处方药,全年平均每天每万人分别有 4.8 和 4.1 个处方。喹诺酮类、四环素类和硝基呋喃类药物的全年平均处方量分别为每万人 1.8、1.0 和 0.54 次。

青霉素类药物的处方率在不同季节的变化幅度最大,其季节性部分的振幅为每 10,000 人中每天增加 1.1 个处方(峰值至平均值)(95% CI,0.96 至 1.3)。其次是大环内酯类(振幅为 0.74;95% CI 为 0.59 至 0.89)、喹诺酮类(振幅为 0.081;95% CI 为 0.044 至 0.12)、硝基呋喃类(振幅为 0.044;95% CI 为 0.024 至 0.064)和四环素类(振幅为 0.031;95% CI 为 0.011 至 0.052)(图 1B)。

不同抗生素类别的处方高峰时间各不相同(图 1B)。大环内酯类和青霉素类药物的使用高峰期分别出现在冬季的 1 月下旬(阶段,1.7 个月;95% CI,1.3 至 2.1;注意阶段指数为 1.0,代表 1 月 1 日)和 2 月上旬(阶段,2.2 个月;95% CI,2.0 至 2.5)。四环素和硝基呋喃类药物的使用高峰期在夏季,分别在 6 月中旬(阶段,6.5 个月;95% CI,4.6 至 8.5)和 8 月下旬(阶段,8.8 个月;95% CI,8.2 至 9.4)。最后,喹诺酮类药物的使用在每年的 1 月初和 7 月初达到两次高峰(分别为 1.0 个月(95% CI,0.6 至 1.5)和 7.0 个月(95% CI,6.6 至 7.5))。

线图应用案例2

这是Yong Liu 课题组2023年发表于NC上的文章,第一作者为Pin Su,题目为:Microbiome homeostasis on rice leaves isregulated by a precursor molecule of ligninbiosynthesis https://doi.org/10.1038/s41467-023-44335-3

图 11 |  WT、KO 和 OE 株系中 4-HCA 的浓度。

(a) 比较 LC-MS 总离子色谱图(TICs (a) WT(野生型,ZH11 水稻品种)、KOline1、KOline2、KOline3、OEline1、OEline2 和 OEline3 植物的叶提取物与 4-HCA 标准品的 LC-MS 总离子色谱图 (TIC)。

结果

不出所料,与 WT 植物相比,三个 OsPAL02-OE 株系的 4-HCA 浓度明显更高,而三个 OsPAL02-KO 株系的 4-HCA 浓度则更低(补充图 11、补充表 6 和补充数据 5)。

线图R语言实战

源代码及测试数据链接:

https://github.com/YongxinLiu/MicrobiomeStatPlot/

或公众号后台回复“MicrobiomeStatPlot”领取

软件包安装

 基于CRAN安装R包,检测没有则安装p_list = c("aplot", "tidyverse", "readxl", "cowplot", "pheatmap", "patchwork",            "ggalt", "scatterplot3d", "reshape2", "ggplot2", "ggthemes", "gridExtra")for(p in p_list){if (!requireNamespace(p)){install.packages(p)}    library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)}
# 加载R包 Load the packagesuppressWarnings(suppressMessages(library(ggplot2)))suppressWarnings(suppressMessages(library(reshape2)))suppressWarnings(suppressMessages(library(readxl)))suppressWarnings(suppressMessages(library(scatterplot3d)))suppressWarnings(suppressMessages(library(ggalt)))suppressWarnings(suppressMessages(library(patchwork)))suppressWarnings(suppressMessages(library(pheatmap)))suppressWarnings(suppressMessages(library(cowplot)))suppressWarnings(suppressMessages(library(tidyverse)))suppressWarnings(suppressMessages(library(aplot)))suppressWarnings(suppressMessages(library(ggthemes)))suppressWarnings(suppressMessages(library(gridExtra)))

1.双坐标折线图Dual-coordinate line chart

参考:

https://mp.weixin.qq.com/s/Rn6Q_xRGcXaNJtc2tn5KWw

参考文献:Haplotype-resolved genome assembly provides insights into evolutionary history of the tea plant Camellia sinensis.https://doi.org/10.1038/s41588-021-00895-y

# 模拟一些数据# Simulate some dataset.seed(123)data <- data.frame(  Position = seq(21.24, 21.32, by = 0.003),  Landrace_CSS = rnorm(27, mean = 1, sd = 0.5),  Elite_CSS = rnorm(27, mean = 0.5, sd = 0.5),  Fst = rnorm(27, mean = 0.04, sd = 0.03))
# 定义阴影区域# Define the region for the shaded areashade_start <- 21.28shade_end <- 21.30
# 田岛的 D 临界值# Tajima's D thresholdtajimas_d_threshold <- 0
# Fst 阈值# Fst thresholdfst_threshold <- 0.05fst_threshold_scaled <- fst_threshold * 15 - 1.5  # Rescale Fst as done for the line
# 绘图# Create the plotp1 <- ggplot(data, aes(x = Position)) +  geom_rect(aes(xmin = shade_start, xmax = shade_end, ymin = -Inf, ymax = Inf),            fill = "#FAE6BE", alpha = 0.3, inherit.aes = FALSE) +  geom_line(aes(y = Landrace_CSS, color = "Landrace CSS"), size = 1.2) +  geom_line(aes(y = Elite_CSS, color = "Elite CSS"), size = 1.2) +  geom_line(aes(y = Fst * 15 - 1.5, color = "Fst"), size = 1, linetype = "dashed") +  geom_hline(yintercept = tajimas_d_threshold, linetype = "dotted", color = "black", size = 0.8) +  geom_hline(yintercept = fst_threshold_scaled, linetype = "dashed", color = "red", size = 0.8) +  scale_color_manual(values = c("Landrace CSS" = "#5ebcc2", "Elite CSS" = "#cd78ac", "Fst" = "#869c51")) +  theme_bw() +  theme(    legend.title = element_blank(),    legend.position = "top",    legend.text = element_text(size = 12),    panel.grid.major = element_blank(),    panel.grid.minor = element_blank(),    axis.text = element_text(size = 12),    axis.title = element_text(size = 14, face = "bold"),    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),    plot.subtitle = element_text(size = 14, hjust = 0.5)  ) +  labs(y = "Tajima's D", x = "Genomic Position") +  ggtitle("BAS1")
# 为 Fst 添加二级 Y 轴# Add a secondary y-axis for Fstp11 <- p1 + scale_y_continuous(sec.axis = sec_axis(~ (. - 1.5) / 15, name = "Fst"))# p11
# 保存图形为PDF# Save as plotggsave("results/Plot_01.pdf", plot = p11, width = 10, height = 7)

2. 3D多分组曲线图3D multi-group curve chart 

参考:https://mp.weixin.qq.com/s/dVEhs700e3k8va9jCeqjBw

参考文献:Su, P., Kang, H., Peng, Q. et al. Microbiome homeostasis on rice leaves is regulated by a precursor molecule of lignin biosynthesis. Nat Commun 15, 23 (2024).

https://doi.org/10.1038/s41467-023-44335-3

# 加载数据# load datadf <- read.table("data/data.txt",header = T, check.names = F)
# 数据格式转换# Data format conversiondf1 <- reshape2::melt(df, id.vars = c("Wavelength"),            measure.vars = c('Amplitude_1','Amplitude_2','Amplitude_3',                     'Amplitude_4','Amplitude_5','Amplitude_6','Amplitude_7'))
# 添加轨道列# Add track columndf1$y <- rep(c(0.5,1.5,2.5,3.5,4.5,5.5,6.5),each = 512)
# 定义颜色# Defining the colourcolors <- c("#00b1c3", "#e274b6", "#5ba880", "#ab89dd", "#eb9f76", "#c28f5e", "#524438")
# 创建 PDF 文件# Creating PDF filespdf("results/Plot_02.pdf", width = 10, height = 6)
# 绘制第 1 组数据 Amplitude_1# Plotting data set 1 Amplitude_1p2 <- scatterplot3d(x = df1$Wavelength[1:512],                   y = df1$y[1:512],                   z = df1$value[1:512],                   type = 'l', lwd = 2.5,                   color = colors[1],                   scale.y = 0.9,                   y.ticklabs = c("", 'Amplitude_1', 'Amplitude_2', 'Amplitude_3',                                  'Amplitude_4', 'Amplitude_5', 'Amplitude_6', 'Amplitude_7'),                   xlim = c(min(df1$Wavelength), max(df1$Wavelength)),                   ylim = c(0, 7), zlim = c(0, 2500),                   y.axis.offset = 0.5,                   box = FALSE, grid = TRUE, angle = 45,                   xlab = paste0(colnames(df1)[1], " (nm)"),                   ylab = '', zlab = 'Amplitude (mV)'                   )
# 绘制其他组数据# Plotting other sets of datafor (i in 2:7) {  p2$points3d(df1$Wavelength[((i-1)*512+1):(i*512)],             df1$y[((i-1)*512+1):(i*512)],             df1$value[((i-1)*512+1):(i*512)],             type = 'l', col = colors[i], lwd = 2.5)}
# 添加图例 - 更高对比度且颜色相互区分清晰# Add legend - higher contrast and colours are clearly distinguishable from each otherlegend('topright', legend = c('Amplitude_1', 'Amplitude_2', 'Amplitude_3',                              'Amplitude_4', 'Amplitude_5', 'Amplitude_6', 'Amplitude_7'),       col = colors, lty = 1, bty = 'n', lwd = 2.5,       cex = 1.2, text.col = "black", inset = c(0.01, 0.01))
# 关闭 PDF 文件设备# Close PDF file devicedev.off()#> png #>   2

3.平滑曲线折线图Smooth curve line chart

参考:https://mp.weixin.qq.com/s/jHFlveANWrHX-HzJMdz9tA 参考论文:A highly conserved core bacterial microbiota with nitrogen-fixation capacity inhabits the xylem sap in maize plants.https://doi.org/10.1038/s41467-022-31113-w

# 读取数据# Load datatop_genus <- read.delim("data/top_genus.txt",                       header = TRUE,                        row.names = 1,                        sep = "\t",                       stringsAsFactors = FALSE,                       check.names = FALSE)
# 设置 Genus 因子级别# Setting Genus Factor Levelstop_genus$Genus <- factor(  top_genus$Genus,  levels = c("Bacillus","Cronobacter", "Unclassified_Enterobacterales",             "Klebsiella", "Pantoea", "Pseudomonas", "Rosenbergiella"),  labels = c("Bacillus", "Cronobacter", "Unclassified_Enterobacterales",             "Klebsiella", "Pantoea", "Pseudomonas", "Rosenbergiella"))
# 设置颜色# Setting coloursphy.cols <- c("#00b1c3", "#e274b6", "#5ba880", "#ab89dd",              "#eb9f76", "#c28f5e", "#524438")
# 普通折线图# Ordinary Line Chartp1 <- ggplot(data = top_genus,             aes(x = Compartment, y = RA, group = Genus, colour = Genus)) +  geom_point(size = 3) +  geom_line() +  scale_x_discrete(limits = c("RS", "RE", "VE", "SE", "LE", "P", "BS")) +  scale_colour_manual(values = phy.cols) +  labs(x = "Compartments", y = "Relative abundance (%)") +  theme_bw() +  theme(axis.text.x = element_text(size = 10),        axis.text.y = element_text(size = 10),        axis.title.x = element_text(size = 12),        axis.title.y = element_text(size = 12),        legend.title = element_text(size = 10),        legend.text = element_text(size = 8),        legend.position = "bottom")
# 平滑曲线图# Smoothing Curve Diagramp2 <- ggplot(data = top_genus,             aes(x = Compartment, y = RA, group = Genus, colour = Genus)) +  geom_point(size = 3) +  geom_xspline(spline_shape = -0.5) +  scale_x_discrete(limits = c("RS", "RE", "VE", "SE", "LE", "P")) +  scale_colour_manual(values = phy.cols) +  labs(x = "Compartments", y = "Relative abundance (%)") +  theme_bw() +  theme(axis.text.x = element_text(size = 10),        axis.text.y = element_text(size = 10),        axis.title.x = element_text(size = 12),        axis.title.y = element_text(size = 12),        legend.title = element_text(size = 10),        legend.text = element_text(size = 8),        legend.position = "bottom")
# 拼接图形# montagefinal_plot <- p1 + p2 + plot_layout(guides = "collect") +  plot_annotation(theme = theme(legend.position = "bottom"))
# 保存为 PDF# Save as PDFggsave("results/Plot_03.pdf", final_plot, width = 10, height = 6)

4.折线图添加置信区间Add confidence interval to line chart 

参考:https://mp.weixin.qq.com/s/MDEo_z1JCPhGvLUKu513_w 参考文献:Metabolic cross-feeding structures the assembly of polysaccharide degrading communities.https://www.science.org/doi/10.1126/sciadv.abk3076

# 加载数据# Load datadesign <- read.csv("data/test_design.csv", row.names = 1)otu <- read.csv("data/test_otu.csv", row.names = 1)df <- otu[1:20, 1:5]
# Z-score 转换# Z-score conversiondf_map <- t(scale(t(df), center = TRUE, scale = TRUE))annotation_row <- data.frame(Group = factor(rep(c("GH", "GT", "CBM", "AA", "CE"), c(4, 4, 4, 4, 4))))rownames(annotation_row) <- rownames(df_map)
# 设置颜色# Setting colourscol = list(Group = c(GH = "#eb9f76", GT = "#00b1c3", CBM = "#e274b6", AA = "#5ba880", CE = "#ab89dd"))
# 颜色优化 - 白色到绿色渐变色# Colour optimisation - white to green gradientheatmap_colors <- colorRampPalette(c("white", "#D8E0BB", "#eb9f76"))(100)
# 绘制热图,调整字体大小和边框# Drawing heat maps, adjusting font size and bordersmfs <- pheatmap(df_map,                annotation_row = annotation_row,                fontsize_number = 10, fontsize = 8,                cluster_cols = FALSE,                annotation_colors = col,                border_color = "white",                color = heatmap_colors)
# 绘制折线图# Plotting line graphsdf <- read.csv("data/test_otu.csv", row.names = 1)plot_data <- df[700:709, c(7, 12)]colnames(plot_data) <- c("Microbes", "error_bar")
# 计算误差条# Calculation error barsplot_data$error_bar <- plot_data$error_bar / 5Days <- as.data.frame(1:10)plot_data <- cbind(plot_data, Days)colnames(plot_data) <- c("Microbes", "Error_bar", "Days")
# 绘制优化折线图,增强视觉效果# Plotting optimised line charts for visual enhancementp2 <- ggplot(plot_data, aes(Days, Microbes)) +  geom_point(color = "#00b1c3", size = 4) +  geom_line(position = position_dodge(0.1), color = "#00b1c3", cex = 1.3) +  geom_ribbon(aes(ymin = Microbes - Error_bar, ymax = Microbes + Error_bar),              fill = "#00b1c3", alpha = 0.2) +  theme_bw(base_line_size = 1.2, base_rect_size = 1.2) +  labs(y = "Frequency", x = "Time (hours)") +  theme(axis.text = element_text(colour = 'black', size = 12),        axis.title = element_text(size = 14, face = "bold"))
# 组合图像并保存为 PDF# Combine images and save as PDFcombined_plot <- cowplot::plot_grid(mfs$gtable, p2, ncol = 2, labels = LETTERS[1:2])
# 保存为 PDF# Save as PDFggsave("results/Plot_04.pdf", plot = combined_plot, width = 10, height = 6)

5.分组折线图展示多个基因组的NxGrouped line chart showing Nx of multiple genomes 

参考:https://mp.weixin.qq.com/s/fjygBXghR79SkN_3r7TyvQ 参考文献:https://www.nature.com/articles/s41586-022-04808-9

# 加载数据# Load datadat <- read_excel("data/41586_2022_4808_MOESM5_ESM.xlsx",                  sheet = "Fig1b",                  skip = 1)
# 绘图# Plotp5 <- ggplot(data = dat, aes(x = Nx, y = Contig, group = Name)) +  geom_line(aes(color = Tech), size = 1.2) +  theme_classic(base_size = 14) +    labs(x = "Nx (%)", y = "Contig Length (MB)") +  scale_y_continuous(labels = function(x) {x / 1e6},                     limits = c(min(dat$Contig), max(dat$Contig))) +  scale_x_continuous(expand = expansion(mult = c(0, 0)),                     labels = function(x) {x * 100},                     breaks = seq(0, 1, 0.2)) +  scale_color_manual(values = c('HiFi' = "#00b1c3",                                'ONT' = "#c28f5e",                                'CLR' = "#e274b6"),                     name = "Platform") +  geom_point(data = dat %>% filter(Name == "SL4.0" | Name == "SL5.0"),             aes(x = Nx, y = Contig, shape = Name, fill = Name),             size = 3, color = "black") +  scale_shape_manual(values = c(21, 24)) +  scale_fill_manual(values = c("gray", "gray")) +  theme(panel.grid = element_blank(),        legend.position = c(0.8, 0.9),        legend.direction = "horizontal",        legend.title = element_text(size = 12),        legend.text = element_text(size = 10),        axis.text = element_text(size = 12),        axis.title = element_text(size = 14, face = "bold"))
# 保存为 PDF# Save as PDFggsave("results/Plot_05.pdf", plot = p5, width = 8, height = 6, units = "in")

6.折线图展示核心基因和非必须基因的数量Line chart showing the number of core genes and non-essential genes 

参考:https://mp.weixin.qq.com/s/nlku124HSEPMoka7od5jwg 参考文献:A chickpea genetic variation map based on the sequencing of 3,366 genomes.https://doi.org/10.1038/s41586-021-04066-1

# 读取数据# Load datadf<-read_excel("data/41586_2021_4066_MOESM13_ESM.xlsx")head(df)
table(df$Repeat)#> #>    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 #> 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 #>   17   18   19   20   21   22   23   24   25   26   27   28   29   30   31   32 #> 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 #>   33   34   35   36   37   38   39   40   41   42   43   44   45   46   47   48 #> 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 #>   49   50   51   52   53   54   55   56   57   58   59   60   61   62   63   64 #> 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 #>   65   66   67   68   69   70   71   72   73   74   75   76   77   78   79   80 #> 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 #>   81   82   83   84   85   86   87   88   89   90   91   92   93   94   95   96 #> 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 2258 #>   97   98   99  100 #> 2258 2258 2258 2258
# 非必须基因绘图# Optional gene mappingp61 <- ggplot() +  stat_summary(data = df, aes(x = `Number of individuals`, y = `Core-genome`, fill = "Core-genome"),               geom = "ribbon", alpha = 0.4,               fun.data = "mean_cl_boot", fun.args = list(conf.int = 0.99)) +  stat_summary(data = df, aes(x = `Number of individuals`, y = `Pan-genome`, fill = "Pan-genome"),               geom = "ribbon", alpha = 0.4,               fun.data = "mean_cl_boot", fun.args = list(conf.int = 0.99)) +  scale_fill_manual(values = c("Core-genome" = "#5e88c2", "Pan-genome" = "#db5983"),                    name = "Genomes") +  theme_minimal(base_size = 14) +  labs(x = NULL, y = "Core-genome / Pan-genome Size") +  theme(panel.grid = element_blank(),        axis.line = element_line(),        axis.text.x = element_blank(),        axis.ticks.x = element_blank(),        axis.ticks.y = element_line(),        legend.position = "bottom")
# 可选基因组(Dispensable-genome)部分# Dispensable-genome (DG) sectionp62 <- ggplot() +  stat_summary(data = df, aes(x = `Number of individuals`, y = `Dispensable-genome`, fill = "Dispensable-genome"),               geom = "ribbon", alpha = 0.4,               fun.data = "mean_cl_boot", fun.args = list(conf.int = 0.99)) +  scale_fill_manual(values = c("Dispensable-genome" = "#51afa8"),                    name = "Genomes") +  theme_minimal(base_size = 14) +  labs(x = "Number of Individuals", y = "Dispensable-genome Size") +  theme(panel.grid = element_blank(),        axis.line = element_line(),        axis.ticks.x = element_line(),        axis.ticks.y = element_line(),        legend.position = "bottom")
# 保存为 PDF# Save as PDFpdf(file = "results/Plot_06.pdf", width = 8, height = 8)(p61 / p62) + plot_layout(heights = c(2, 1)) +  plot_annotation(theme = theme(plot.margin = margin(10, 10, 10, 10)))dev.off()#> png #>   2

7.折线图给节点添加误差棒Add error bars to nodes in line chart 

参考:https://mp.weixin.qq.com/s/fsd73VyBonTALRzfpKF7LA

# 载入数据# Load datadf <- read.csv("data/test_otu.csv",row.names = 1)plot_data <- df[2021:2056,c(7,12)]colnames(plot_data) <- c("Microbes","error_bar")
# 除以10倍作为误差棒,没有实际意义# Dividing by 10 is used as the error bar, which has no practical significanceplot_data$error_bar = plot_data$error_bar/10
# 随机选择1-10天作为纵坐标时间# Randomly select 1-10 days as the vertical axis timeDays <- as.data.frame(c(1:12,1:12,1:12))plot_data <- cbind(plot_data,Days)
# 添加分组# Add groupsGroup <- as.data.frame(c(rep("CK1",12),rep("CK2",12),rep("CK3",12)))plot_data <- cbind(plot_data,Group)
# 数据命名# Name the datacolnames(plot_data) <- c("Microbes","Error_bar","Days","Group")
# 调整数据范围# Adjust the data rangeplot_data[1:12,]$Microbes <- plot_data[1:12,]$Microbes+2plot_data[13:24,]$Microbes <- plot_data[13:24,]$Microbes+2
# 添加小型柱状图# Add a small bar chartplot_data2 <- plot_data[c(8, 22, 25), ]p71 <- ggplot() + geom_bar(data = plot_data2, aes(x = Group, y = Microbes, fill = Group), position = "dodge", stat = "identity", width = 0.75) + geom_errorbar(data = plot_data2, aes(x = Group, ymin = Microbes - Error_bar, ymax = Microbes + Error_bar, color = Group), width = 0.15, size = 1.1) + scale_fill_manual(values = c("#5e88c2", "#d66356", "#31985f")) + scale_color_manual(values = c("#5e88c2", "#d66356", "#31985f")) + theme_minimal() + theme(legend.position = 'none', axis.text = element_text(colour = 'black', size = 8), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(color = "black"), # 添加黑色轴线 panel.border = element_rect(color = "black", fill = NA, size = 1))
# 绘制散点图,添加折线和平均值点,可自行调整和添加分组, annotation_custom添加左上角小图。# Draw a scatter plot, add a line and average point, you can adjust and add groups by yourself, annotation_custom adds a small picture in the upper left corner.pdf("results/Plot_07.pdf", width = 8, height = 6)ggplot(plot_data, aes(Days, Microbes, color = Group, fill = Group)) + geom_point(size = 3, shape = 21, color = "gray2") + geom_line(position = position_dodge(0.1), cex = 0.6) + geom_errorbar(aes(ymin = Microbes - Error_bar, ymax = Microbes + Error_bar), width = 0.5, cex = 0.6) + scale_color_manual(values = c("#5e88c2", "#d66356", "#31985f")) + scale_fill_manual(values = c("#5e88c2", "#d66356", "#31985f")) + labs(x = "Days", y = "Relative abundance (%)") + theme_minimal() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(color = "black"), panel.border = element_rect(color = "black", fill = NA, size = 1), axis.text = element_text(colour = 'black', size = 10)) + # 调整小图位置 # Adjusting the position of the small picture annotation_custom(grob = ggplotGrob(p71),ymin = 7, ymax = 12, xmin = 3, xmax = 8 ) dev.off()#> png #> 2

8.折线图给节点添加误差棒,并添加组间差异比较Add error bars to nodes in line chart and add comparison of differences between groups 

参考:https://mp.weixin.qq.com/s/WQkkqRA4Zx7cSnpTAVaqPg 参考文献:Gut microbiota modulates weight gain in mice after discontinued smoke exposure.https://doi.org/10.1038/s41586-021-04194-8

# 载入数据# Load datadata <- read.csv('data/2022.1.15.CSV',header = T)head(data)  #View(data)
# 数据整理# Data processingdata1 <- data[,c(1:5)] %>%  reshape2::melt(id=c('treat','Day'))
# 设置绘制误差棒errorbar时用到的函数# Set the function used to draw the error bar errorbartopbar <- function(x){  return(mean(x)+sd(x)/sqrt(length(x))) #误差采用了mean+-sem(The error is expressed as mean+-sem)}
bottombar <- function(x){  return(mean(x)-sd(x)/sqrt(length(x)))}
# 绘图# Plotp81 <- ggplot(data1, aes(Day, value, color = treat)) +  geom_rect(aes(xmin = 21, xmax = 40, ymin = (-Inf), ymax = Inf),            fill = 'grey90', color = 'grey90') +  geom_vline(xintercept = 21, linetype = 2, cex = 1.2) +  stat_summary(geom = 'line', fun = 'mean', cex = 2.5) +  stat_summary(geom = 'errorbar',               fun.min = bottombar, fun.max = topbar,               width = 1, cex = 0.8, aes(color = treat)) +  stat_summary(geom = 'point', fun = 'mean', aes(fill = treat),               size = 5, pch = 21, color = 'black') +  theme_classic(base_size = 15) +  theme(legend.position = 'none')
# 折线图添加细节# Add details to the line chartp82 <- p81 +  scale_color_manual(values = c('#5494cc', '#0d898a', '#e18283', '#f9cc52')) +  scale_fill_manual(values = c('#5494cc', '#0d898a', '#e18283', '#f9cc52')) +  scale_y_continuous(breaks = seq(0, 60, 20), expand = c(0, 0)) +  scale_x_continuous(breaks = seq(0, 40, 10), expand = c(0, 0)) +  labs(y = 'Weight change(%)') +  theme(axis.line = element_line(size = 1),        axis.text = element_text(color = 'black'),        axis.ticks = element_line(size = 1, color = 'black')) +  ggplot2::annotate(geom = 'segment', x = 36.2, xend = 36.2, y = 18, yend = 26, cex = 1.2) +  ggplot2::annotate(geom = 'text', label = '***', x = 37.5, y = 22, size = 7, angle = 90) +  ggplot2::annotate(geom = 'segment', x = 38, xend = 38, y = 18, yend = 40, cex = 1.2) +  ggplot2::annotate(geom = 'text', label = '****', x = 39.5, y = 29, size = 7, angle = 90)
# 添加图例(add legend)# 折线图最终效果(The final effect of the line chart)linechart <- p82 +  coord_cartesian(clip = 'off', ylim = c(0, 60), xlim = c(0, 40)) +  theme(plot.margin = margin(1.5, 0.5, 0.5, 0.5, 'cm')) +  geom_rect(aes(xmin = 0, xmax = 3, ymin = 63, ymax = 65), fill = '#00bae1', color = 'black') +  geom_rect(aes(xmin = 10, xmax = 13, ymin = 63, ymax = 65), fill = '#4bb395', color = 'black') +  geom_rect(aes(xmin = 20, xmax = 23, ymin = 63, ymax = 65), fill = '#f47aab', color = 'black') +  geom_rect(aes(xmin = 30, xmax = 33, ymin = 63, ymax = 65), fill = '#c29c5e', color = 'black') +  ggplot2::annotate('text', x = 6, y = 64, label = 'Non-SMK', size = 4) +  ggplot2::annotate('text', x = 15.2, y = 64, label = 'SMK', size = 4) +  ggplot2::annotate('text', x = 26.4, y = 64, label = 'Non-SMK+\nabx', size = 4) +  ggplot2::annotate('text', x = 36, y = 64, label = 'SMK+abx', size = 4)
# 添加柱状图(Add a bar chart)# 数据处理(data processing)data2 <- data[,6:9] %>%  gather(key = treat)
data3 <- data[,10:13] %>%  gather(key = treat)
# 左边柱状图(Left column chart)leftchart <- ggplot(data2, aes(factor(treat, levels = c('SMK.abx', 'Non_SMK.abx', 'SMK', 'Non_SMK')), value)) +  stat_summary(geom = 'bar', fun = 'mean', fill = 'white', color = 'black', width = 0.7, cex = 1) +  stat_summary(geom = 'errorbar', fun.min = bottombar, fun.max = topbar, width = 0.3, cex = 0.8, color = 'black') +  geom_jitter(aes(color = factor(treat, levels = c('SMK.abx', 'Non_SMK.abx', 'SMK', 'Non_SMK'))),              width = 0.1, size = 1.5) +  scale_color_manual(values = c('#c29c5e', '#4bb395', '#f47aab', '#00bae1')) +  labs(x = NULL, y = NULL) +  scale_y_continuous(limits = c(-40, 600), expand = c(0, 0)) +  geom_hline(yintercept = 0, cex = 1) +  theme_classic(base_size = 15) +  theme(axis.ticks.y = element_blank(),        axis.text.y = element_blank(),        legend.position = 'none',        axis.line = element_line(size = 1),        axis.text = element_text(color = 'black'),        axis.ticks = element_line(size = 1, color = 'black')) +  coord_flip()
# 右边柱状图(Right bar graph)rightchart <- ggplot(data3, aes(factor(treat, levels = c('SMK.abx.1', 'Non_SMK.abx.1', 'SMK.1', 'Non_SMK.1')), value)) +  stat_summary(geom = 'bar', fun = 'mean', fill = 'white', color = 'black', width = 0.7, cex = 1) +  stat_summary(geom = 'errorbar', fun.min = bottombar, fun.max = topbar, width = 0.3, cex = 0.8, color = 'black') +  geom_jitter(aes(color = factor(treat, levels = c('SMK.abx.1', 'Non_SMK.abx.1', 'SMK.1', 'Non_SMK.1'))),              width = 0.1, size = 1.5) +  scale_color_manual(values = c('#c29c5e', '#4bb395', '#f47aab', '#00bae1')) +  labs(x = NULL, y = NULL) +  scale_y_continuous(limits = c(-40, 500), expand = c(0, 0)) +  geom_hline(yintercept = 0, cex = 1) +  theme_classic(base_size = 15) +  theme(axis.ticks.y = element_blank(),        axis.text.y = element_blank(),        legend.position = 'none',        axis.line = element_line(size = 1),        axis.text = element_text(color = 'black'),        axis.ticks = element_line(size = 1, color = 'black')) +  coord_flip()
# 叠加绘图# Added plotfinal_plot <- linechart +  annotation_custom(ggplotGrob(leftchart), xmin = 0, xmax = 20.5, ymin = 40, ymax = 57) +  annotation_custom(ggplotGrob(rightchart), xmin = 21, xmax = 39.5, ymin = 40, ymax = 57) +  ggplot2::annotate('text', label = 'iAUC: Exposure', x = 10.5, y = 58.5, size = 7) +  ggplot2::annotate('text', label = 'iAUC: Cessation', x = 31, y = 58.5, size = 7)
# 保存为PDF# Save plotggsave('results/Plot_08.pdf', final_plot, width = 8, height = 6)

使用此脚本,请引用下文:

Yong-Xin Liu, Lei Chen, Tengfei Ma, Xiaofang Li, Maosheng Zheng, Xin Zhou, Liang Chen, Xubo Qian, Jiao Xi, Hongye Lu, Huiluo Cao, Xiaoya Ma, Bian Bian, Pengfan Zhang, Jiqiu Wu, Ren-You Gan, Baolei Jia, Linyang Sun, Zhicheng Ju, Yunyun Gao, Tao Wen, Tong Chen. 2023. EasyAmplicon: An easy-to-use, open-source, reproducible, and community-based pipeline for amplicon data analysis in microbiome research. iMeta 2: e83. https://doi.org/10.1002/imt2.83

Copyright 2016-2024 Defeng Bai baidefeng@caas.cn, Chuang Ma 22720765@stu.ahau.edu.cn, Jiani Xun 15231572937@163.com, Yong-Xin Liu liuyongxin@caas.cn

宏基因组推荐
本公众号现全面开放投稿,希望文章作者讲出自己的科研故事,分享论文的精华与亮点。投稿请联系小编(微信号:yongxinliu 或 meta-genomics)

猜你喜欢

iMeta高引文章 fastp 复杂热图 ggtree 绘图imageGP 网络iNAP
iMeta网页工具 代谢组MetOrigin 美吉云乳酸化预测DeepKla
iMeta综述 肠菌菌群 植物菌群 口腔菌群 蛋白质结构预测

10000+:菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature

系列教程:微生物组入门 Biostar 微生物组  宏基因组

专业技能:学术图表 高分文章 生信宝典 不可或缺的人

一文读懂:宏基因组 寄生虫益处 进化树 必备技能:提问 搜索  Endnote

扩增子分析:图表解读 分析流程 统计绘图

16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun

生物科普:  肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流快速解决科研困难,我们建立了“宏基因组”讨论群,己有国内外6000+ 科研人员加入。请添加主编微信meta-genomics带你入群,务必备注“姓名-单位-研究方向-职称/年级”。高级职称请注明身份,另有海内外微生物PI群供大佬合作交流。技术问题寻求帮助,首先阅读《如何优雅的提问》学习解决问题思路,仍未解决群内讨论,问题不私聊,帮助同行。

点击阅读原文



宏基因组
宏基因组/微生物组是当今世界科研最热门的研究领域之一,为加强本领域的技术交流与传播,推动中国微生物组计划发展,中科院青年科研人员创立“宏基因组”公众号,目标为打造本领域纯干货技术及思想交流平台。
 最新文章