如何获得~1000种蛋白质的成对"序列相似性得分"？

Question

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

我想获得每对蛋白质的成对序列相似性得分.

R中的任何包装都可用于获得蛋白质序列的爆炸相似性评分？

Answer 1

根据蔡斯的建议,bioconductor确实是要走的路,特别是Biostrings包裹.要安装后者,我建议安装核心bioconductor库:

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

这样您就可以涵盖所有依赖项.现在,对准2个蛋白质序列或为此事,你将需要使用任何两个序列pairwiseAlignment的Biostrings.给定protseq.fasta2个序列的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我建议你阅读本文.

另一种方法是进行多序列比对而不是成对比对.您可以使用bio3d并从那里使用seqaln函数来对齐fasta文件中的所有序列.