Convert/transform R 中的丰度 (OTU) table/data.frame（到 fasta 文件）

Question

我目前正在处理一个大型数据集，到目前为止，我可以通过无数次 google 搜索和长时间的尝试和错误会话很好地解决我的所有 ideas/problems。我已经成功地使用 plyr 和 reshape 函数对我的不同数据集进行了一些转换并学到了很多东西，但我认为我已经达到了我目前的 R 知识不再对我有帮助的地步。

即使我的问题听起来很具体（即 OTU table 和 fasta 文件），我想我的尝试是跨许多不同领域（而不仅仅是生物信息学）的通用 R 应用程序。

现在，我合并了一个参考序列文件，丰度table，我想根据这个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 命令在 mothur wont work, because it needs a fasta file as input (which I obviously don't have). It could be very possible that there is something on Bioconductor（OTUbase？），但老实说，我不知道从哪里开始，我很高兴能得到任何帮助。我真的很想在 R 中这样做，因为我喜欢使用它，但也非常欢迎任何其他想法。

//小编辑：

下面的两个答案都非常有效（请参阅我的评论）- 我还发现了两个可能的不太优雅和非 R 的解决方法（尚未测试）：

因为我已经有一个分类文件和大量的 OTU table，我认为 mothur 命令 make.biom could be used to create a biom-format file。我还没有使用过 biom 文件，但我认为有一些工具和脚本可用于将 biom 文件数据再次保存为 fasta
convert Qiime files to oligotyping format - 这还需要一个分类文件和一个 Otu table

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

Answer 1

试试这个，它逐行遍历数据帧并连接序列的重复：

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,])获取序列的命名向量。

Convert/transform R 中的丰度 (OTU) table/data.frame（到 fasta 文件）

Convert/transform an abundance (OTU) table/data.frame (to a fasta file) in R

r

fasta

bioconductor

mothur