在R中转换/转换丰度(OTU)表/data.framework(转换为fasta文件)

Question

目前我正在处理一个大型数据集，到目前为止，我可以通过无数的google搜索和长期的尝试与错误会话来解决我的所有想法/问题。我已经成功地使用plyr和重塑函数对我的不同数据集进行了一些转换，并且学到了很多，但我认为我已经到了一个我现在的R知识将不再对我有帮助的地步。

即使我的问题听起来非常具体(即OTU表和fasta文件)，我想我的尝试是在许多不同的领域(而不仅仅是生物信息学)的一个常见的R应用程序。

现在，我已经将一个引用序列文件与一个丰度表合并，我想根据这个data.frame的信息生成一个特定的文件--一个fasta文件。

我的df现在看起来有点像这样：

repSeq     sw.1.102 sw.3.1021 sw.30.101 sw.5.1042 ...
ACCT-AGGA  3        0         1         0
ACCT-AGGG  1        1         2         0
ACTT-AGGG  0        1         0         25
...

生成的文件应该如下所示：

>sw.1.102_1
ACCT-AGGA
>sw.1.102_2
ACCT-AGGA
>sw.1.102_3
ACCT-AGGA
>sw.1.102_4
ACCT-AGGG
>sw.3.1021_1
ACCT-AGGG
>sw.3.1021_2
ACTT-AGGG
>sw.30.101_1
ACCT-AGGA
>sw.30.101_2
ACCT-AGGG
...

如您所见，我希望使用有关每个示例(即sw.n)的(引用)序列数的信息来创建一个(fasta)文件。

我对R中的循环没有经验(我只在简单的处理尝试中使用基本循环)，但我认为这在这里可以做到。我从SeqinR包中找到了write.fasta函数，但在那里找不到任何解决方案。deunique.seqs命令在蛀虫中无法工作，因为它需要一个fasta文件作为输入(显然我没有)。在生物导体 (OTUbase?)上可能有一些东西，但老实说，我不知道从哪里开始，我很高兴得到任何帮助。我真的很想在R里做这个，因为我喜欢和它一起工作，但是任何其他的想法也是非常受欢迎的。

//小编辑：

以下两个答案都非常有效(见我的评论)--我还发现了两种可能的不那么优雅和非R的解决办法(尚未测试)：

• 因为我已经有了一个分类法文件和一个丰富的OTU表，所以我认为可以使用mothur命令make.biom来创建一个biom格式文件。我还没有使用biom文件，但我认为有一些工具和脚本可以再次将biom文件数据保存为fasta。
• 将奇姆文件转换成少类型格式 -这还需要一个分类法文件和一个Otu表。

不确定这两种方法是否有效--因此，如果我错了，请纠正我。

Answer 1

尝试一下，它逐行遍历dataframe，并将重复的序列连在一起：

fasta_seq<-apply(df,1,function(x){
        p<-x[1]
        paste(unlist(mapply(function(x,y,z){
                if(as.numeric(y)>0) {paste(">",x,"_",(z+1):(z+y),"\n",p,"\n",sep="")}
        },colnames(df)[-1],as.numeric(x[-1]),c(0,lag(cumsum(as.numeric(x[-1])))[-1]),USE.NAMES=F)),collapse="")                
        })

write(paste(fasta_seq,collapse=""),"your_file.txt")

Answer 2

这里是您的数据，强制使用矩阵(对于同构类型的矩形数据来说，这是一种更自然的表示)。

df <- read.delim(textConnection(
    "repSeq     sw.1.102 sw.3.1021 sw.30.101 sw.5.1042
     ACCT-AGGA  3        0         1         0
     ACCT-AGGG  1        1         2         0
     ACTT-AGGG  0        1         0         25"
    ), sep="", row.names=1)
m <- as.matrix(df)

棘手的部分是弄清楚如何对重复的列名条目进行编号。我通过创建适当长度的序列和取消列表来做到这一点。然后，我创建了一个包含两行的矩阵，第一行(根据原始矩阵中条目的要求复制colnames() )是id，第二行是序列。

csum <- colSums(m)
idx <- unlist(lapply(csum, seq_len), use.names=FALSE)
res <- matrix(c(sprintf(">%s_%d", rep(colnames(m), csum), idx), # id
                rep(rownames(m)[row(m)], m)),                   # sequence
              nrow=2, byrow=TRUE)

使用writeLines(res, "your.fasta")写出结果，或使用setNames(res[2,], res[1,])获取序列的命名向量。