如何在R中进行基本的多序列比对？

Question

（我曾尝试在BioStars上提问，但为了让文本挖掘的人认为有更好的解决方案，我也在这里重新发布）

我要完成的任务是对齐多个序列。

我没有一个基本的匹配模式。我只知道“真实”模式应该是长度为“30”，并且我所拥有的序列在随机点引入了缺失值。

以下是这些序列的示例。左侧显示缺失值的实际位置，右侧显示我们能够观察到的序列。

我的目标是仅使用右列中的序列（基于每个位置上的许多字母相同的事实）来重构左列。

                     Real_sequence           The_sequence_we_see
1   CGCAATACTAAC-AGCTGACTTACGCACCG CGCAATACTAACAGCTGACTTACGCACCG
2   CGCAATACTAGC-AGGTGACTTCC-CT-CG   CGCAATACTAGCAGGTGACTTCCCTCG
3   CGCAATGATCAC--GGTGGCTCCCGGTGCG  CGCAATGATCACGGTGGCTCCCGGTGCG
4   CGCAATACTAACCA-CTAACT--CGCTGCG   CGCAATACTAACCACTAACTCGCTGCG
5   CGCACGGGTAAGAACGTGA-TTACGCTCAG CGCACGGGTAAGAACGTGATTACGCTCAG
6   CGCTATACTAACAA-GTG-CTTAGGC-CTG   CGCTATACTAACAAGTGCTTAGGCCTG
7   CCCA-C-CTAA-ACGGTGACTTACGCTCCG   CCCACCTAAACGGTGACTTACGCTCCG

这里有一个示例代码，可以复制上面的示例:

ATCG <- c("A","T","C","G")
set.seed(40)
original.seq <- sample(ATCG, 30, T)
seqS <- matrix(original.seq,200,30, T)
change.letters <- function(x, number.of.changes = 15, letters.to.change.with = ATCG) 
{
    number.of.changes <- sample(seq_len(number.of.changes), 1)
    new.letters <- sample(letters.to.change.with , number.of.changes, T)
    where.to.change.the.letters <- sample(seq_along(x) , number.of.changes, F)
    x[where.to.change.the.letters] <- new.letters
    return(x)
}
change.letters(original.seq)
insert.missing.values <- function(x) change.letters(x, 3, "-") 
insert.missing.values(original.seq)

seqS2 <- t(apply(seqS, 1, change.letters))
seqS3 <- t(apply(seqS2, 1, insert.missing.values))

seqS4 <- apply(seqS3,1, function(x) {paste(x, collapse = "")})
require(stringr)
# library(help=stringr)
all.seqS <- str_replace(seqS4,"-" , "")

# how do we allign this?
data.frame(Real_sequence = seqS4, The_sequence_we_see = all.seqS)

我理解，如果我只有一个字符串和一个模式，我将能够使用

library(Biostrings)
pairwiseAlignment(...)

但在我提出的情况下，我们需要将许多序列彼此对齐（而不是将它们对齐到一个模式）。

是否有一种已知的方法可以在R中完成这个任务？

- Tal Galili

Answer 1

虽然这是一个相当老的帖子，但我不想错过提到自Bioconductor 3.1以来，有一个名为'msa'的软件包，它实现了三种不同的多序列比对算法的接口：ClustalW、ClustalOmega和MUSCLE。该软件包可在所有主要平台（Linux/Unix、Mac OS和Windows）上运行，并且是自包含的，意味着您不需要安装任何外部软件。更多信息可以在http://www.bioinf.jku.at/software/msa/和http://www.bioconductor.org/packages/release/bioc/html/msa.html找到。

- UBod

---

感谢您更新信息。我丢失了被接受的答案，但我并不介意。毕竟，那个答案已经拖延了7年。 - Joris Meys

Answer 2

在R中编写对齐算法对我来说似乎不是个好主意，但MUSCLE算法在bio3d包中有一个R接口（函数seqaln()）。请注意，您必须先安装此算法。

或者，您可以使用任何可用的算法（例如ClustalW, MAFFT, T-COFFEE），并使用bioconductor功能将多个序列对齐导入到R中。例如，请参见此处。

- Joris Meys

---

嗨，Joris。很高兴知道这样的接口存在（我在搜索中找不到它）。我不确定你所说的“你必须先安装此算法”是什么意思（你是指软件包吗？）-但我想我很快就会弄清楚:)谢谢你的帮助:) Tal - Tal Galili

这不是一个软件包，而是一个独立的可执行文件：http://www.drive5.com/muscle/downloads.htm - Ben Bolker

正如Ben所说，MUSCLE是一个广为人知且表现良好的多序列比对算法。因此，您需要先安装可执行文件，然后尝试使用bio3d包。这可能有点奇怪，因为bio3d主要用于蛋白质建模，但MUSCLE也可以比对多个DNA序列。 - Joris Meys

请注意，现在有新的软件包可用，包括 msa，它可以处理 MUSCLE、CLUSTALW 和 CLUSTALOMEGA 算法。另外，它不依赖于任何外部软件。 - Roman Luštrik

Answer 3

您可以使用DECIPHER包在R中执行多重对齐操作。根据您的示例，它应该是这样的：

library(DECIPHER)
dna <- DNAStringSet(all.seqS)
aligned_DNA <- AlignSeqs(dna)

它非常快速，并且至少与此处列出的其他方法一样准确（详见论文）。希望这能帮到您！

- Erik Wright

---

DECIPHER是否假定读入的数据是FASTA格式？我没有看到任何关于字符串为原始字符的示例。 - Vass

Answer 4

您正在寻找多个序列的全局比对算法。

在提问之前，您是否查看了维基百科？首先了解什么是全局比对，然后再寻找多序列比对。

维基百科并没有提供很多关于算法的细节，但这篇论文更好一些。

- Jules Olléon

---

@Jules 这并没有告诉他如何在 R 中实现。我认为 Tal 并不会因为收到一些他可以轻松通过谷歌搜索得到的链接而受益。而且很可能在他提出这个问题之前就已经做过了。 - Joris Meys

你好，朱尔斯。感谢你的良好意图和分享演示文稿的链接。祝好，Tal - Tal Galili

@Joris 鉴于他提出问题的方式，看起来 Tal 不知道这是一个标准的生物信息学问题，所以我认为一些关于这个领域的背景知识可能会对他有所帮助...但你说得对，我应该把它作为评论而不是答案发布。 - Jules Olléon

@Jules：没有恶意。你不知道Tal是一位生物统计学家，对R有相当深入的了解，所以我在评论中可能有点严厉。 - Joris Meys

嗨，朱尔斯。你说得对，这是我第一次尝试使用R来完成这样的任务，因此我才会提出问题。正如你所说，建议阅读背景资料非常友善。乔里斯 - 你和我都知道很少有人拥有深入的R知识，而像我今天这样的状态，我不敢自认为其中之一。我仍然非常感激你的好意 :) 你的，塔尔 - Tal Galili

@Jules：鉴于Tal的友好评论，我取消了我的负评。为了能够这样做，我不得不在你的帖子中添加了一个空格。 - Joris Meys