r - 如何获得约 1000 个蛋白质的成对 "sequence similarity score"？

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 作为替换矩阵来全局比对这 2 个序列，则打开一个空位的惩罚为 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 文件中的所有序列。