Используя HMMs для анализа профиля семейства белков
В этом примере показано, как профили HMM используются, чтобы охарактеризовать семейства белков. Анализ профиля является ключевым инструментом в биоинформатике. Общие попарные методы сравнения обычно не чувствительны и достаточно специфичны для анализа отдаленно связанных последовательностей. В отличие от этого профили Скрытой модели Маркова (HMM) обеспечивают лучшую альтернативу, чтобы связать последовательность запроса со статистическим описанием семейства последовательностей. Профили HMM используют специфичную для положения систему выигрыша, чтобы получить информацию о степени сохранения в различных положениях в нескольких выравнивание этих последовательностей. Анализ профиля HMM может использоваться для нескольких, упорядочивают выравнивание, для поиска в базе данных, чтобы анализировать сегментацию состава и шаблона последовательности, и предсказать структуру белка и определить местоположение генов путем предсказания открытых рамок считывания.
Доступ базы данных PFAM
Начните этот пример с уже созданного HMM семейства белков. Получите модель для известного 7-кратного трансмембранного приемника от базы данных Sanger Institute. Номер ключа PFAM является PF00002. Также получите предварительно выровненные последовательности, используемые, чтобы обучить эту модель. Больше информации о базе данных PFAM может быть найдено в http://pfam.xfam.org/.
hmm_7tm = gethmmprof(2); seed_seqs = gethmmalignment(2,'type','seed');
Для вашего удобства ранее загруженные последовательности включены в MAT-файл. Обратите внимание на то, что данные в общедоступных репозиториях часто курируются и обновляются; поэтому результаты этого примера могут немного отличаться, когда вы используете актуальные наборы данных.
load('gpcrfam.mat','hmm_7tm','seed_seqs')
Модели и выравнивания могут также храниться и проанализированный в позже непосредственно из файлов с помощью pfamhmmread, fastaread и multialignread функции.
Отобразите имена и содержимое первых трех загруженных последовательностей с помощью seqdisp команда.
seqdisp(seed_seqs([1 2 3]),'row',70)
ans =
23x81 char array
'>VIPR2_HUMAN/123-371 '
' 1 YILVKAIYTL GYSVS.LMSL ATGSIILCLF .RKLHCTR.N YIHLNLFLSF ILRAISVLVK .DDVLYSSS.'
' 71 GTLHCPD... .......... .......... ....QPSSW. ..V.GCKLSL VFLQYCIMAN FFWLLVEGLY'
'141 LHTLLVA... ...MLPP.RR CFLAYLLIGW GLPTVCIGAW TAAR...... .........L YLED......'
'211 ......TGC. WDTN.DHSVP W....WVIRI PILISIIVNF VLFISIIRIL LQKLT..... .SPDVGGNDQ'
'281 SQY....... .......... .......... ....KRLAKS TLLLIPLFGV HYMV..FAVF PISI...S.S'
'351 KYQILFELCL GSF....QGL VV '
' '
'>VIPR_CARAU/100-348 '
' 1 FRSVKIGYTI GHSVS.LISL TTAIVILCMS .RKLHCTR.N YIHMHLFVSF ILKAIAVFVK .DAVLYDVIQ'
' 71 ESDNCS.... .......... .......... .....TASV. ....GCKAVI VFFQYCIMAS FFWLLVEGLY'
'141 LHALLAVS.. ...FFSE.RK YFWWYILIGW GGPTIFIMAW SFAK...... .........A YFND......'
'211 ......VGC. WDIIENSDLF W....WIIKT PILASILMNF ILFICIIRIL RQKIN..... .CPDIGRNES'
'281 NQY....... .......... .......... ....SRLAKS TLLLIPLFGI NFII..FAFI PENI...K.T'
'351 ELRLVFDLIL GSF....QGF VV '
' '
'>VIPR1_RAT/140-386 '
' 1 YNTVKTGYTI GYSLS.LASL LVAMAILSLF .RKLHCTR.N YIHMHLFMSF ILRATAVFIK .DMALFNSG.'
' 71 EIDHCS.... .......... .......... .....EASV. ....GCKAAV VFFQYCVMAN FFWLLVEGLY'
'141 LYTLLAVS.. ...FFSE.RK YFWGYILIGW GVPSVFITIW TVVR...... .........I YFED......'
'211 ......FGC. WDTI.INSSL W....WIIKA PILLSILVNF VLFICIIRIL VQKLR..... .PPDIGKNDS'
'281 SPY....... .......... .......... ....SRLAKS TLLLIPLFGI HYVM..FAFF PDNF...K.A'
'351 QVKMVFELVV GSF....QGF VV '
Больше информации относительно того, как хранить информацию HMM профиля в структуре MATLAB®, найдено в справке для hmmprofstruct.
Профилируйте выравнивание HMM
Чтобы протестировать инструмент выравнивания HMM профиля, можно перестроить последовательности от нескольких выравнивание к модели HMM. Сначала сотритесь, периоды в последовательностях раньше форматировали загруженные выровненные последовательности. Выполнение этого удаляет информацию о выравнивании из последовательностей.
seqs = strrep({seed_seqs.Sequence},'.','');
names = {seed_seqs.Header};
Теперь выровняйте все белки к профилю HMM.
fprintf('Aligning sequences ') scores = zeros(numel(seqs),1); aligned_seqs = cell(numel(seqs),1); for sn=1:numel(seqs) fprintf('.') [scores(sn),aligned_seqs{sn}]=hmmprofalign(hmm_7tm,seqs{sn}); end fprintf('\n')
Aligning sequences ................................
Затем отправьте результаты в веб-браузер, чтобы лучше исследовать новое несколько выравнивание. Столбцы, отмеченные * в нижней части, указывают, когда модель была в "соответствии", или "удалите" состояние.
hmmprofmerge(aligned_seqs,names,scores)
Можно также исследовать выравнивание из командного окна; hmmprofmerge функция с одним выходным аргументом помещает выровненные последовательности в массив символов.
str = hmmprofmerge(aligned_seqs); str(1:10,1:80)
ans =
10x80 char array
'YILVKAIYTLGYSVS.LMSLATGSIILCLF.RKLHCTR.NYIHLNLFLSFILRAISVLVK.DDVLYSSSG-TLH......'
'FRSVKIGYTIGHSVS.LISLTTAIVILCMS.RKLHCTR.NYIHMHLFVSFILKAIAVFVK.DAVLYDVIQESDN......'
'YNTVKTGYTIGYSLS.LASLLVAMAILSLF.RKLHCTR.NYIHMHLFMSFILRATAVFIK.DMALFNSG-EIDH......'
'FGAIKTGYTIGHSLS.LISLTAAMIILCIF.RKLHCTR.NYIHMHLFMSFIMRAIAVFIK.DIVLFESG-ESDH......'
'YLSVKALYTVGYSTS.LVTLTTAMVILCRF.RKLHCTR.NFIHMNLFVSFMLRAISVFIK.DWILYAEQD-SSH......'
'FSTVKIIYTTGHSIS.IVALCVAIAILVAL.RRLHCPR.NYIHTQLFATFILKASAVFLK.DAAIFQGDS-TDH......'
'LSTLKQLYTAGYATS.LISLITAVIIFTCF.RKFHCTR.NYIHINLFVSFILRATAVFIK.DAVLFSDET-QNH......'
'FDRLGMIYTVGYSVS.LASLTVAVLILAYF.RRLHCTR.NYIHMHLFLSFMLRAVSIFVK.DAVLYSGATLDEA......'
'FERLYVMYTVGYSIS.FGSLAVAILIIGYF.RRLHCTR.NYIHMHLFVSFMLRATSIFVK.DRVVHAHIGVKEL......'
'ALNLFYLTIIGHGLS.IASLLISLGIFFYF.KSLSCQR.ITLHKNLFFSFVCNSVVTIIH.LTAVANNQALVAT......'
Поиск подобия со сравнением последовательности
При наличии профиля HHM, который описывает это семейство, имеет несколько преимуществ перед простым сравнением последовательности. Предположим, что у вас есть новый олигонуклеотид, что вы хотите относиться к 7-трансмембранному семейству приемников. В данном примере получите последовательность белка от NCBI и извлеките последовательность аминокислоты.
mousegpcr = getgenpept('NP_783573');
Bai3 = mousegpcr.Sequence;
Эта последовательность также обеспечивается в MAT-файле gpcrfam.mat.
load('gpcrfam.mat','mousegpcr') Bai3 = mousegpcr.Sequence; seqdisp(Bai3,'row',70)
ans =
22x82 char array
' 1 MKAVRNLLIY IFSTYLLVMF GFNAAQDFWC STLVKGVIYG SYSVSEMFPK NFTNCTWTLE NPDPTKYSIY'
' 71 LKFSKKDLSC SNFSLLAYQF DHFSHEKIKD LLRKNHSIMQ LCSSKNAFVF LQYDKNFIQI RRVFPTDFPG'
' 141 LQKKVEEDQK SFFEFLVLNK VSPSQFGCHV LCTWLESCLK SENGRTESCG IMYTKCTCPQ HLGEWGIDDQ'
' 211 SLVLLNNVVL PLNEQTEGCL TQELQTTQVC NLTREAKRPP KEEFGMMGDH TIKSQRPRSV HEKRVPQEQA'
' 281 DAAKFMAQTG ESGVEEWSQW SACSVTCGQG SQVRTRTCVS PYGTHCSGPL RESRVCNNTA LCPVHGVWEE'
' 351 WSPWSLCSFT CGRGQRTRTR SCTPPQYGGR PCEGPETHHK PCNIALCPVD GQWQEWSSWS HCSVTCSNGT'
' 421 QQRSRQCTAA AHGGSECRGP WAESRECYNP ECTANGQWNQ WGHWSGCSKS CDGGWERRMR TCQGAAVTGQ'
' 491 QCEGTGEEVR RCSEQRCPAP YEICPEDYLI SMVWKRTPAG DLAFNQCPLN ATGTTSRRCS LSLHGVASWE'
' 561 QPSFARCISN EYRHLQHSIK EHLAKGQRML AGDGMSQVTK TLLDLTQRKN FYAGDLLVSV EILRNVTDTF'
' 631 KRASYIPASD GVQNFFQIVS NLLDEENKEK WEDAQQIYPG SIELMQVIED FIHIVGMGMM DFQNSYLMTG'
' 701 NVVASIQKLP AASVLTDINF PMKGRKGMVD WARNSEDRVV IPKSIFTPVS SKELDESSVF VLGAVLYKNL'
' 771 DLILPTLRNY TVVNSKVIVV TIRPEPKTTD SFLEIELAHL ANGTLNPYCV LWDDSKSNES LGTWSTQGCK'
' 841 TVLTDASHTK CLCDRLSTFA ILAQQPREIV MESSGTPSVT LIVGSGLSCL ALITLAVVYA ALWRYIRSER'
' 911 SIILINFCLS IISSNILILV GQTQTHNKSI CTTTTAFLHF FFLASFCWVL TEAWQSYMAV TGKIRTRLIR'
' 981 KRFLCLGWGL PALVVATSVG FTRTKGYGTD HYCWLSLEGG LLYAFVGPAA AVVLVNMVIG ILVFNKLVSR'
'1051 DGILDKKLKH RAGQMSEPHS GLTLKCAKCG VVSTTALSAT TASNAMASLW SSCVVLPLLA LTWMSAVLAM'
'1121 TDKRSILFQI LFAVFDSLQG FVIVMVHCIL RREVQDAFRC RLRNCQDPIN ADSSSSFPNG HAQIMTDFEK'
'1191 DVDIACRSVL HKDIGPCRAA TITGTLSRIS LNDDEEEKGT NPEGLSYSTL PGNVISKVII QQPTGLHMPM'
'1261 SMNELSNPCL KKENTELRRT VYLCTDDNLR GADMDIVHPQ ERMMESDYIV MPRSSVSTQP SMKEESKMNI'
'1331 GMETLPHERL LHYKVNPEFN MNPPVMDQFN MNLDQHLAPQ EHMQNLPFEP RTAVKNFMAS ELDDNVGLSR'
'1401 SETGSTISMS SLERRKSRYS DLDFEKVMHT RKRHMELFQE LNQKFQTLDR FRDIPNTSSM ENPAPNKNPW'
'1471 DTFKPPSEYQ HYTTINVLDT EAKDTLELRP AEWEKCLNLP LDVQEGDFQT EV '
Во-первых, использование локального выравнивания сравнивает новую последовательность с одной из последовательностей в нескольких выравнивание. Например, используйте первую последовательность, в этом случае человеческий белок 'VIPR2'. Алгоритм Смита-лодочника (swalign) может использовать выигрыш матриц. Выигрыш матриц может получить вероятность замены символов. Последовательности в этом примере, как известно, только отдаленно связаны, таким образом, BLOSUM30 является хорошим выбором для матрицы выигрыша.
VIPR2 = seqs{1};
[sc_aa_affine, alignment] = swalign(Bai3,VIPR2,'ScoringMatrix',...
'blosum30','gapopen',5,'extendgap',3,'showscore',true);
sc_aa_affine
sc_aa_affine = 69.6000
Путем рассмотрения пробела выигрыша, по-видимому, связаны обе последовательности. Однако это отношение не могло быть выведено из точечной диаграммы.
Bai3_aligned_region = strrep(alignment(1,:),'-',''); seqdotplot(VIPR2,Bai3_aligned_region,7,2) ylabel('VIPR2'); xlabel('Bai3');
Имеет любой эти два примера достаточно доказательства, чтобы подтвердить, что эти последовательности связаны? Один способ протестировать это состоит в том, чтобы случайным образом создать поддельную последовательность с тем же распределением аминокислот и видеть, как это выравнивается семейству. Заметьте, что счет локального выравнивания между поддельной последовательностью и белком VIPR2 не значительно ниже, чем счет выравнивания между Bia3 и белками VIPR2. Чтобы гарантировать воспроизводимость результатов этого примера, мы сбрасываем глобальный случайный генератор.
rng(0,'twister'); fakeSeq = randseq(1000,'FROMSTRUCTURE',aacount(VIPR2)); sc_fk_affine = swalign(fakeSeq,VIPR2,'ScoringMatrix','blosum30',... 'gapopen',5,'extendgap',3,'showscore',true)
sc_fk_affine = 60.4000
В отличие от этого, когда вы выравниваете обе последовательности семейству, использующему обученный HMM профиля, счет выравнивания целевой последовательности к профилю семейства значительно больше, чем счет выравнивания поддельной последовательности.
sc_aa_hmm = hmmprofalign(hmm_7tm,Bai3) sc_fk_hmm = hmmprofalign(hmm_7tm,fakeSeq)
sc_aa_hmm = 214.5286 sc_fk_hmm = -49.1624
Исследование опций выравнивания HMM профиля
Так же к swalign функция выравнивания, когда вы используете выравнивания профиля, можно визуализировать пробел выигрыша с помощью showscore опция к hmmprofalign функция.
Отобразите Bai3, выровненный к 7tm_2 семейство.
hmmprofalign(hmm_7tm,Bai3,'showscore',true); title('log-odds score for best path: Bai3');
Отобразите "поддельную" последовательность, выровненную к 7tm_2 семейство.
hmmprofalign(hmm_7tm,fakeSeq,'showscore',true); title('log-odds score for best path: fake sequence');
Отобразите Bai3, глобально выровненный к 7tm_2 семейство.
[sc_aa_hmm,align,ptrs] = hmmprofalign(hmm_7tm,Bai3); Bai3_hmmaligned_region = Bai3(min(ptrs):max(ptrs)); hmmprofalign(hmm_7tm,Bai3_hmmaligned_region,'showscore',true); title('log-odds score for best path: Bai3 aligned globally');
Выровняйтесь tandemly повторил области.
naa = numel(Bai3_hmmaligned_region); repeats = randseq(1000,'FROMSTRUCTURE',aacount(Bai3)); %artificial example repeats(200+(1:naa)) = Bai3_hmmaligned_region; repeats(500+(1:naa)) = Bai3_hmmaligned_region; repeats(700+(1:naa)) = Bai3_hmmaligned_region; hmmprofalign(hmm_7tm,repeats,'showscore',true); title('log-odds score for best path: Bai3 tandem repeats');
Поиск областей фрагмента
В MATLAB® можно искать области фрагмента путем ручной активации B->M и M->E вероятности перехода модели HMM.
hmm_7tm_f = hmm_7tm; hmm_7tm_f.BeginX(3:end)=.002; hmm_7tm_f.MatchX(1:end-1,4)=.002;
Создайте случайную последовательность или модель фрагмента, с небольшой вставкой белка Bai3:
fragment = randseq(1000,'FROMSTRUCTURE',aacount(Bai3));
fragment(501:550) = Bai3_hmmaligned_region(101:150);
Попытайтесь выровнять случайную последовательность со вставленным пептидом к обеим моделям, глобальной переменной и модели фрагмента:
hmmprofalign(hmm_7tm,fragment,'showscore',true); title('log-odds score for best path: PF00002 global '); hmmprofalign(hmm_7tm_f,fragment,'showscore',true); title('log-odds score for best path: PF00002 fragment domains');
Исследование профиля HMMs
Функциональный showhmmprof интерактивный инструмент должен исследовать HMM профиля. Попробуйте правые и левые клики мыши по фигурам модели. Существует три графика для каждой модели: (1) вероятности эмиссии символа в Match состояния, (2) вероятности эмиссии символа в Insert состояния, и (3) Transition вероятности.
showhmmprof(hmm_7tm,'scale','logodds')
Альтернативный метод, чтобы исследовать HMM профиля путем создания логотипа последовательности из нескольких выравнивание. Логотип последовательности отображает частоту основ, найденных в каждом положении в данной области, обычно для связывающего участка. Используя hmm_7tm последовательности, считайте фрагмент Паращитовидной железы связанным с гормоном приемником пептида (предшественник) найденный в n-конечной-остановке последовательности PTRR_Human. seqlogo позволяет быстрое визуальное сравнение того, как хорошо эта область сохраняется через 7tm семейство.
seqlogo(str,'startat',1,'endat',20,'alphabet','AA')
Профилируйте оценку
HMMs профиля может также быть оценен от выравнивания кратного. Когда новые последовательности, связанные с семейством, найдены, возможно повторно оценить параметры модели.
hmm_7tm_new = hmmprofestimate(hmm_7tm,str)
hmm_7tm_new =
struct with fields:
Name: '7tm_2'
PfamAccessionNumber: 'PF00002.19'
ModelDescription: '7 transmembrane receptor (Secretin family)'
ModelLength: 243
Alphabet: 'AA'
MatchEmission: [243x20 double]
InsertEmission: [243x20 double]
NullEmission: [0.0768 0.0418 0.0396 0.0305 0.0201 0.0378 ... ]
BeginX: [244x1 double]
MatchX: [242x4 double]
InsertX: [242x2 double]
DeleteX: [242x2 double]
FlankingInsertX: [2x2 double]
LoopX: [2x2 double]
NullX: [2x1 double]
В случае, если ваши последовательности не предварительно выравниваются, можно также использовать multialign функция прежде, чем оценить новый профиль HMM. Возможно совершенствовать профиль HMM путем перестройки последовательностей к модели и переоценки модели итеративно, пока вы не сходитесь к локально оптимальной модели.
aligned_seqs = multialign(seqs); hmm_7tm_ma = hmmprofestimate(hmmprofstruct(270),aligned_seqs) showhmmprof(hmm_7tm_ma,'scale','logodds') close; close; % close insertion emission prob. and transition prob.
hmm_7tm_ma =
struct with fields:
ModelLength: 270
Alphabet: 'AA'
MatchEmission: [270x20 double]
InsertEmission: [270x20 double]
NullEmission: [0.0500 0.0500 0.0500 0.0500 0.0500 0.0500 0.0500 ... ]
BeginX: [271x1 double]
MatchX: [269x4 double]
InsertX: [269x2 double]
DeleteX: [269x2 double]
FlankingInsertX: [2x2 double]
LoopX: [2x2 double]
NullX: [2x1 double]
Выровняйте все последовательности к новой модели.
fprintf('Aligning sequences ') scores = zeros(numel(seqs),1); aligned_seqs = cell(numel(seqs),1); for sn=1:numel(seqs) fprintf('.') [scores(sn),aligned_seqs{sn}]=hmmprofalign(hmm_7tm_ma,seqs{sn}); end fprintf('\n') str = hmmprofmerge(aligned_seqs); str(1:10,1:80)
Aligning sequences ................................
ans =
10x80 char array
'YILVKAIYTLGYSVSLMSLATGSIILCLF.RKLHCTRNYIHLNLFLSFILRAISVLVKDDVLYSS---SGTLHCP-....'
'FRSVKIGYTIGHSVSLISLTTAIVILCMS.RKLHCTRNYIHMHLFVSFILKAIAVFVKDAVLYDVIQ--ESDNCS-....'
'YNTVKTGYTIGYSLSLASLLVAMAILSLF.RKLHCTRNYIHMHLFMSFILRATAVFIKDMALFNS---GEIDHCS-....'
'FGAIKTGYTIGHSLSLISLTAAMIILCIF.RKLHCTRNYIHMHLFMSFIMRAIAVFIKDIVLFES---GESDHCH-....'
'YLSVKALYTVGYSTSLVTLTTAMVILCRF.RKLHCTRNFIHMNLFVSFMLRAISVFIKDWILYAE---QDSSHCF-....'
'FSTVKIIYTTGHSISIVALCVAIAILVAL.RRLHCPRNYIHTQLFATFILKASAVFLKDAAIFQG---DSTDHCS-....'
'LSTLKQLYTAGYATSLISLITAVIIFTCF.RKFHCTRNYIHINLFVSFILRATAVFIKDAVLFSD---ETQNHCL-....'
'FDRLGMIYTVGYSVSLASLTVAVLILAYF.RRLHCTRNYIHMHLFLSFMLRAVSIFVKDAVLYSGATLDEAERLTE....'
'FERLYVMYTVGYSISFGSLAVAILIIGYF.RRLHCTRNYIHMHLFVSFMLRATSIFVKDRVVHAHIGVKELESLIM....'
'ALNLFYLTIIGHGLSIASLLISLGIFFYF.KSLSCQRITLHKNLFFSFVCNSVVTIIHLTAVANNQALVATNP---....'
Покажите выровненные последовательности в Браузере документации.
hmmprofmerge(aligned_seqs,names,scores)