Используя 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: [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)

Документация