Этот пример показывает, как вычислить отношения 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/