我有大量 fasta 格式的蛋白质序列。
我想获得每对蛋白质的成对序列相似性分数。
R 中的任何包都可用于获得蛋白质序列的爆炸相似性分数？
最佳答案
根据蔡斯的建议， bioconductor 确实是要走的路，特别是 Biostrings 包。要安装后者，我建议像这样安装核心bioconductor库:
source("http://bioconductor.org/biocLite.R")
biocLite()
这样，您将涵盖所有依赖项。现在，要对齐 2 个蛋白质序列或任何两个序列，您需要使用
pairwiseAlignment 中的 Biostrings 。给定一个由 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 文件中的所有序列。
关于r - 如何获得约 1000 种蛋白质的成对 "sequence similarity score"？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/6529675/