编程模板-R语言脚本写作：最简单的统计与绘图，包安装、命令行参数解析、文件读取、表格和矢量图输出

不常用R的朋友，可能并不能意识到好的编辑模板带来的长期效率提高和代码重用性，非计算机专业科研人员且需要使用R语言，此文会对您有帮助。

个人认为：是否能熟悉使用Shell(项目流程搭建)+R(数据统计与可视化)+Perl/Python等(胶水语言，数据格式转换，软件间衔接)三门语言是一位合格生物信息工程师的标准。

之前分享过我个人的《Shell语言》和《Perl语言》脚本写作模板(蓝色字为链接直达)，今天再分享一下我的R语言模板，一次性解决困扰新手的众多问题，如包安装、命令行参数解析、文件读取、Anova组间统计和箱线图展示、表格和矢量图输出等。

希望对新手有所帮助！老司机尽情拍砖，指导我写出更好用、更规范的脚本。

模板主要内容

此模板主要以下分为5部分：

1. 程序功能描述和分析思路

每次写脚本，一定要写清楚程序的功能、实现的主要步骤，每个参数的详细说明和使用示例。如果英语描述不方便理解，建议中英文全写，方便自己和同行快速理解。千万要认识写清楚这些，现点养成好习惯，将来节省的时间会更多。

2. 依赖关系检查、安装和加载

R己有过万的可用包，很少需要我们开发新功能。每次需要加载需要的包，但没有会报错，直接安装又不会提醒导致重复安装浪费时间。我添加判断R包是否安装的代码，当R包缺失时会自己动从CRAN安装。对于CRAN没有的包，也提供了Bioconductor和Github安装代码，只需手动修改包名安装即可。

此步还包括加载了optparse包后，对命令行的解析，以及在屏幕上显示参数，让用户校验是否正确。

3. 读取输入文件

R用于作图，读入最多的是制表符分割的表格。我提供了使用rnorm生成的测试数据，并输入了测试文件，其实这种方式在程序测试、教程和实验分析对照中也是很常用的。

同时提供了read.table读取文件，或read.table(file.choose())弹窗选择输入文件两种方式。

4. 统计与绘图

这其实才是正文部分，用户根据自己的需求编写代码分析数据，完成工作。本模板以anova组间比较，和箱线图展示为例，方法大家测试工作过程。

表1. Anova分析测试数据结果

图1. 箱线图展示组间比较，且按组分配不同的颜色和形状(方便多组时区分)，并添加统计p值。

5. 保存图表

分析的结果最重要的工作是保存图表，用于进一步分析和论文写作。本模板提供了输出表格的代码，并解决了存在行名输出时，文本文件在Excel中打开列名位置串一格的问题(R表输出时行/列名交叉点为空)。图片建议一律输出PDF格式，主要控制图片的长宽即可，不存辨率的问题，后期AI修改组合非常方便。杂志也更喜欢，因为他们接收后修改方便。

模板使用方法

保存文末的全部代码，推荐在Rstudio中保存文件名为template.r，可以逐行运行学习，体会常用功能。

命令行中使用

# 方法1. 直接使用Rscript执行
Rscript template.r

# 调置输入数据文件，输出图表文件名前缀
Rscript template.r -i data_table.txt -o output

# 方法2. 添加可执行权限和环境变量
chmod +x template.r
# 链接至有权限写入的环境变量目录，如我的是~/bin/
ln `pwd`/template.r ~/bin/
# 任何地方可以直接调用此脚本
template.r -i data_table.txt -o output

程序运行结果如下：

[1] "The input file is data_table.txt"
[1] "The output file prefix is output"
[1] "The output table is output.txt"
[1] "The output figure is output.pdf"

文中把每个段落用if条件语句控制是否执行。方便读者开、关功能代码段。

用户还可根据自己的需求添加段落，再也不要为找常用语句的写法到处搜索了，有了模板，常用语句段落直接复制、修改，快速实现自己的工作！！！

模板全部代码及讲解

#!/usr/bin/env Rscript

# 1. 程序功能描述和主要步骤

# 程序功能：Anova组间统计和箱线图展示
# Anova Script functions: Calculate pvalue of groups by aov and TukeyHSD
# Main steps: 
# - Reads data table input.txt
# - Calculate pvalue and save in output.txt
# - Draw boxplot and save in output.pdf

# 程序使用示例
# USAGE
# Default
# anova.r     -i data_table.txt
#                        -o otuput filename prefix for output directory name 

