Этот пример показывает, как профили HMM используются для характеристики семейств белков. Профильный анализ является ключевым инструментом в биоинформатике. Общие методы попарного сравнения обычно не чувствительны и достаточно специфичны для анализа отдаленно связанных последовательностей. В противоположность этому профили скрытой модели Маркова (HMM) обеспечивают лучшую альтернативу для связи последовательности запроса со статистическим описанием семейства последовательностей. Профили HMM используют зависящую от позиции систему оценки для сбора информации о степени сохранения в различных положениях во множественном выравнивании этих последовательностей. Анализ профиля HMM может быть использован для выравнивания нескольких последовательностей, для поиска в базе данных, для анализа состава последовательности и сегментации шаблона, а также для прогнозирования структуры белка и определения местоположения генов путем прогнозирования открытых кадров считывания.
Начните этот пример с уже построенного HMM семейства белков. Извлеките модель для хорошо известного 7-кратного трансмембранного рецептора из базы данных Института Сэнгера. Номер ключа 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. Сначала удалите периоды в последовательностях, используемых для форматирования загруженных выровненных последовательностей. При этом информация о выравнивании удаляется из последовательностей.
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 ................................
Затем отправьте результаты в Web Browser, чтобы лучше изучить новое несколько трасс. Столбцы, помеченные * внизу, указывают, когда модель находилась в состоянии «совпадение» или «удаление».
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
Аналогично 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');
Выравнивание повторяемых совместно доменов.
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');
Функция 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')
Параметры HMM профиля также могут быть оценены по нескольким трассам. По мере обнаружения новых последовательностей, связанных с семейством, можно повторно оценить параметры модели.
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: [1x20 double]
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: [1x20 double]
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)