usearch标准注释文件转化为标准8列注释文件 sintax2taxonomy

在易生信的某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')

因为懒,没有写成函数调用形式。

  • 0
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值