在易生信的某ppt中,出现了标准两列注释文件和八列注释文件,八列注释文件在之后某一步用到,但是ppt并未给出从usearch标准注释文件otus.sintax转换到标准八列注释文件的代码,所以自己尝试手动写了一下
但是shell命令不熟悉,就用R写了。
比如我手中的sintax第一行是:
OTU_1 k:Bacteria(1.00),p:Proteobacteria(1.00),c:Betaproteobacteria(1.00),o:Burkholderiales(1.00),f:Oxalobacteraceae(1.00),g:Herbaspirillum(0.66),s:Sclerotinia_homoeocarpa(0.66) + k:Bacteria,p:Proteobacteria,c:Betaproteobacteria,o:Burkholderiales,f:Oxalobacteraceae
只要能把Bacteria、Proteobacteria、Betaproteobacteria、Burkholderiales、Oxalobacteraceae、Herbaspirillum、Sclerotinia_homoeocarpa分别提取出来即可。思路是先使用read.table读取sintax文件,得到数据框:
V1为otu名,V2相比V4信息更加全面:
> input$V4[1]
[1] "k:Bacteria,p:Proteobacteria,c:Betaproteobacteria,o:Burkholderiales,f:Oxalobacteraceae"
> input$V2[1]
[1] "k:Bacteria(1.00),p:Proteobacteria(1.00),c:Betaproteobacteria(1.00),o:Burkholderiales(1.00),f:Oxalobacteraceae(1.00),g:Herbaspirillum(0.66),s:Sclerotinia_homoeocarpa(0.66)"
所以选择V2来进行字符串分割。首先将逗号作为分隔符进行分割,得到:
[1] "k:Bacteria(1.00)" "p:Proteobacteria(1.00)" "c:Betaproteobacteria(1.00)"
[4] "o:Burkholderiales(1.00)" "f:Oxalobacteraceae(1.00)" "g:Herbaspirillum(0.66)"
[7] "s:Sclerotinia_homoeocarpa(0.66)"
再对每个字符串,例如"k:Bacteria(1.00)" ,使用substr进行截取第3(包括)到倒数第6(不包括)个字符:
substr(str,3,nchar(str)-6)
[1] "Bacteria"
nchar()是求字符串长度的函数。
最后将每个otu裁剪的字符串放在一个列表中,作为一行的taxonomy,最后将每列进行rbind拼接并输出即可。
注意,因为有些otu没用s和g等,所以如果分割出来的字符串数量少于7个,将使用'na'进行补充。
最后代码如下:
filepath='D:/Rdata/50/result' #your file path,please change before use
setwd(filepath)
# transform otus.sintax to taxnomy.txt
input = read.table("otus.sintax", header=F, sep="\t", comment.char="", stringsAsFactors=F)#read the sintax
title=paste('OTUID','Kingdom','Phylum','Class','Order','Family','Genus','Species',sep='\t')#creatw the header
output=title
for(i in 1:nrow(input)){
name=input$V1[i]#otu name
seed=input$V2[i]
seed_group=strsplit(seed,',')#string group
g=1#answer group for every row
for (j in 1:length(seed_group[[1]])) {
str=seed_group[[1]][j]
g[j]=substr(str,3,nchar(str)-6) #intercept string
}
for (j in length(g)+1:7){
g[j]='na' #fill with "na"
}
data=paste(name,g[1],g[2],g[3],g[4],g[5],g[6],g[7],sep='\t')
output=rbind(output,data)#splicing each row
}
write(output,'taxonomy.txt')
因为懒,没有写成函数调用形式。