如何将R中'msa'包的结果输出到fasta答案

【问题标题】：How to output results of 'msa' package in R to fasta如何将R中'msa'包的结果输出到fasta
【发布时间】：2020-07-02 22:23:02
【问题描述】：

我正在使用R package msa，一个核心Bioconductor 包，用于多序列比对。在msa 中，我使用 MUSCLE 比对算法来比对蛋白质序列。

library(msa)
myalign <- msa("test.fa", method=c("Muscle"), type="protein",verbose=FALSE)

test.fa 文件是一个标准的 fasta 文件如下（为了简洁，被截断）：

>sp|P31749|AKT1_HUMAN_RAC
MSDVAIVKEGWLHKRGEYIKTWRPRYFLL
>sp|P31799|AKT1_HUMAN_RAC
MSVVAIVKEGWLHKRGEYIKTWRFLL

当我在文件上运行代码时，我得到：

 MUSCLE 3.8.31   

Call:
   msa("test.fa", method = c("Muscle"), type = "protein", verbose = FALSE)

MsaAAMultipleAlignment with 2 rows and 480 columns
    aln 
[1] MSDVAIVKEGWLHKRGEYIKTWRPRYFLL
[2] MSVVAIVKEGWLHKRGEYIKTWR---FLL
Con MS?VAIVKEGWLHKRGEYIKTWR???FLL

如您所见，非常合理的对齐方式。

我想将有间隙的比对写入 fasta 文件，最好不使用共有序列（例如 Con 行）。所以，我想要：

>sp|P31749|AKT1_HUMAN_RAC
MSDVAIVKEGWLHKRGEYIKTWRPRYFLL
>sp|P31799|AKT1_HUMAN_RAC
MSVVAIVKEGWLHKRGEYIKTWR---FLL

我查看了msa 帮助，并且该软件包似乎没有内置方法可以写入任何文件类型，fasta 或其他文件。

seqinr package 看起来很有希望，因为它可能会将此输出读取为 msf 格式，尽管它很奇怪。但是，seqinr 似乎需要一个文件作为起点读入。我什至无法使用write(myalign, ...) 保存它。

【问题讨论】：

你真的应该停止使用大写的包名。

标签： r alignment bioinformatics fasta

【解决方案1】：

我写了一个函数：

alignment2Fasta <- function(alignment, filename) {
  sink(filename)

  n <- length(rownames(alignment))
  for(i in seq(1, n)) {
    cat(paste0('>', rownames(alignment)[i]))
    cat('\n')
    the.sequence <- toString(unmasked(alignment)[[i]])
    cat(the.sequence)
    cat('\n')  
  }

  sink(NULL)
}

用法：

mySeqs <- readAAStringSet('test.fa')
myAlignment <- msa(mySeqs)
alignment2Fasta(myAlignment, 'out.fasta')

【讨论】：

【解决方案2】：

我认为您应该首先遵循帮助页面中显示具有特定 read 函数的输入的示例，然后使用对齐方式：

mySeqs <- readAAStringSet("test.fa")
myAlignment <- msa(mySeqs)

然后rownames 函数将传递序列名称：

 rownames(myAlignment)
[1] "sp|P31749|AKT1_HUMAN_RAC" "sp|P31799|AKT1_HUMAN_RAC"

（不是您要求的，但将来可能有用。）然后如果您执行：

detail(myAlignment)   #function actually in Biostrings

....您会在交互模式下获得一个可以保存的文本文件

2 29
sp|P31749|AKT1_HUMAN_RAC   MSDVAIVKEG WLHKRGEYIK TWRPRYFLL
sp|P31799|AKT1_HUMAN_RAC   MSVVAIVKEG WLHKRGEYIK TWR---FLL

如果你想尝试破解一个函数，你可以得到一个用代码编写的文件，然后查看正在使用的 Biostrings detail 函数代码

> showMethods( f= 'detail')
Function: detail (package Biostrings)
x="ANY"
x="MsaAAMultipleAlignment"
    (inherited from: x="MultipleAlignment")
x="MultipleAlignment"

showMethods( f= 'detail', classes='MultipleAlignment', includeDefs=TRUE)
Function: detail (package Biostrings)
x="MultipleAlignment"
function (x, ...) 
{
    .local <- function (x, invertColMask = FALSE, hideMaskedCols = TRUE) 
    {
        FH <- tempfile(pattern = "tmpFile", tmpdir = tempdir())
        .write.MultAlign(x, FH, invertColMask = invertColMask, 
            showRowNames = TRUE, hideMaskedCols = hideMaskedCols)
        file.show(FH)
    }
    .local(x, ...)
}

【讨论】：

我永远感激不尽！第一个选项将完全符合我的要求。随着这项工作的推进，我还将努力更好地理解第二个选择。最后，我将来会更加注意包名称的细节，因为我对这些名称是由包作者特别选择的这一事实很敏感。
根据您对 Biostrings 和内置函数的帮助，我发现还有一个内置函数可以保存为 phylip 格式。不出所料，它是 write.phylip(x, ...)。去搞清楚！再次感谢！
我在阅读文档时看到了该功能，但它似乎并不完全符合您的要求。在本页底部：rdrr.io/bioc/Biostrings/api，您可以从 Biostrings 中找到更完整的 write 函数列表。

【解决方案3】：

您可以使用bio2mds 库中的export.fasta 函数。

# reading of the multiple sequence alignment of human GPCRS in FASTA format: aln <- import.fasta(system.file("msa/human_gpcr.fa", package = "bios2mds")) export.fasta(aln)

【讨论】：

【解决方案4】：

您可以先将您的 msa 对齐（“AAStringSet”）转换为“align”对象，然后导出为 fasta，如下所示：

library(msa)

library(bios2mds)

mysequences <-readAAStringSet("test.fa")

alignCW  <- msa(mysequences)

#https://rdrr.io/bioc/msa/man/msaConvert.html

alignCW_as_align <- msaConvert(alignCW, "bios2mds::align")

export.fasta(alignCW_as_align, outfile = "test_alignment.fa", ncol = 60, open = "w")

【讨论】：