# 参数说明
# Options
# -i/--input     输入数据表文件 input.txt
# -o/--output     输出结果文件名前缀 output_prefix, 通常会有统计表txt和矢量图pdf
options(warn = -1)

# 2. 依赖关系检查、安装和加载
# See whether these packages exist on comp. If not, install.
package_list <- c("optparse","reshape2","ggplot2","ggpubr")

for(p in package_list){
  if(!suppressWarnings(suppressMessages(require(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))){
    install.packages(p, repos="http://cran.r-project.org")
    suppressWarnings(suppressMessages(library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))
  }
}
# 另两种常见R包安装方法
if (FALSE){
  # Bioconductor安装
  source("https://bioconductor.org/biocLite.R")
  biocLite(c("reshape2"))
  # Github安装
  install.packages("devtools", repo="http://cran.us.r-project.org")
  library(devtools)
  install_github("kassambara/ggpubr")
}

# 清理工作环境 clean enviroment object
rm(list=ls()) 

# 加载依赖关系 Load essential packages
library(optparse)
library(reshape2)
library(ggplot2)
library(ggpubr)

# 解析命令行
if (TRUE){
  option_list <- list(
  make_option(c("-i", "--input"), type="character", default="input.txt",
              help="Input table file to read [default %default]"),
  make_option(c("-o", "--output"), type="character", default="output",
              help="output directory or prefix [default %default]")
)
opts <- parse_args(OptionParser(option_list=option_list))

# 显示输入输出确认是否正确
print(paste("The input file is ", opts$input,  sep = ""))
print(paste("The output file prefix is ", opts$output, sep = ""))
}

# 3. 读取输入文件
# 需要使用哪种方式，将其设置为TRUE，其它为FALSE即可

# 产生两组数据比较
if (TRUE){
  # 产生两种各10个正态分布数值，分别以均值为1和2，标准差为均值的0.5倍，colnames修改组名称
  set.seed(123)
  A = rnorm(10, mean = 1, sd = 0.5)
  B = rnorm(10, mean = 2, sd = 1)
  dat=as.data.frame(cbind(A,B))
  colnames(dat)=c("GroupA","GroupB")
  # 保存数据文件方便下面演示
  write.table(dat, file="opts$input", append = F, quote = F, sep="\t", eol = "\n", na = "NA", dec = ".", row.names = F, col.names = T)
}

# 从文件中读取
if (FALSE){
  dat = read.table("opts$input", header=T, row.names = NULL, sep="\t")
}

# 弹出窗口选择文件
if (FALSE){
  dat = read.table(file.choose(), header=T, row.names = NULL, sep="\t")
}

# 4. 统计与绘图
if (TRUE){
  # 将宽表格将换为长表格(变量位于同列方便操作)
  dat_melt=melt(dat, id.vars = 0)
  # anova比较各组variable的观测值value
  anova <- aov(value ~ variable, data = dat_melt)
  # 计算Tukey显著性差异检验
  Tukey_HSD <- TukeyHSD(anova, ordered = TRUE, conf.level = 0.95)
  # 提取比较结果
  Tukey_HSD_table <- as.data.frame(Tukey_HSD$variable) 
  # 展示组间差异：四列分别为均值差异，差异最小值，差异最大值和校正的P值
  Tukey_HSD_table
  # 采用ggpubr包中的ggboxplot展示组间比较箱线图，并添加Wilcoxon整体差异P值
  # 增加了jitter点，color和shape按分组variable变化，添加统计值 
  p <- ggboxplot(dat_melt, x="variable", y="value", color = "variable", 
                 add = "jitter", shape="variable") +stat_compare_means()
  p
}

# 5. 保存图表
if (TRUE){
  # 保存一个制表符，解决存在行名时，列名无法对齐的问题
  write.table("\t", file=paste(opts$output,".txt",sep=""),append = F, quote = F, eol = "", row.names = F, col.names = F)
  # 保存统计结果，有waring正常
  write.table(Tukey_HSD_table, file=paste(opts$output,".txt",sep=""), append = T, quote = F, sep="\t", eol = "\n", na = "NA", dec = ".", row.names = T, col.names = T)
  print(paste("The output table is ", opts$output, ".txt",  sep = ""))
  # 保存图片至文件，pdf方便AI修改成出版级图片
  ggsave(file=paste(opts$output,".pdf",sep=""), p, width = 5, height = 3)
  print(paste("The output figure is ", opts$output, ".pdf",  sep = ""))
}

猜你喜欢

10000+：菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature Cell专刊 肠道指挥大脑

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

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

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

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

16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板: Shell  R Perl

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

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外2400+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读