多个(非生物、离散状态)序列的比对
Alignment of multiple (non-biological, discrete state) sequences
我有一些数据描述了一组有序的离散事件(或状态)。有 34 种可能的状态,它们可以以任何顺序出现并且可以重复。每个事件序列可以包含任意数量的事件,关键是事件序列要多于 2 个。我的最终目标是将这些序列聚类成相似的子集,但我的直觉是这没有意义,除非这些序列对齐,使得等效事件在所有序列中占据相同的位置。
我非常熟悉生物序列的多重比对,但我遇到的所有软件(MUSCLE、MAFFT、T-COFFEE、Clustal* 等)都需要 DNA、RNA 或 AA 序列,而且我的状态比任何一个都多,所以我无法让它们工作。
我发现了 pairwise 对齐算法的各种实现,例如 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 中实现的东西。
假设我们需要匹配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")
我建议使用 MAFFT sequence alignment。通常,这用于对齐生物序列,但它可以选择使用 --anysymbol 对齐文本。请注意,MAFFT 是一个 bash 脚本,需要一个 input/output 文件。
输入文件(mafft_anysymbol_input.txt):
>Seq1
ABCDEFG
>Seq2
ACDGH
>Seq3
BDEFEGI
>Seq4
AH
>Seq5
DEGHI
R 代码到 运行 bash 脚本:
#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 and --textmatrix。它需要 ascii 代码输入(额外的数据类型转换),但您可以选择通过分数关联相似的字母(无论您选择如何定义相似)。例如,您可以关联大小写字母或字母 with/without 重音符号。
我有一些数据描述了一组有序的离散事件(或状态)。有 34 种可能的状态,它们可以以任何顺序出现并且可以重复。每个事件序列可以包含任意数量的事件,关键是事件序列要多于 2 个。我的最终目标是将这些序列聚类成相似的子集,但我的直觉是这没有意义,除非这些序列对齐,使得等效事件在所有序列中占据相同的位置。
我非常熟悉生物序列的多重比对,但我遇到的所有软件(MUSCLE、MAFFT、T-COFFEE、Clustal* 等)都需要 DNA、RNA 或 AA 序列,而且我的状态比任何一个都多,所以我无法让它们工作。
我发现了 pairwise 对齐算法的各种实现,例如 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 中实现的东西。
假设我们需要匹配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")
我建议使用 MAFFT sequence alignment。通常,这用于对齐生物序列,但它可以选择使用 --anysymbol 对齐文本。请注意,MAFFT 是一个 bash 脚本,需要一个 input/output 文件。
输入文件(mafft_anysymbol_input.txt):
>Seq1
ABCDEFG
>Seq2
ACDGH
>Seq3
BDEFEGI
>Seq4
AH
>Seq5
DEGHI
R 代码到 运行 bash 脚本:
#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 and --textmatrix。它需要 ascii 代码输入(额外的数据类型转换),但您可以选择通过分数关联相似的字母(无论您选择如何定义相似)。例如,您可以关联大小写字母或字母 with/without 重音符号。