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

Question

我正在使用 R package msa, a core Bioconductor 包，用于多序列比对。在 msa 内，我使用 MUSCLE 比对算法来比对蛋白质序列。

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

test.fa文件是一个标准的fasta，如下所示（t运行cated，为简洁起见）：

>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, ...) 保存它。

Answer 1

我认为您应该首先按照显示特定 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, ...)
}

Answer 2

我写了一个函数：

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')

Answer 3

您可以使用 export.fasta function from bio2mds 库。

# 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)

Answer 4

您可以先将您的msa对齐（“AAStringSet”）先转换为“对齐”对象，然后导出为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")

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

How to output results of 'msa' package in R to fasta

r

bioinformatics

alignment

fasta