12.1 创建和运行一个脚本

12.1.1 问题

你想要创建和运行一个脚本。

12.1.2 方案

R 脚本通常是以 .R 为文件拓展名的纯文本文件。因此创建R脚本可以是任意的文本编辑器，R 最佳的集成开发环境是 RStudio，它开源并有免费的桌面版本和服务器版本发布，推荐读者下载、安装和使用。

在类 Unix 系统中，除了编辑器，我们还可以使用终端命令。例如，创建一个输出 Hello world! 的 R 脚本。

  1. echo 'cat("Hello world!")' > test.R

运行 R 脚本可以在 R 控制台使用 source()函数。

  1. source("test.R", print.eval = TRUE)
  2. #> Hello world!

也可以在终端中使用 RScript 执行。

  1. RScript test.R
  2. #> Using library: /Users/wsx/R_library
  3. #> 载入需要的程辑包：pacman
  4. #> Hello world!

12.1.3 案例：计算细菌基因组核心蛋白相似性

  • 应用场景

细菌分类学研究中，需要借助基因组水平的相似度来界定是否属于新物种，是否是一个未发现的新属水平或者新科水平，乃至更高的分类学单元（界/门/纲/目/科/属/种）。

在基因组的核酸水平研究中，有诸如 dDDH（数字化 DNA 分子杂交）、核苷酸平均相似度（Average Nucleotide Identity，ANI）等指标来界定是否属于新物种；而在基因组蛋白质水平相类似的指标较少，比如氨基酸平均相似度（Average Amino acid Identity，AAI）和保守蛋白比率（percentage of conserved proteins，POCP）等。

  • 简要过程

两两比对细菌基因组的蛋白序列，互为参考数据库进行 blastp 比对（A 作数据库，B 查询；B 作数据库，A 查询），数据筛选的标准是：一致度大于 40%，查询片段的长度大于原片段长度的 50%，e 值小于 1e-5。

  • 参考文献

Qin, Q. L., Xie, B. B., Zhang, X. Y., Chen, X. L., Zhou, B. C., Zhou, J., … & Zhang, Y. Z. (2014). A proposed genus boundary for the prokaryotes based on genomic insights. Journal of bacteriology, 196(12), 2210-2215.

