В этом примере показано, как извлечь некоторые последовательности из 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=3×2 struct
Start
Stop
Теперь посмотрите на ORFs в мыши ген HEXA. В этом случае ORF находится также на первой системе координат.
mouseORFs = seqshoworfs(mouseHEXA.Sequence)
mouseORFs=3×2 struct
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 - 1
humanStop = 599
mouseStop = find(mouseProtein(mouseStart:end)=='*',1) + mouseStart - 1
mouseStop = 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'