Этот пример показывает, как вычислить отношения Ka/Ks для восьми генов в H5N1 и вирусных геномах H2N3, и выполнить филогенетический анализ гена HA от вируса H5N1, изолированного от цыплят через Африку и Азию. Для филогенетического анализа вы восстановите соединяющее соседа дерево и создадите 3-D график расстояний последовательности с помощью многомерного масштабирования. Наконец, вы сопоставите географические точки, где каждая последовательность HA была найдена на региональной карте. Последовательности, используемые в этом примере, были выбраны из тематического исследования птичьего гриппа на Вычислительном Веб-сайте Геномики [1]. Примечание: итоговый раздел в этом примере требует Mapping Toolbox™.
Существует три типа вируса гриппа: Тип A, B и C. Все геномы гриппа состоят из восьми сегментов или генов что код для полимеразы B2 (PB2), полимеразы B1 (PB1), полимераза A (PA), гемагглютинин (HA), нуклеобелок (NP), нейраминидаза (NA), матрица (M1) и неструктурные белки (NS1). Примечание: у вируса Типа C есть эстераза гемагглютинина (HE), гомолог к HA.
Из трех типов гриппа Тип A имеет потенциал, чтобы быть самым разрушительным. Это влияет на птиц (его естественное водохранилище), люди и другие млекопитающие и было главной причиной глобальных эпидемий гриппа. Тип B влияет только на людей, вызывающих локальные эпидемии, и Тип C не имеет тенденцию вызывать тяжелую болезнь.
Введите, гриппы далее классифицируются в различные подтипы согласно изменениям в последовательностях аминокислот HA (H1-16) и NA (N1-9) белки. Оба белка расположены за пределами вируса. HA присоединяет вирус к клетке - хозяину, затем помогает в процессе вируса, сплавляемого в к ячейке. NA отсекает недавно созданный вирус от клетки - хозяина, таким образом, это может идти дальше к здоровой новой ячейке. Различие в составе аминокислоты в белке и recombination различного HA и белках NA способствует, чтобы Ввести способность гриппов перейти разновидности хоста (т.е. птица людям) и широкий спектр серьезности. Много новых наркотиков разрабатываются, чтобы предназначаться для HA и белков NA [2,3,4].
В 1 997, подтип H5N1 вируса птичьего гриппа, Тип вирус гриппа, сделал неожиданный скачок людям в Гонконге, вызывающем смертельные случаи шести человек. Чтобы управлять быстро распространяющейся болезнью, вся домашняя птица в Гонконге была уничтожена. Анализ последовательности вируса H5N1 показывают здесь [2,4].
Расследование отношений Ka/Ks для каждого сегмента гена вируса H5N1 обеспечит некоторое понимание, как каждый изменяется в зависимости от времени. Ka/Ks является отношением несинонимичных изменений в синонимичном в последовательности. Для более подробного объяснения отношений Ka/Ks смотрите Анализирующие Синонимичные и Несинонимичные Уровни Замены. Чтобы вычислить Ka/Ks, вам нужна копия гена от двух моментов времени. Можно использовать вирус H5N1, изолированный от цыплят в Гонконге в 1 997 и 2001. Для сравнения можно включать вирус H2N3, изолированный от крякв в Альберте в 1 977 и 1985 [1].
В целях этого примера данные о последовательности обеспечиваются в четырех структурах MATLAB®, которые были созданы genbankread.
Загрузите H5N1 и данные о последовательности H2N3.
load('birdflu.mat','chicken1997','chicken2001','mallard1977','mallard1985')
Данные в общедоступных репозиториях часто курируются и обновляются. Можно получить актуальные наборы данных при помощи функции getgenbank. Обратите внимание на то, что, если данные действительно изменились, результаты этого примера могут немного отличаться, когда вы используете актуальные наборы данных.
chicken1997 = arrayfun(@(x)getgenbank(x{:}),{chicken1997.Accession});
chicken2001 = arrayfun(@(x)getgenbank(x{:}),{chicken2001.Accession});
mallard1977 = arrayfun(@(x)getgenbank(x{:}),{mallard1977.Accession});
mallard1985 = arrayfun(@(x)getgenbank(x{:}),{mallard1985.Accession});
Можно извлечь только фрагмент кодирования последовательностей нуклеотида с помощью функции featureparse. Функция featureparse возвращает структуру с полями, содержащими информацию от раздела Features в файле GenBank включая с полем Sequence, которое содержит только последовательность кодирования.
for ii = 1:numel(chicken1997) ntSeq97{ii} = featureparse(chicken1997(ii),'feature','cds','sequence',true); ntSeq01{ii} = featureparse(chicken2001(ii),'feature','cds','sequence',true); ntSeq77{ii} = featureparse(mallard1977(ii),'feature','cds','sequence',true); ntSeq85{ii} = featureparse(mallard1985(ii),'feature','cds','sequence',true); end ntSeq97{1}
ans =
struct with fields:
Location: '<1..>2273'
Indices: [1 2273]
UnknownFeatureBoundaries: 1
gene: 'PB2'
codon_start: '1'
product: 'PB2 protein'
protein_id: 'AAF02361.1'
db_xref: 'GI:6048850'
translation: 'RIKELRDLMSQSRTREILTKTTVDHMAIIKKYTSGRQEKNPALRMKWMMAMKYPITADKRIMEMIPERNEQGQTLWSKTNDAGSDRVMVSPLAVTWWNRNGPTTSTVHYPKVYKTYFEKVERLKHGTFGPVHFRNQVKIRRRVDMNPGHADLSAKEAQDVIMEVVFPNEVGARILTSESQLTITKEKREELKNCNISPLMVAYMLERELVRKTRFLPVAGGTSSVYIEVLHLTQGTCWEQMYTPGGEVRNDDVDQSLIIAARNIVRRATVSADPLASLLEMCHSTQIGGVRMVDILKQNPTEEQAVDICKAAMGLRISSSFSFGGFTFKRTKGFSVKREEEVLTGNLQTLKIQVHEGYEEFTMVGRRATAILRKATRRMIQLIVSGRDEQSIAEAIIVAMVFSQEDCMIKAVRGDLNFVNRANQRLNPMHQLLRHFQKDAKVLFQNWGIEPIDNVMGMIGILPDMTPSTEMSLRGVRVSKMGVDEYSSTERVVVSIDRFLRVRDQQGNVLLSPEEVSETQGMEKLTITYSSSMMWEINGPESVLVNTYQWIIRNWETVKIQWSQEPTMLYNKMEFEPFQSLVPKAARSQYSGFVRTLFQQMRDVLGTFDTVQIIKLLPFAAAPPEQSRMQFSSLTVNVRGSGMRILVRGNSPAFNYNKTTKRLTILGKDAGAITEDPDEGAAGVESAVLRGFLILGKEDKRYGPALSINELSNLTKGEKANVLIGQGDVVLVMKRKRDSSILTDSQTATKRIRMAIN'
Sequence: 'agaataaaagaactaagagatttgatgtcgcaatctcgcacacgcgagatactgacaaaaaccactgtggatcatatggccataattaagaagtacacatcaggaagacaggagaagaaccccgctcttagaatgaaatggatgatggcgatgaaatacccgatcacagctgacaaaagaataatggagatgatccctgaaaggaatgagcaaggtcaaactctttggagcaaaacaaatgacgctggatcagacagggtaatggtatcacctctggctgtaacgtggtggaacagaaatggaccaacaacaagtacagtccattatccaaaggtgtataaaacctactttgaaaaggttgaaagattaaaacacggaacctttggccctgttcatttccggaatcaagtcaaaatacgccgcagggttgacatgaaccctggccatgcagatctcagcgctaaagaagcacaagatgtcatcatggaggtcgttttcccaaatgaagttggagccaggatattgacatcagagtcacagctgacaataacaaaggaaaagagggaggaactcaagaattgtaatatttctcctttaatggtggcatatatgttggaaagagaattggttcgcaagaccagattcctaccagtggctggtgggacaagcagcgtatatatagaagtattgcatttgacccaaggaacctgctgggagcagatgtacacaccaggaggggaggtaagaaatgatgatgttgaccaaagtttaatcattgctgctaggaacattgtcaggagagcaacagtatcagcagacccattggcttcactcttggagatgtgccatagcacacaaattggcggagtaagaatggtagacatccttaaacaaaacccaacagaagagcaagctgtagatatatgcaaggcagcaatgggtttaagaatcagctcatccttcagctttggagggttcactttcaaaagaacaaaggggttttctgtcaaaagagaggaagaagtgcttacaggcaacctccaaacattgaagatacaagtacatgaaggatatgaggaattcacaatggttggacgaagagcaacagccattctaagaaaagcaaccagaaggatgatccaactgatagtcagcgggagggacgagcaatcaattgctgaggcaattattgtagcaatggtgttctcacaagaagattgcatgataaaggcagtccgaggtgatttgaatttcgtaaacagagcaaatcaacgactgaaccccatgcaccaactcctgagacacttccaaaaggatgcaaaggtgctgtttcaaaactggggaattgaaccaatcgacaatgtcatggggatgattggaatattgcctgacatgacccccagcacggaaatgtcactaagaggagtgagagttagtaaaatgggggtggatgaatattctagcactgaaagagtggtcgtgagcattgaccgtttcttaagggtccgagatcagcaaggaaatgtactcctatcccctgaagaagttagtgagacacagggaatggaaaagttgacgataacttattcatcgtctatgatgtgggaaattaacggcccagaatcagttctagttaacacataccaatggatcattaggaattgggagactgtaaagatccaatggtcccaagaacccaccatgctatacaataagatggagtttgaaccatttcaatctttagtaccaaaggctgccagaagccaatatagtggatttgtgagaacgctattccagcagatgcgtgatgttttgggaacatttgacactgttcaaataatcaaactactaccatttgcagcagccccacctgaacagagtaggatgcaattttcttctctgactgtgaatgtgaggggatcaggaatgagaatacttgtgagaggtaactcccctgcgtttaactacaacaagacaactaagaggcttacaatacttgggaaggacgcaggtgcaattacagaggacccagatgaaggagcagcaggggtagagtctgcagtattgagagggtttctaattctcggcaaagaagacaaaagatatggaccagcattgagcatcaatgaactgagcaatcttacgaaaggggagaaagctaatgtattgatagggcaaggagacgtggtgttggtaatgaaacggaaacgggactctagcatacttactgacagccagacagcgaccaaaagaattcggatggccatcaatta'
Визуальный осмотр структур последовательности показал, что некоторые гены имеют варианты соединения встык, представленные в файлах GenBank. Поскольку этот анализ находится только на PB2, PB1, PA, HA, NP, НА, M1 и генах NS1, необходимо удалить любые варианты соединения встык.
Удалите варианты соединения встык из 1 997 H5N1
ntSeq97{7}(1) = [];% M2
ntSeq97{8}(1) = [];% NS2
Удалите варианты соединения встык из 1 977 H2N3
ntSeq77{2}(2) = [];% PB1-F2
ntSeq77{7}(1) = [];% M2
ntSeq77{8}(1) = [];% NS2
Удалите варианты соединения встык из 1 985 H2N3
ntSeq85{2}(2) = [];% PB1-F2
ntSeq85{7}(1) = [];% M2
ntSeq85{8}(1) = [];% NS2
Необходимо выровнять последовательности нуклеотида, чтобы вычислить отношение Ka/Ks. Выровняйте последовательности белка для каждого гена (доступный в поле 'перевода') использование функции nwalign, затем вставьте разрывы в последовательность нуклеотида с помощью seqinsertgaps. Используйте функциональный dnds, чтобы вычислить несинонимичные и синонимичные уровни замены для каждого из этих восьми генов в вирусных геномах. Если вы интересуетесь наблюдением выравниваний последовательности, установите 'многословную' опцию на истину при использовании dnds.
Названия генов гриппа
proteins = {'PB2','PB1','PA','HA','NP','NA','M1','NS1'};
Вирус H5N1
for ii = 1:numel(ntSeq97) [sc,align] = nwalign(ntSeq97{ii}.translation,ntSeq01{ii}.translation,'alpha','aa'); ch97seq = seqinsertgaps(ntSeq97{ii}.Sequence,align(1,:)); ch01seq = seqinsertgaps(ntSeq01{ii}.Sequence,align(3,:)); [dn,ds] = dnds(ch97seq,ch01seq); H5N1.(proteins{ii}) = dn/ds; end
Вирус H2N3
for ii = 1:numel(ntSeq77) [sc,align] = nwalign(ntSeq77{ii}.translation,ntSeq85{ii}.translation,'alpha','aa'); ch77seq = seqinsertgaps(ntSeq77{ii}.Sequence,align(1,:)); ch85seq = seqinsertgaps(ntSeq85{ii}.Sequence,align(3,:)); [dn,ds] = dnds(ch77seq,ch85seq); H2N3.(proteins{ii}) = dn/ds; end H5N1 H2N3
H5N1 =
struct with fields:
PB2: 0.0226
PB1: 0.0240
PA: 0.0307
HA: 0.0943
NP: 0.0517
NA: 0.1015
M1: 0.0460
NS1: 0.3010
H2N3 =
struct with fields:
PB2: 0.0048
PB1: 0.0021
PA: 0.0089
HA: 0.0395
NP: 0.0071
NA: 0.0559
M1: 0
NS1: 0.1954
Примечание: результаты отношения Ka/Ks могут отличаться от показанных на [1] из-за вариантов соединения встык последовательности.
Визуализируйте отношения Ka/Ks в 3-D гистограмме.
H5N1rates = cellfun(@(x)(H5N1.(x)),proteins);
H2N3rates = cellfun(@(x)(H2N3.(x)),proteins);
bar3([H2N3rates' H5N1rates']);
ax = gca;
ax.XTickLabel = {'H2N3','H5N1'};
ax.YTickLabel = proteins;
zlabel('Ka/Ks');
view(-115,16);
title('Ka/Ks Ratios for H5N1 (Chicken) and H2N3 (Mallard) Viruses');
NS1, HA и NA имеют больше несинонимичный с синонимичными отношениями по сравнению с другими генами и в H5N1 и в H2N3. Изменения последовательности белка в этих генах были приписаны увеличению патогенности H5N1. В частности, изменения в гене HA могут предоставить вирусу способность передать в разновидности других около птиц [2,3].
Вирус H5N1 присоединяет к ячейкам в желудочно-кишечном тракте птиц и дыхательных путях людей. Изменения в белке HA, который помогает связать вирус со здоровой ячейкой и упрощает ее присоединение к ячейке, то, что позволяет вирусу влиять на различные органы в тех же и различных разновидностях. Это может предоставить ему способность спрыгнуть с птиц людям [2,3]. Можно выполнить филогенетический анализ белка HA от вируса H5N1, изолированного от цыплят в разное время (годы) в различных областях Азии и Африки, чтобы исследовать их отношение друг другу.
Загрузите данные о последовательности аминокислот HA из 16 областей/времен из MAT-файла, обеспеченного birdflu.mat, или получите актуальные данные о последовательности от репозитория NCBI с помощью функции getgenpept.
load('birdflu.mat','HA')
HA = arrayfun(@(x)getgenpept(x{:}),{HA.Accession});
Создайте новый массив структур, содержащий поля, соответствующие последовательности аминокислот (Последовательность) и информация об источнике (Заголовок). Можно извлечь информацию об источнике от HA с помощью featureparse, затем анализируют с regexp.
for ii = 1:numel(HA) source = featureparse(HA(ii),'feature','source'); strain = regexp(source.strain,'A/[Cc]hicken/(\w+\s*\w*).*/(\d+)','tokens'); proteinHA(ii).Header = sprintf('%s_%s',char(strain{1}(1)),char(strain{1}(2))); proteinHA(ii).Sequence = HA(ii).Sequence; end proteinHA(1)
ans =
struct with fields:
Header: 'Nigeria_2006'
Sequence: 'mekivllfaivslvksdqicigyhannsteqvdtimeknvtvthaqdilekthngklcdldgvkplilrdcsvagwllgnpmcdeflnvpewsyivekinpandlcypgnfndyeelkhllsrinhfekiqiipksswsdheassgvssacpyqgrssffrnvvwlikkdnayptikrsynntnqedllvlwgihhpndaaeqtrlyqnpttyisvgtstlnqrlvpkiatrskvngqsgrmeffwtilkpndainfesngnfiapenaykivkkgdstimkseleygncntkcqtpigainssmpfhnihpltigecpkyvksnrlvlatglrnspqgerrrkkrglfgaiagfieggwqgmvdgwygyhhsneqgsgyaadkestqkaidgvtnkvnsiidkmntqfeavgrefnnlerrienlnkkmedgfldvwtynaellvlmenertldfhdsnvknlydkvrlqlrdnakelgngcfefyhrcdnecmesvrngtydypqyseearlkreeisgvklesigtyqilsiystvasslalaimvaglslwmcsngslqcrici'
Выровняйте последовательности аминокислот HA с помощью multialign и визуализируйте выравнивание с seqalignviewer.
alignHA = multialign(proteinHA); seqalignviewer(alignHA);
Вычислите расстояния между последовательностями с помощью seqpdist с методом Jukes-Cantor. Используйте seqneighjoin, чтобы восстановить филогенетическое дерево с помощью соединяющего соседа метода. Seqneighjoin возвращает объект phytree.
distHA = seqpdist(alignHA,'method','Jukes-Cantor','alpha','aa'); HA_NJtree = seqneighjoin(distHA,'equivar',alignHA);
Используйте метод view, сопоставленный с объектами phytree открыть дерево в Инструменте Phylogenetic Tree.
view(HA_NJtree);
Другой способ визуализировать отношение между последовательностями состоит в том, чтобы использовать многомерное масштабирование (MDS) с расстояниями, вычисленными для филогенетического дерева. Эта функциональность обеспечивается функцией cmdscale в Statistics and Machine Learning Toolbox™. Для получения дополнительной информации о cmdscale смотрите Классическое Многомерное Масштабирование.
[Y,eigvals] = cmdscale(distHA);
Можно использовать собственные значения, возвращенные cmdscale, чтобы помочь вести решение о том, использовать ли первые два или три измерения в графике.
sigVecs = [1:3;eigvals(1:3)';eigvals(1:3)'/max(abs(eigvals))]; report = ['Dimension Eigenvalues Normalized' ... sprintf('\n %d\t %1.4f %1.4f',sigVecs)]; display(report);
report =
'Dimension Eigenvalues Normalized
1 0.0062 1.0000
2 0.0028 0.4462
3 0.0014 0.2209'
Первые две размерности представляют значительную часть данных, но третье все еще содержит информацию, которая может помочь разрешить кластеры в данных о последовательности. Можно создать трехмерный график рассеивания с помощью функции plot3.
locations = {proteinHA(:).Header};
figure
plot3(Y(:,1),Y(:,2),Y(:,3),'ok');
text(Y(:,1)+0.002,Y(:,2),Y(:,3)+0.001,locations,'interpreter','no');
title('MDS Plot of HA Sequences');
view(-21,12);
Кластеры, кажется, соответствуют группировкам в филогенетическом дереве. Найдите последовательности, принадлежащие каждому кластеру с помощью метода subtree phytree. Одни из необходимых входных параметров subtree являются номером узла (количество листов + количество ответвлений), который будет корневым узлом нового поддерева. Для вашего примера кластер, содержащий Хэбэй и Гонконг в графике MDS, эквивалентен поддереву, корневой узел которого является Ответвлением 14, который является Узлом 30 (16 листов + 14 ответвлений).
cluster1 = get(subtree(HA_NJtree,30),'LeafNames'); cluster2 = get(subtree(HA_NJtree,21),'LeafNames'); cluster3 = get(subtree(HA_NJtree,19),'LeafNames');
Получите индекс для последовательностей, принадлежащих каждому кластеру.
[cl1,cl1_ind] = intersect(locations,cluster1);
[cl2,cl2_ind] = intersect(locations,cluster2);
[cl3,cl3_ind] = intersect(locations,cluster3);
[cl4,cl4_ind] = setdiff(locations,{cl1{:} cl2{:} cl3{:}});
Измените цветные и символы маркера на графике MDS соответствовать каждому кластеру.
h = plot3(Y(cl1_ind,1),Y(cl1_ind,2),Y(cl1_ind,3),'^',... Y(cl2_ind,1),Y(cl2_ind,2),Y(cl2_ind,3),'o',... Y(cl3_ind,1),Y(cl3_ind,2),Y(cl3_ind,3),'d',... Y(cl4_ind,1),Y(cl4_ind,2),Y(cl4_ind,3),'v'); numClusters = 4; col = autumn(numClusters); for i = 1:numClusters h(i).MarkerFaceColor = col(i,:); end set(h(:),'MarkerEdgeColor','k'); text(Y(:,1)+0.002,Y(:,2),Y(:,3),locations,'interpreter','no'); title('MDS Plot of HA Sequences'); view(-21,12);
Для более подробной информации об использовании отношений Ka/Ks phylogenetics и MDS для анализа последовательности видят Кристьанини и Хана [5].
ПРИМЕЧАНИЕ: Вам нужен Mapping Toolbox, чтобы произвести следующую фигуру.
Используя инструменты от Mapping Toolbox, можно построить местоположение, где каждый вирус был изолирован на карте Африки и Азии. Для этого вам нужны широта и долгота для каждого местоположения. Для получения информации о нахождении картографических данных в Интернете смотрите Картографические данные Доступа в Интернете для Mapping Toolbox. Широта и долгота для столицы каждой географической области, где вирусы были изолированы, обеспечиваются для этого примера.
Создайте структуру геоstruct, regionHA, который содержит географическую информацию для каждой функции или последовательность, чтобы быть отображенным. Геоstruct требуется, чтобы иметь Геометрию, Лэта и поля Лона, которые задают тип функции, широту и долготу. Эта информация используется путем отображения функций в Mapping Toolbox, чтобы отобразить картографические данные.
[regionHA(1:16).Geometry] = deal('Point'); [regionHA(:).Lat] = deal(9.10, 34.31, 15.31, 39.00, 39.00, 39.00, 55.26,... 15.56, 34.00, 33.14, 34.20, 23.00, 37.35, 44.00,... 22.11, 22.11); [regionHA(:).Lon] = deal(7.10, 69.08, 32.35, 116.00, 116.00, 116.00,... 65.18, 105.48, 114.00, 131.36, 131.40, 113.00,... 127.00, 127.00, 114.14, 114.14);
Геоstruct может также иметь поля атрибута, которые содержат дополнительную информацию о каждой функции. Добавьте поля Name и Cluster атрибута в структуру regionHA. Поле Cluster содержит кластерный номер последовательности, который вы будете использовать, чтобы идентифицировать кластерное членство последовательностей.
[regionHA(:).Name] = deal(proteinHA.Header); [regionHA(cl1_ind).Cluster] = deal(1); [regionHA(cl2_ind).Cluster] = deal(2); [regionHA(cl3_ind).Cluster] = deal(3); [regionHA(cl4_ind).Cluster] = deal(4); regionHA(1)
ans =
struct with fields:
Geometry: 'Point'
Lat: 9.1000
Lon: 7.1000
Name: 'Nigeria_2006'
Cluster: 3
Создайте структуру с помощью функции makesymbolspec, которая будет содержать маркер и спецификации цветов для каждого маркера, который будет отображен на карте. Вы передадите эту структуру функции geoshow. Маркеры символа и цвета собираются соответствовать кластерам в графике MDS.
clusterSymbols = makesymbolspec('Point',... {'Cluster',1,'Marker', '^'},... {'Cluster',2,'Marker', 'o'},... {'Cluster',3,'Marker', 'd'},... {'Cluster',4,'Marker', 'v'},... {'Cluster',[1 4],'MarkerFaceColor',autumn(4)},... {'Default','MarkerSize', 6},... {'Default','MarkerEdgeColor','k'});
Загрузите информацию об отображении и используйте функцию geoshow, чтобы построить вирусные местоположения на карте.
load coast load topo figure fig = gcf; fig.Renderer = 'zbuffer'; worldmap([-45 85],[0 160]) setm(gca,'mapprojection','robinson',... 'plabellocation',30,'mlabelparallel',-45,'mlabellocation',30) plotm(lat, long) geoshow(topo, topolegend, 'DisplayType', 'texturemap') demcmap(topo) brighten(.60) geoshow(regionHA,'SymbolSpec',clusterSymbols); title('Geographic Locations of HA Sequence in Africa and Asia')
ПРИМЕЧАНИЕ: Вам нужен Mapping Toolbox, чтобы экспортировать данные в KML-отформатированный файл.
Используя функцию kmlwrite от Mapping Toolbox, можно записать местоположение и информацию об аннотации для каждой последовательности к KML-отформатированному файлу. Google Earth отображает географические данные из файлов KML в его Наземном браузере. Функция kmlwrite Mapping Toolbox переводит геоstruct, такой как regionHA, в KML-отформатированный файл, который будет использоваться Google Earth. Для получения дополнительной информации о kmlwrite смотрите Экспортирующие Векторные Данные о Точке к KML.
Можно далее аннотировать каждую последовательность информацией от раздела Features файла GenBank с помощью функции featureparse. Можно затем использовать эту информацию, чтобы заполнить геоstruct, regionHA, и отобразить его в табличной форме как тег описания для каждого placemark в браузере Google Earth. В геоstruct обязательными полями является Геометрия, Лэт и поле Лона. Все другие поля считаются атрибутами placemark.
for i = 1:numel(HA) feats = featureparse(HA(i),'Feature','source'); regionHA(i).Strain = feats.strain; if isfield(feats,'country') regionHA(i).Country = feats.country; else regionHA(i).Country = 'N/A'; end year = regexp(regionHA(i).Name,'\d+','match'); regionHA(i).Year = year{1}; % Create a link to GenPept record through the accession number regionHA(i).AccessionNumber = ... ['<a href="http://www.ncbi.nlm.nih.gov/sites/entrez?db=Protein&cmd=search&term=',... HA(i).Accession,'">',HA(i).Accession,'</a>']; end [regionHA.SequenceLength] = deal(HA.LocusSequenceLength);
Создайте структуру атрибута с помощью функции makeattribspec, которую вы будете использовать, чтобы отформатировать таблицу описания для каждого маркера. Структура атрибута диктует порядок и форматирование каждого атрибута. Можно также использовать его, чтобы не отобразить один из атрибутов в геоstruct, regionHA.
attribStruct = makeattribspec(regionHA);
Удалите Поле имени и переупорядочьте поля в структуре атрибута.
attribStruct = rmfield(attribStruct,'Name'); attribStruct = orderfields(attribStruct,{'AccessionNumber','Strain',... 'SequenceLength','Country','Year','Cluster'}); regionHA = orderfields(regionHA,{'AccessionNumber','Strain',... 'SequenceLength','Country','Year','Cluster','Geometry','Lon','Lat',... 'Name'});
Переформатируйте метки атрибута для отображения в таблице.
attribStruct.AccessionNumber.AttributeLabel = '<b>Accession Number</b>'; attribStruct.Strain.AttributeLabel = '<b>Viral Strain</b>'; attribStruct.SequenceLength.AttributeLabel = '<b>Sequence Length</b>'; attribStruct.Country.AttributeLabel = '<b>Country of Origin</b>'; attribStruct.Year.AttributeLabel = '<b>Year Isolated</b>'; attribStruct.Cluster.AttributeLabel = '<b>Cluster Membership</b>';
Запишите геоstruct regionHA в KML-отформатированный файл во временной директории.
kmlDirectory = tempdir; filename = fullfile(kmlDirectory,'HA_geographic_locations.kml'); kmlwrite(filename,regionHA,'Description',attribStruct,'Name',{regionHA.Strain},... 'Icon','http://maps.google.com/mapfiles/kml/shapes/arrow.png','iconscale',1.5);
Можно отобразить файл KML в браузере Google Earth [6]. Google Earth должен быть установлен в системе. На платформах Windows® отобразите файл KML с:
winopen(filename)
Для пользователей UNIX и пользователей Mac, отобразите файл KML с:
cmd = 'googleearth ';
fullfilename = fullfile(pwd, filename);
system([cmd fullfilename])
В данном примере файл KML был ранее отображен с помощью Google Earth Pro. Изображение Google Earth было затем сохраненным использованием меню "File-> Save Image" Google Earth. Это - то, как данные в вашем файле KML изучают, когда загружено Google Earth. Получить это представление move вокруг и увеличить масштаб области по Азии.
Кликните по placemark, чтобы просмотреть информацию о последовательности. Инвентарный номер в каждой таблице данных является гиперссылкой на файл последовательности GenPept в Базе данных Белка NCBI.
Опционально, удалите новый файл KML из своего KML директория вывода. (Этот пример должен вымыться после себя; в действительном приложении вы, вероятно, хотели бы не использовать этот шаг.)
delete(filename)
[1] https://computationalgenomics.blogs.bristol.ac.uk/case_studies/birdflu_demo
[2] Лейвер, W.G., Bischofberger, N. и Уэбстер, R.G., "разоружая вирусы гриппа", научный американец, 280 (1):78-87, 1999.
[3] Suzuki, Y. и Masatoshi, N., "Источник и эволюция вирусных генов гемагглютинина гриппа", молекулярная биология и эволюция, 19 (4):501-9, 2002.
[4] Gambaryan, A., и др., "Эволюция приемника обязательный фенотип гриппа вирусы (H5)", Вирусология, 344 (2):432-8, 2006.
[5] Cristianini, N. и Хан, M.W., "Введение в вычислительную геномику: подход тематических исследований", издательство Кембриджского университета, 2007.
[6] Изображения Google Earth были получены с помощью Google Earth Pro. Для получения дополнительной информации о Google Earth и Google Earth Pro, посетите http://earth.google.com/