我们有2个DNA序列(字符串):
>1
ATGCAT
135198
>2
ATCAT
预期输出:首先,我们需要对齐这两个字符串,然后按索引获取相关注释:
ATGCAT
AT-CAT
13-198
第一部分可以使用Biostrings包:
library(Biostrings)
p <- DNAString("ATCAT")
s <- DNAString("ATGCAT")
s_annot <- "135198"
x <- pairwiseAlignment(pattern = p, subject = s)
aligned(x)
# A DNAStringSet instance of length 1
# width seq
# [1] 6 AT-CAT
as.character(x)
# [1] "AT-CAT"
as.matrix(x)
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] "A" "T" "-" "C" "A" "T"
第二部分的当前解决方法:
annot <- unlist(strsplit(s_annot, ""))
annot[ which(c(as.matrix(x)) == "-") ] <- "-"
# [1] "1" "3" "-" "1" "9" "8"
效果很好,但是我想知道是否有Biostrings方式(或任何其他包),也许通过将注释保留在metadata插槽中,然后在对齐后得到匹配元数据中匹配碱基的注释,如下所示:
getSlots("DNAString")
# shared offset length elementMetadata metadata
# "SharedRaw" "integer" "integer" "DataTable_OR_NULL" "list"
# just an idea, non-working code
s@metadata <- unlist(strsplit(s_annot , ""))
x <- pairwiseAlignment(pattern = p, subject = s)
metadata(x)
# [[1]]
# [1] "1" "3" "-" "1" "9" "8"
注意:
我们有2个DNA序列(字符串):> 1 ATGCAT 135198> 2 ATCAT预期输出:首先,我们需要对齐这2个字符串,然后按索引获取相关注释:ATGCAT AT-CAT 13-198第一部分...
可能的解决方案:
根据要求,是一个Biopython解决方案: