Использование HMM для профильного анализа белкового семейства
Этот пример показывает, как профили HMM используются для характеристики белковых семейств. Профильный анализ является ключевым инструментом в биоинформатике. Общие методы парного сравнения обычно не чувствительны и достаточно специфичны для анализа отдаленно связанных последовательностей. Напротив, профили Hidden Markov Model (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 функция с одним выходным аргументом помещает выровненные последовательности в массив char.
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');
Выравнивание тандемно повторяющихся областей.
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');
Изучение профилей HMM
Функция showhmmprof является интерактивным инструментом для исследования профиля HMM. Попробуйте правый и левый клики мыши над рисунками модели. Для каждой модели существует три графика: (1) вероятности излучения символа в Match состояния, (2) вероятности излучения символа в Insert состояния и (3) Transition вероятности.
showhmmprof(hmm_7tm,'scale','logodds')
Альтернативным методом исследования HMM профиля является создание логотипа последовательности из нескольких выравниваний. Логотип последовательности отображает частоту основ, обнаруженных в каждой позиции в заданной области, обычно для сайта привязки. Используя hmm_7tm последовательности, примите во внимание фрагмент пептидного рецептора, связанного с паратиреоидными гормонами (предшественник), обнаруженную на n-конце PTRR_Human последовательности. The 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)