以下 R 脚本都是在 windows 操作平台上进行的。

  1. # 下载所分析的基因组数据（蛋白序列）
  2. # 存放于 Rawdata 文件夹中
  3. if (!dir.exists('Rawdata')) {
  4.   dir.create('Rawdata')
  5. }
  7. # 示例-1: Pseudomonas aeruginosa
  8. # ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/006/765/GCF_000006765.1_ASM676v1/GCF_000006765.1_ASM676v1_protein.faa.gz
  10. # 示例-2: Acinetobacter baumannii
  11. # ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/746/645/GCF_000746645.1_ASM74664v1/GCF_000746645.1_ASM74664v1_protein.faa.gz
  13. # 使用 R.utils 中的 gunzip 解压缩
  14. library(R.utils)
  16. download.file("ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/006/765/GCF_000006765.1_ASM676v1/GCF_000006765.1_ASM676v1_protein.faa.gz",
  17.               destfile = "Rawdata/Pseudomonas_aeruginosa.faa.gz")
  18. gunzip("Rawdata/Pseudomonas_aeruginosa.faa.gz")
  20. download.file("ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/746/645/GCF_000746645.1_ASM74664v1/GCF_000746645.1_ASM74664v1_protein.faa.gz",
  21.               destfile = "Rawdata/Acinetobacter_baumannii.faa.gz")
  22. gunzip("Rawdata/Acinetobacter_baumannii.faa.gz")
  1. # 使用 dbplyr 对数据框中的某列去重复
  2. library(dbplyr)
  3. # 使用 seqinr 格式化 fasta 格式的序列
  4. library(seqinr)
  6. # 检查存放中间过程文件的文件夹是否存在
  7. if (!dir.exists('Database')) {
  8.   dir.create('Database')
  9. }
  11. if (!dir.exists('Result')) {
  12.   dir.create('Result')
  13. }
  15. # 获取所有待分析基因组文件名
  16. genome.files <- list.files('Rawdata')
  18. # 对所有的待分析基因组建库
  19. for (gn in genome.files) {
  20.   header.file <- strsplit(gn,'.',fixed = T)[[1]][1]
  21.   commond.makedb <- paste0('diamond.exe makedb --in Rawdata/',
  22.                            gn, ' --db Database/', header.file)
  23.   system(commond.makedb)
  24. }
  26. # 获取多基因组的两两比对的组合数据集
  27. genome.comn <- combn(genome.files,2)
  29. # 计算核心蛋白相似性的骨架命令
  30. blast.comm1 <- 'diamond.exe blastp -q Rawdata/'
  31. blast.comm2 <- ' -d Database/'
  32. blast.comm3 <- ' -e 1e-5 --id 40 -o Result/'
  34. # 建立新变量，保存运算结果
  35. pocp.vector <- c()
  37. for (i in (1:dim(genome.comn)[2]) ) {
  38.   a.genome <- genome.comn[,i][1]
  39.   b.genome <- genome.comn[,i][2]
  41.   a.header <- strsplit(a.genome,'.',fixed = T)[[1]][1]
  42.   b.header <- strsplit(b.genome,'.',fixed = T)[[1]][1]
  44.   a.genome.seq <- read.fasta(paste0('Rawdata/', a.genome),'AA')
  45.   b.genome.seq <- read.fasta(paste0('Rawdata/', b.genome),'AA')
  47.   a.total <- length(a.genome.seq)
  48.   b.total <- length(b.genome.seq)
  50.   str(a.genome.seq)
  51.   str(b.genome.seq)
  53.   a.seq.list <- names(a.genome.seq)
  54.   b.seq.list <- names(b.genome.seq)
  55.   a.seq.length <- c()
  57.   for (nm in a.seq.list) {
  58.     tmp.len <- length(a.genome.seq[[which(a.seq.list == nm)]])
  59.     a.seq.length <- append(a.seq.length, tmp.len)
  60.   }
  61.   b.seq.length <- c()
  63.   for (nm in b.seq.list) {
  64.     tmp.len <- length(b.genome.seq[[which(b.seq.list == nm)]])
  65.     b.seq.length <- append(b.seq.length, tmp.len)
  66.   }
  68.   a.seq.df <- data.frame(a.seq.list, a.seq.length)
  69.   colnames(a.seq.df) <- c('V1','length')
  70.   b.seq.df <- data.frame(b.seq.list, b.seq.length)
  71.   colnames(b.seq.df) <- c('V1','length')
  73.   print(paste0('-- Blasting: ',a.header,' - VS - ',b.header))
  75.   # 「正向」-- A 为查询，B 为参考数据库
  76.   result.forward <- paste0(a.header,'_VS_',b.header,'.tab')
  77.   system(paste0(blast.comm1, a.genome,
  78.                 blast.comm2, b.header,
  79.                 blast.comm3, result.forward))
  80.   df.forward <- read.table(paste0('Result/',result.forward),
  81.                            header = F,sep = '\t',
  82.                            stringsAsFactors = F)
  83.   df.forward <- df.forward  %>%  distinct(V1,.keep_all = T)
  84.   df.forward <- merge(df.forward, a.seq.df, by = 'V1', all.x = T)
  85.   df.forward$align <- df.forward$V4 / df.forward$length
  86.   df.forward <- df.forward[which(df.forward$V3 > 40 & df.forward$align > 0.5 & df.forward$V11 < 1e-5),]
  87.   C1 <- dim(df.forward)[1]
  89.   # 「反向」-- B 为查询，A 为参考数据库
  90.   result.backward <- paste0(b.header,'_VS_',a.header,'.tab')
  91.   system(paste0(blast.comm1, b.genome,
  92.                 blast.comm2, a.header,
  93.                 blast.comm3, result.backward))
  94.   df.backward <- read.table(paste0('Result/',result.backward),
  95.                            header = F,sep = '\t',
  96.                            stringsAsFactors = F)
  97.   df.backward <- df.backward %>% distinct(V1,.keep_all = T)
  98.   df.backward <- merge(df.backward, b.seq.df, by = 'V1', all.x = T)
  99.   df.backward$align <- df.backward$V4 / df.backward$length
  100.   df.backward <- df.backward[which(df.backward$V3 > 40 & df.backward$align > 0.5 & df.backward$V11 < 1e-5),]
  101.   C2 <- dim(df.backward)[1]
  103.   pocp <- (C1 + C2)/(a.total + b.total)
  104.   pocp.vector <- append(pocp.vector, paste0(a.header,'\t',b.header,'\t',pocp))
  106.   print(paste0('-- Pair blast done: ',a.header,' - VS - ',b.header))
  107.   print(paste0('-- The POCP : ', pocp))
  108.   print('----------------------------------')
  109. }
  111. write(pocp.vector, 'resultPOCP.txt')
  113. # 删除分析过程中的冗余文件
  114. unlink("Database", recursive = TRUE)
  115. unlink("Result", recursive = TRUE)
  117. # 重新创建新文件夹
  118. dir.create('Database')
  119. dir.create('Result')

12.1.3.1 提示

更多关于 POCP 计算的相关技巧，请点击这里阅读。