Согласно предложению Чейза, bioconductor - это действительно путь, в частности пакет Biostrings.Для установки последнего я бы предложил установить базовую библиотеку bioconductor следующим образом:
source("http://bioconductor.org/biocLite.R")
biocLite()
Таким образом, вы охватите все зависимости.Теперь, чтобы выровнять 2 последовательности белка или любые две последовательности по этому вопросу, вам нужно будет использовать pairwiseAlignment из Biostrings.Имеется файл fasta protseq.fasta из 2 последовательностей, который выглядит следующим образом:
>protein1
MYRALRLLARSRPLVRAPAAALASAPGLGGAAVPSFWPPNAAR
MASQNSFRIEYDTFGELKVPNDKYYGAQTVRSTMNFKIGGVTE
RMPTPVIKAFGILKRAAAEVNQDYGLDPKIANAIMKAADEVAE
GKLNDHFPLVVWQTGSGTQTNMNVNEVISNRAIEMLGGELGSK
IPVHPNDHVNKSQ
>protein2
MRSRPAGPALLLLLLFLGAAESVRRAQPPRRYTPDWPSLDSRP
LPAWFDEAKFGVFIHWGVFSVPAWGSEWFWWHWQGEGRPYQRF
MRDNYPPGFSYADFGPQFTARFFHPEEWADLFQAAGAKYVVLT
TKHHEGFTNW*
Если вы хотите глобально выровнять эти 2 последовательности, используя, скажем, 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.