该整理代码封个R包了-来自复兴号的学习笔记

本次基于扩增子的培训代码全部的环境在terminal中。这将在无论是个体电脑还是服务器上都可以顺利布置，更加高效，灵活，把握核心内容。我目前的习惯还是基于rstudio的Rmd。我考虑到的是可视化工作需要不断修改代码以及生成报告。但是分析工作需要面对的代码必然会更多。

今天启程南京，在高铁上有将近五个小时时间，闲来无聊，我来总结一下并且学点东西。

为了进一步提高自己对于平台和分析工具的使用，必须加快对这几项技能的学习

第一项技能：R包构建：

由于之前尚未对许多自编功能写成函数，当然更没有写成R包，显然随着分析工作内容的增加，我面对的代码早已冗余繁杂，为了更加精简代码，让思路游走于代码，我必须将这些个性化的功能编写成函数，不断的加入到R包中。今天我就来创建这个R包，相信随着代码的增加，这个R包也可以帮到大家。

刘老师在这次培训中使用了一个R包，在win下的分析流程运行完成之后，我便开始在Ubuntu上布置分析流程，这个名为amplicon的R包确实装了好半天。以这个R包为模板，我开始创建。

第一步：R包成型

打开rdstudio，创建新的项目，选择构建R包，选择自己要保存的位置，这样就可以开始R包的制作了。

我的目的很明确：当我载入R包的时候可以使用我常用的那几个功能即可。

我在项目中的R文件夹中添加R脚本即可。当然我没有检查依赖，将功能部署完成，并且在功能前面添加载入R包即可。反正我的都安装了这些R包了。

点击打包R包并安装，系统自动重新启动R

输入代码，进行测试 这个函数的目的是对不同研究的微生物群落数据进行合并，有多种模式可选。这是方便将不同研究群落一起分析的一个有用的函数，是基于phyloseq对象编写的。后面附有源代码；

library(easyMicrobiome)

ps1 = readRDS("./data/ps1.rds")
ps2 = readRDS("./data/ps2.rds")
ps_merge = merge_ps(ps1,ps2,model = 1)

ps_merge

当我输入

?merge_ps

提示找不到对应的文件。我没有编写帮助文件哎。弄一下

帮助文件编写

我为merge_ps函数编写了一个Rd文件，保存到man文件夹中。

\name{merge_ps}
\alias{merge_ps}
\title{merge_ps}
\usage{

}
\description{
当前我的工作集中搭建在phyloeq文件上，为了将不同的OTU表格进行合并分析，我手动写了这个功能。基于phyloseq对象合并OTU表格，tax表格和mapping表格。

共有三种模式
有s三种模式，model =1 是取共有部分进行合并
model =2 是保留x中的全部OTU，y中如果有其他OTU，全部舍弃，如果没有x中的这些otu使用0填充。
注意进化树可不能合并，必须重新运算，这里没有必要花费很长的时间来运算这个。
model = 3,xy全部OTU进行合并

mapping文件：两个ps对象mapping文件必须包含SampleTypeeh和 Desciption列，其他内容随便。

}
\examples{
#example1
ps1 = data(/data/ps1.rds)
ps2 = data(/data/ps2.rds)

merge_ps = (ps1,ps2,model = 1)

}

这次运行？merge_ps就可以成功了，自己用中文写的帮助文件

如果使用帮助文件，那么我就要内置数据集，方便调用数据。虽然我建立了一个文件件名为data，里面保存了两份phyloseq的对象。但是我载入这个包的时候却没有自动载入这个数据集

对easyMicrobiome包添加数据集

刚开始我直接将两个phyloseq对象（后缀名：.rds）添加到data文件中，但是无法载入。这个时候刘老师的amplicon包中对应位置的文件却都是RDA格式的，这是r语言的环境，使用save保存的对象。

于是我做了一下改动

ps1 = readRDS("D:/Shared_Folder/my_R_packages/easy_microbiome/easyMicrobiome/data/ps1.rds")
ps2 = readRDS("D:/Shared_Folder/my_R_packages/easy_microbiome/easyMicrobiome/data/ps2.rds")
save(ps2, file = "D:/Shared_Folder/my_R_packages/easy_microbiome/easyMicrobiome/data/ps2.rda")
save(ps1, file = "D:/Shared_Folder/my_R_packages/easy_microbiome/easyMicrobiome/data/ps1.rda")

重新打包R包后再次运行。

这次R包就带有数据了。

library(easyMicrobiome)
merge_ps(ps1,ps2,model = 1)

我注意到我的R包帮助文件好丑，特别是和刘老师的进行一个比对，发现特别丑。缺少参数解释，选项解释.

完善?merge_ps帮助文件

参照刘老师帮助文件，我使用中文编写了一下

\name{merge_ps}
\alias{merge_ps}
\title{merge_ps}
\usage{
merge_ps (ps1,ps2,model = 1)
}
\arguments{
\item{ps1}{两个合并phyloseq对象中的一个，这个对象至少包含OTU表格，mapping文件，注释TAX文件;}

\item{ps2}{两个合并phyloseq对象中的第二个，注意有顺序，这个对象至少包含OTU表格，mapping文件，注释TAX文件;}

\item{model}{这是合并模式的选择，可以选择1，2，3，这三种模式;}

}
\value{
model 1,2,3.分别代表了三种合并方式，model =1 是取共有部分进行合并
；model =2 是保留x中的全部OTU，y中如果有其他OTU，全部舍弃，如果没有x中的这些otu使用0填充。
；model = 3,xy全部OTU进行合并。注意进化树可不能合并，必须重新运算，这里没有必要花费很长的时间来运算这个。
}

