Этот пример показывает, как извлечь некоторые последовательности от GenBank®, найдите открытые рамки считывания (ORFs), и затем выровняйте последовательности с помощью глобальных и локальных алгоритмов выравнивания.
Один из многих захватывающих разделов веб-сайта NCBI является разделом Генов и болезней. Этот раздел обеспечивает всестороннее введение в медицинскую генетику.
В этом примере вы будете смотреть на гены, сопоставленные с болезнью Тея-Сакса. Болезнь Тея-Сакса является автосомальной удаляющейся болезнью, вызванной мутациями в обеих аллелях гена (HEXA, который коды для альфа-подблока hexosaminidase A) на хромосоме 15.
Ссылочная последовательность NCBI для HEXA имеет инвентарный номер NM_000520. Можно использовать функцию getgenbank, чтобы получить информацию о последовательности из репозитория данных NCBI и загрузить его в MATLAB®.
humanHEXA = getgenbank('NM_000520');
Путем выполнения BLAST ищут или путем поиска в геноме мыши, можно найти ортогональный ген, AK080777 является инвентарным номером для мыши hexosaminidase ген.
mouseHEXA = getgenbank('AK080777');
Для вашего удобства ранее загруженные последовательности включены в MAT-файл. Обратите внимание на то, что данные в общедоступных репозиториях часто курируются и обновляются; поэтому результаты этого примера могут немного отличаться, когда вы используете актуальные наборы данных.
load('hexosaminidase.mat','humanHEXA','mouseHEXA')
Можно использовать функциональный seqshoworfs, чтобы искать ORFs в последовательности для человеческого гена HEXA. Заметьте, что самый длинный ORF находится на кадре первого чтения. Выходное значение в переменной humanORFs является структурой, дающей положение запуска и кодонов остановки для всего ORFs на каждой рамке считывания.
humanORFs = seqshoworfs(humanHEXA.Sequence)
humanORFs = 1x3 struct array with fields:
Start
Stop
Теперь посмотрите на ORFs в мыши ген HEXA. В этом случае ORF находится также на первом кадре.
mouseORFs = seqshoworfs(mouseHEXA.Sequence)
mouseORFs = 1x3 struct array with fields:
Start
Stop
Первый шаг должен использовать глобальное выравнивание последовательности, чтобы искать общие черты между этими последовательностями. Вы могли посмотреть на выравнивание между последовательностями нуклеотида, но это обычно более поучительно, чтобы посмотреть на выравнивание между последовательностями белка в этом примере, мы знаем, что последовательности кодируют последовательности. Используйте функцию nt2aa, чтобы преобразовать последовательности нуклеотида в соответствующие последовательности аминокислот. Заметьте, что ген HEXA происходит в первом кадре для обеих последовательностей, в противном случае необходимо использовать входной параметр Frame, чтобы задать альтернативный кадр кодирования.
humanProtein = nt2aa(humanHEXA.Sequence); mouseProtein = nt2aa(mouseHEXA.Sequence);
Один из самых легких способов искать подобие между последовательностями с точечной диаграммой.
seqdotplot(mouseProtein,humanProtein)
Warning: Match matrix has more points than available screen pixels.<br> Scaling image by factors of 1 in X and 2 in Y.
xlabel('Human hexosaminidase A');ylabel('Mouse hexosaminidase A');
С настройками по умолчанию точечную диаграмму немного трудно интерпретировать, таким образом, можно попробовать немного более строгую точечную диаграмму.
seqdotplot(mouseProtein,humanProtein,4,3)
Warning: Match matrix has more points than available screen pixels.<br> Scaling image by factors of 1 in X and 2 in Y.
xlabel('Human hexosaminidase A');ylabel('Mouse hexosaminidase A');
Диагональная строка указывает, что существует, вероятно, хорошее выравнивание, таким образом, можно теперь смотреть на глобальное выравнивание с помощью функционального nwalign, который использует алгоритм Needleman-Wunsch.
[score, globalAlignment] = nwalign(humanProtein,mouseProtein);
Функциональный showalignment отображает выравнивание в окне фигуры MATLAB с соответствием и подобными остатками, подсвеченными в различных цветах.
showalignment(globalAlignment);
Выравнивание очень хорошо за исключением терминальных сегментов. Например, заметьте разреженные совпадающие пары в первых положениях. Это происходит, потому что глобальное выравнивание пытается обеспечить соответствие полностью к концам и существует точка, где штраф за открытие новых разрывов сопоставим со счетом соответствия с остатками. В некоторых случаях желательно удалить штраф разрыва, добавленный в концах глобального выравнивания; это позволяет вам лучше совпадать с этой парой последовательностей. Этот метод обычно известен как 'полуглобальное' выравнивание или 'glocal' выравнивание.
[score, globalAlignment] = nwalign(humanProtein,mouseProtein,'glocal',true);
showalignment(globalAlignment);
Другой способ совершенствовать ваше выравнивание только при помощи последовательностей белка. Заметьте, что выровненная область разграничена запуском (M-метионин), и остановитесь (*) аминокислоты в последовательностях. Если последовательность сокращена так, чтобы только переведенные области были рассмотрены, то кажется вероятным, что вы получите лучшее выравнивание. Используйте команду find, чтобы искать индекс аминокислоты запуска в каждой последовательности:
humanStart = find(humanProtein == 'M',1)humanStart = 70
mouseStart = find(mouseProtein == 'M',1)mouseStart = 11
Точно так же используйте команду find, чтобы искать индекс первой остановки, происходящей после запуска перевода. Необходимо соблюдать особую осторожность, потому что существует также остановка в самом начале последовательности humanProtein.
humanStop = find(humanProtein(humanStart:end)=='*',1) + humanStart - 1humanStop = 599
mouseStop = find(mouseProtein(mouseStart:end)=='*',1) + mouseStart - 1mouseStop = 539
Используйте эти индексы, чтобы обрезать последовательности.
humanSeq = humanProtein(humanStart:humanStop); humanSeqFormatted = seqdisp(humanSeq)
humanSeqFormatted = 9x70 char array
' 1 MTSSRLWFSL LLAAAFAGRA TALWPWPQNF QTSDQRYVLY PNNFQFQYDV SSAAQPGCSV'
' 61 LDEAFQRYRD LLFGSGSWPR PYLTGKRHTL EKNVLVVSVV TPGCNQLPTL ESVENYTLTI'
'121 NDDQCLLLSE TVWGALRGLE TFSQLVWKSA EGTFFINKTE IEDFPRFPHR GLLLDTSRHY'
'181 LPLSSILDTL DVMAYNKLNV FHWHLVDDPS FPYESFTFPE LMRKGSYNPV THIYTAQDVK'
'241 EVIEYARLRG IRVLAEFDTP GHTLSWGPGI PGLLTPCYSG SEPSGTFGPV NPSLNNTYEF'
'301 MSTFFLEVSS VFPDFYLHLG GDEVDFTCWK SNPEIQDFMR KKGFGEDFKQ LESFYIQTLL'
'361 DIVSSYGKGY VVWQEVFDNK VKIQPDTIIQ VWREDIPVNY MKELELVTKA GFRALLSAPW'
'421 YLNRISYGPD WKDFYIVEPL AFEGTPEQKA LVIGGEACMW GEYVDNTNLV PRLWPRAGAV'
'481 AERLWSNKLT SDLTFAYERL SHFRCELLRR GVQAQPLNVG FCEQEFEQT* '
mouseSeq = mouseProtein(mouseStart:mouseStop); mouseSeqFormatted = seqdisp(mouseSeq)
mouseSeqFormatted = 9x70 char array
' 1 MAGCRLWVSL LLAAALACLA TALWPWPQYI QTYHRRYTLY PNNFQFRYHV SSAAQAGCVV'
' 61 LDEAFRRYRN LLFGSGSWPR PSFSNKQQTL GKNILVVSVV TAECNEFPNL ESVENYTLTI'
'121 NDDQCLLASE TVWGALRGLE TFSQLVWKSA EGTFFINKTK IKDFPRFPHR GVLLDTSRHY'
'181 LPLSSILDTL DVMAYNKFNV FHWHLVDDSS FPYESFTFPE LTRKGSFNPV THIYTAQDVK'
'241 EVIEYARLRG IRVLAEFDTP GHTLSWGPGA PGLLTPCYSG SHLSGTFGPV NPSLNSTYDF'
'301 MSTLFLEISS VFPDFYLHLG GDEVDFTCWK SNPNIQAFMK KKGFTDFKQL ESFYIQTLLD'
'361 IVSDYDKGYV VWQEVFDNKV KVRPDTIIQV WREEMPVEYM LEMQDITRAG FRALLSAPWY'
'421 LNRVKYGPDW KDMYKVEPLA FHGTPEQKAL VIGGEACMWG EYVDSTNLVP RLWPRAGAVA'
'481 ERLWSSNLTT NIDFAFKRLS HFRCELVRRG IQAQPISVGC CEQEFEQT* '
Если вы выровняете эти две последовательности и затем просмотрите их, то вы будете видеть очень хорошее глобальное выравнивание.
[score, alignment] = nwalign(humanSeq,mouseSeq); showalignment(alignment);
Открытая информация о рамке считывания также доступна от вывода команды seqshoworfs, но индексы основаны на последовательностях нуклеотида. Используйте эти индексы, чтобы обрезать исходные последовательности нуклеотида и затем перевести их в аминокислоты.
humanPORF = nt2aa(humanHEXA.Sequence(humanORFs(1).Start(1):humanORFs(1).Stop(1))); mousePORF = nt2aa(mouseHEXA.Sequence(mouseORFs(1).Start(1):mouseORFs(1).Stop(1))); [score, ORFAlignment] = nwalign(humanPORF,mousePORF); showalignment(ORFAlignment);
Также можно использовать информацию об области кодирования (CDS) от структуры данных GenBank, чтобы найти область кодирования генов.
idx = humanHEXA.CDS.indices; humanCodingRegion = humanHEXA.Sequence(idx(1):idx(2)); idx = mouseHEXA.CDS.indices; mouseCodingRegion = mouseHEXA.Sequence(idx(1):idx(2));
Можно также получить перевод областей кодирования от этой структуры.
humanTranslatedRegion = humanHEXA.CDS.translation; mouseTranslatedRegion = mouseHEXA.CDS.translation;
Вместо того, чтобы обрезать последовательности, чтобы искать лучшее выравнивание, альтернативный подход должен использовать локальное выравнивание. Функциональный swalign выполняет локальное выравнивание с помощью алгоритма Смита-лодочника. Это показывает очень хорошее выравнивание для целой области кодирования и разумное подобие для нескольких остатков вне в обоих концы гена.
[score, localAlignment] = swalign(humanProtein,mouseProtein); showalignment(localAlignment);
Все функции выравнивания последовательности, обеспеченные в MATLAB, могут быть настроены. Например, путем изменения строк и столбцов матрицы выигрыша можно выровнять последовательности дополнением а не идентичностью. В этом случае можно переупорядочить матрицу выигрыша NUC44; положительный счет дан для дополнений, в то время как отрицательный счет дан в противном случае. Первые 30 нуклеотидов от мыши ген HEXA будут выровнены к его дополнению.
[M, info] = nuc44; map = nt2int(seqcomplement(info.Order))
map = 1x15 uint8 row vector
4 3 2 1 6 5 8 7 9 10 14 13 12 11 15
Mc = M(:,map)
Mc = 15×15
-4 -4 -4 5 -4 1 1 -4 -4 1 -1 -1 -1 -4 -2
-4 -4 5 -4 1 -4 1 -4 1 -4 -1 -1 -4 -1 -2
-4 5 -4 -4 -4 1 -4 1 1 -4 -1 -4 -1 -1 -2
5 -4 -4 -4 1 -4 -4 1 -4 1 -4 -1 -1 -1 -2
-4 1 -4 1 -4 -1 -2 -2 -2 -2 -1 -3 -1 -3 -1
1 -4 1 -4 -1 -4 -2 -2 -2 -2 -3 -1 -3 -1 -1
1 1 -4 -4 -2 -2 -4 -1 -2 -2 -3 -3 -1 -1 -1
-4 -4 1 1 -2 -2 -1 -4 -2 -2 -1 -1 -3 -3 -1
-4 1 1 -4 -2 -2 -2 -2 -1 -4 -1 -3 -3 -1 -1
1 -4 -4 1 -2 -2 -2 -2 -4 -1 -3 -1 -1 -3 -1
⋮
[score, compAlignment] = nwalign(mouseHEXA.Sequence(1:30), ... seqcomplement(mouseHEXA.Sequence(1:30)), 'SCORINGMATRIX', ... Mc, 'ALPHABET', 'NT')
score = 150
compAlignment = 3x30 char array
'GCTGCTGGAAGGGGAGCTGGCCGGTGGGCC'
'::::::::::::::::::::::::::::::'
'CGACGACCTTCCCCTCGACCGGCCACCCGG'