如何获得约1000种蛋白质的成对“序列相似性分数”？

Question

我有大量fasta格式的蛋白质序列。

我想要得到每对蛋白质的成对序列相似性得分。

R中的任何包都可以用来获得蛋白质序列的blast相似性分数？

Answer 1

根据Chase的建议，bioconductor确实是可行的方法，尤其是Biostrings包。要安装后者，我建议按如下方式安装核心bioconductor库：

source("http://bioconductor.org/biocLite.R")
biocLite()

这样你就可以覆盖所有的依赖项了。现在，要比对2个蛋白质序列或任何两个序列，您将需要使用Biostrings的pairwiseAlignment。给定一个包含2个序列的fasta protseq.fasta文件，如下所示：

>protein1
MYRALRLLARSRPLVRAPAAALASAPGLGGAAVPSFWPPNAAR
MASQNSFRIEYDTFGELKVPNDKYYGAQTVRSTMNFKIGGVTE
RMPTPVIKAFGILKRAAAEVNQDYGLDPKIANAIMKAADEVAE
GKLNDHFPLVVWQTGSGTQTNMNVNEVISNRAIEMLGGELGSK
IPVHPNDHVNKSQ
>protein2
MRSRPAGPALLLLLLFLGAAESVRRAQPPRRYTPDWPSLDSRP
LPAWFDEAKFGVFIHWGVFSVPAWGSEWFWWHWQGEGRPYQRF
MRDNYPPGFSYADFGPQFTARFFHPEEWADLFQAAGAKYVVLT
TKHHEGFTNW*

如果你想将这两个序列全局比对，假设使用BLOSUM100作为你的替换矩阵，0表示打开一个空位，-5表示扩展空位，那么：

require("Biostrings")
data(BLOSUM100)
seqs <- readFASTA("~/Desktop/protseq.fasta", strip.descs=TRUE)
alm <- pairwiseAlignment(seqs[[1]]$seq, seqs[[2]]$seq, substitutionMatrix=BLOSUM100, gapOpening=0, gapExtension=-5)

这样做的结果是(删除了一些对齐以节省空间)：

> alm
Global PairwiseAlignedFixedSubject (1 of 1)
pattern: [1] MYRALRLLARSRPLVRA-PAAALAS....
subject: [1] M-R-------SRP---AGPALLLLL.... 
score: -91

要仅提取每个对齐的分数，请执行以下操作：

> score(alm)
[1] -91

鉴于此，您现在可以使用一些非常简单的循环逻辑轻松地完成所有的成对对齐。为了更好地掌握使用bioconductor进行成对对齐的技巧，我建议您阅读this。

另一种方法是进行多序列比对，而不是成对比对。您可以使用bio3d和seqaln函数对fasta文件中的所有序列进行比对。

Answer 2

6年后，但是：

protr包刚刚问世，它有一个并行化的成对相似性评分函数parGOsim()。它可以获取蛋白质序列的列表，因此不需要编写循环。