\description{
当前我的工作集中搭建在phyloeq文件上，为了将不同的OTU表格进行合并分析，我手动写了这个功能。基于phyloseq对象合并OTU表格，tax表格和mapping表格。

共有三种模式
有s三种模式，model =1 是取共有部分进行合并
model =2 是保留x中的全部OTU，y中如果有其他OTU，全部舍弃，如果没有x中的这些otu使用0填充。
注意进化树可不能合并，必须重新运算，这里没有必要花费很长的时间来运算这个。
model = 3,xy全部OTU进行合并

mapping文件：两个ps对象mapping文件必须包含SampleTypeeh和 Desciption列，其他内容随便。

}
\examples{
#example1
merge_ps(ps1,ps2,model = 1)
#example2
merge_ps(ps1,ps2,model = 2)

#example3
merge_ps(ps1,ps2,model = 3)

}

好像现在就像回事了，一个R包就是这样，多少个R脚本，多少个帮助文件。

小结

之前学习过编写R包的相关知识，但是一直没有主动，建一个R包，本次借鉴刘老师的R包amplicon，我在回南京的途中实战熟悉了一下R包的构建流程。比较不上档次的是整个文件注释部分都是用中文编写的，我仔细想了想，R包似乎都是 英文编写的，中文的很少。

这个包可以安装了

运行下面代码安装library("devtools")
install_github("taowenmicro/easyMicrobiome")


源代码

merge_ps = function(ps1,ps2,model = 1){
library("phyloseq")
library(tidyverse)

vegan_otu <- function(physeq){
OTU <- otu_table(physeq)
if(taxa_are_rows(OTU)){
OTU <- t(OTU)
}
return(as(OTU,"matrix"))
}

vegan_tax <- function(physeq){
tax <- tax_table(physeq)

return(as(tax,"matrix"))
}

#otu_table
if(model == 1){
otu_table1 = as.data.frame(t(vegan_otu(ps1)))
dim(otu_table1)
otu_table2 = as.data.frame(t(vegan_otu(ps2)))
dim(otu_table2)
otu_table3 = merge(otu_table1,otu_table2,by = "row.names",all = T)
dim(otu_table3)
row.names(otu_table3) = otu_table3$Row.names
otu_table3$Row.names = NULL
otu_table3 = as.matrix(otu_table3)
otu_table3[is.na(otu_table3)] <- 0
#tax_table

tax_table1 = as.data.frame((vegan_tax(ps1)))
dim(tax_table1)
tax_table2 = as.data.frame((vegan_tax(ps2)))
dim(tax_table2)
head(tax_table1)
head(tax_table2)
tax_table3 = rbind(tax_table1,tax_table2)
dim(tax_table3)

}
if(model == 2){
otu_table1 = as.data.frame(t(vegan_otu(ps1)))
dim(otu_table1)
otu_table2 = as.data.frame(t(vegan_otu(ps2)))
dim(otu_table2)
otu_table3 = merge(otu_table1,otu_table2,by = "row.names",all.x = T)
dim(otu_table3)
row.names(otu_table3) = otu_table3$Row.names
otu_table3$Row.names = NULL
otu_table3 = as.matrix(otu_table3)
otu_table3[is.na(otu_table3)] <- 0
#tax_table

tax_table1 = as.data.frame((vegan_tax(ps1)))
dim(tax_table1)
tax_table3 = tax_table1
}
if(model == 3){
otu_table1 = as.data.frame(t(vegan_otu(ps1)))
dim(otu_table1)
otu_table2 = as.data.frame(t(vegan_otu(ps2)))
dim(otu_table2)
otu_table3 = merge(otu_table1,otu_table2,by = "row.names",all = T)
dim(otu_table3)
row.names(otu_table3) = otu_table3$Row.names
otu_table3$Row.names = NULL
otu_table3 = as.matrix(otu_table3)
otu_table3[is.na(otu_table3)] <- 0
#tax_table
tax_table1 = as.data.frame((vegan_tax(ps1)))
dim(tax_table1)
tax_table2 = as.data.frame((vegan_tax(ps2)))
dim(tax_table2)
head(tax_table1)
head(tax_table2)
tax_table3 = rbind(tax_table1,tax_table2)
dim(tax_table3)
}

mapping1 = as.data.frame(sample_data(ps1))
head(mapping1)
mapping1 = mapping1[,c("SampleType","Description")]
mapping2 = as.data.frame(sample_data(ps2))
head(mapping2)
# mapping2$BarcodeSequence = NULL
# mapping2$LinkerPrimerSequence = NULL
mapping2$SampleType = paste(mapping2$fianl_SampleType,mapping2$zone,sep = "_")
mapping2 = mapping2[,c("SampleType","Description")]
mapping3 = rbind(mapping1,mapping2)
head(mapping3)

ps_add_out =phyloseq(otu_table(as.matrix(otu_table3),taxa_are_rows = T),
sample_data(mapping3),
tax_table(as.matrix(tax_table3)))
return(ps_add_out)
}