2

我有一些数据描述了一组有序的离散事件（或状态）。有 34 种可能的状态，它们可能以任何顺序出现并可能重复。每个事件序列可以包含任意数量的事件，关键是有超过 2 个事件序列。我的最终目标是将这些序列聚集成相似的子集，但我的直觉是，除非这些序列对齐使得等效事件在所有序列中占据相同位置，否则这将毫无意义。

我非常熟悉生物序列的多重比对，但我遇到的所有软件（MUSCLE、MAFFT、T-COFFEE、Clustal* 等）都需要 DNA、RNA 或 AA 序列，而且我有更多的状态比任何这些，所以我不能让他们工作。

我在 R 中找到了成对比对算法的各种实现，例如 Needleman-Wunsch，但到目前为止还没有遇到任何多序列比对算法的任何通用（非生物）实现。

例如，假设我的数据如下所示：

1: ABCDEFG
2: ACDGH
3: BDEFEGI
4: AH
5: DEGHI

我的目标是让它看起来像这样：

1: ABCDEF-G--
2: A-CD---GH-
3: -B-DEFE--I
4: A-------H-
5: ---DE--GHI

其中-符号表示此序列中没有事件。这是一个简化的示例，实际上，我正在寻找以-与生物序列 MSA 算法相同的方式惩罚打开间隙 ( ) 的方法。

我发现的唯一一个似乎可能做到这一点的软件是 Alphamalig ( http://alggen.lsi.upc.es/recerca/align/alphamalig/intro-alphamalig.html ) 但它很旧而且我不能让它在我的机器上工作。理想情况下，我想要可以在 R 中实现的东西。

4

2 回答 2

2

假设我们需要匹配LETTERS，一个选项是str_match，然后将 更改NA为-，paste

library(stringr)
library(dplyr)
f1 <- Vectorize(function(x) str_match(x, LETTERS))
out1 <- f1(v1)
do.call(paste0, as.data.frame(t(replace_na(out1[!!rowSums(!is.na(out1)),], '-'))))
#[1] "ABCDEFG--" "A-CD--GH-" "-B-DEFG-I" "A------H-" "---DE-GHI"

match拆分后也可以

lst <- strsplit(v1, "")
mx <- match(max(sapply(lst, tail, 1)), LETTERS)
sapply(lst, function(x) paste(replace_na(x[match(LETTERS[seq_len(mx)], 
           x)], '-'), collapse=""))

数据

v1 <- c("ABCDEFG", "ACDGH", "BDEFEGI", "AH", "DEGHI")
于 2019-04-20T17:48:18.017 回答
2

我建议使用MAFFT 序列比对。通常，这用于对齐生物序列，但它可以选择使用 --anysymbol 对齐文本。请注意，MAFFT 是一个 bash 脚本，需要一个输入/输出文件。

输入文件（mafft_anysymbol_input.txt）：

>Seq1
ABCDEFG
>Seq2
ACDGH
>Seq3
BDEFEGI
>Seq4
AH
>Seq5
DEGHI

运行 bash 脚本的 R 代码：

#Be sure that input/output and R files share the same path, otherwise you'll have to specify the path in the mafft script call.
x <- 'mafft --anysymbol mafft_anysymbol_input.txt > mafft_anysymbol_output.txt'
system(x)

输出文件的内容（mafft_anysymbol_output.txt）：

>Seq1
ABCDEFG--
>Seq2
-ACDGH---
>Seq3
--BDEFEGI
>Seq4
----AH---
>Seq5
---DEGHI-

编辑 - 我现在看到您熟悉生物比对工具。如果你想为你的文本对齐创建一个自定义的评分矩阵，请查看 mafft 选项--text 和 --textmatrix。它需要 ascii 代码输入（额外的数据类型转换），但您可以选择按分数关联相似的字母（但是您选择定义相似的）。例如，您可以关联大写和小写字母，或带/不带重音符号的字母。

于 2019-04-24T15:45:58.070 回答