从 R 中的局部比对中识别氨基酸取代

Question

我想确定我的序列中特定感兴趣区域的位置和氨基酸变化，并将该信息存储在表格中。

是否可以通过使用 bioconductor 包在 R 中做这样的事情？我已经设法使用 DECIPHER 包中的 AlignSeqs() 对序列进行了简单的比对，但我无法自动提取序列中的差异。我从 FASTA 文件开始。

我想得到这样的结果：

Isolate ID    Reference_AA   Sample_AA   Pos
1             S              T           254
2             T              D           200
3             L              A           230

我有一个 74 AA 长的参考序列，我想查看与参考相比查询序列（比参考长得多）的差异，并在表格中列出结果。位置列与参考序列中的位置有关，与查询序列中的位置无关。我希望参考序列中的第一个 AA 从 68 开始，而不是 1。

我发现很难为此添加示例序列，因为它们往往很长，但这里有一些更短的东西可以使用（与上表无关）：

>ref
VGRALPDVR

>query1
KSSYLDYAMSVIVGTALPDVRDGLKPVHRRVLY

>query2
ELKSSYLDYAMSVIVGRAAPDVRDGLKPV

预期输出：

ID       Reference_AA    Sample_AA   Pos
query1   R               T           70
query2   A               L           72

Answer 1

由于您想要参考中的位置，您可以对参考序列使用一系列成对的对齐方式。 biostrings 包括一个 mismatchTable 函数，它将为您提供包含所需信息的数据框。使用dplyr 重新格式化：

library(Biostrings)
library(dplyr)

seqs<-readAAStringSet("test.fa")

mismatches <- function(query, ref) {
    pairwiseAlignment(ref, query, substitutionMatrix = "BLOSUM50",
                      gapOpening = 3, gapExtension = 1) %>%
      mismatchTable() %>%
      mutate(ID=names(query), 
             Pos=PatternStart+67, 
             Reference_AA=as.character(PatternSubstring),
             Sample_AA=as.character(SubjectSubstring)) %>% 
             select(ID, Reference_AA, Sample_AA, Pos)
}  

bind_rows(mismatches(seqs[2], seqs[1]), mismatches(seqs[2], seqs[1]))

#>      ID Reference_AA Sample_AA Pos
#>1 query1            R         T  70
#>2 query2            L         A  72

编辑

以下是使用 lapply 循环输入的方法：

bind_rows(lapply(seq_along(seqs[-1]), function(i) mismatches(seqs[i+1], seqs[1])))
#>   ID Reference_AA Sample_AA Pos
#>1 ref            R         T  70
#>2 ref            L         A  72