Этот пример показывает, как рассчитать отношения Ka/Ks для восьми генов в геномах вируса H5N1 и H2N3 и выполнить филогенетический анализ гена HA из вируса H5N1 выделенного из цыплят в Африке и Азии. Для филогенетического анализа необходимо реконструировать дерево соединения соседей и создать 3-D график расстояний последовательностей с помощью многомерного масштабирования. Наконец, необходимо отобразить географические местоположения, в которых каждая последовательность HA была найдена на региональной карте. Последовательности, использованные в этом примере, были выбраны из исследования случая птичьего гриппа на веб-сайте вычислительной геномики [1]. Примечание.Для последнего раздела в этом примере требуется Toolbox™ мэппинга.
Есть три типа вируса гриппа: Тип A, B и C. Все геномы гриппа состоят из восьми сегментов или генов, которые кодируют для полимеразы B2 (PB2), полимеразы B1 (PB1), полимераза A (PA), гемагглютинин (ХА), нуклеобелок (NP), нейраминидаза (NA), матрица (M1) и неструктурные белки (NS1). Примечание: Вирус типа С имеет гемагглютинин-эстеразу (HE), гомолог НА.
Из трех типов гриппа тип А может быть наиболее разрушительным. Он поражает птиц (его естественный водоем), людей и других млекопитающих и является основной причиной глобальных эпидемий гриппа. Тип B поражает только людей, вызывающих местные эпидемии, и тип C не имеет тенденции вызывать серьезные заболевания.
Инфлюензы типа А далее классифицируют на различные подтипы в соответствии с вариациями аминокислотных последовательностей белков HA (H1-16) и NA (N1-9). Оба белка расположены на внешней стороне вируса. HA присоединяет вирус к клетке-хозяину, затем помогает в процессе слияния вируса с клеткой. NA отсекает вновь созданный вирус из клетки-хозяина, чтобы он мог перейти к здоровой новой клетке. Различие в аминокислотном составе в белке и рекомбинация различных белков HA и NA способствуют способности influenzas типа А прыгать виды-хозяева (т.е. птица человеку) и широкому диапазону тяжести. Многие новые препараты разрабатываются для нацеливания на белки HA и NA [2,3,4].
В 1997 году H5N1 подтип вируса птичьего гриппа, вирус гриппа типа А, совершил неожиданный скачок к людям в Гонконге, вызвав смерть шести человек. Для борьбы с быстро распространяющейся болезнью вся птица в Гонконге была уничтожена. Анализ последовательности вируса H5N1 показан здесь [2,4].
Исследование соотношений Ка/Кс для каждого сегмента гена вируса H5N1 даст некоторое представление о том, как каждый меняется с течением времени. Ка/К - отношение несинонимических изменений к синонимам в последовательности. Более подробное объяснение коэффициентов Ка/К см. в разделе Анализ синонимических и несинонимических коэффициентов замещения. Чтобы вычислить Ka/Ks, нужна копия гена из двух временных точек. Можно использовать вирус H5N1 выделенный из цыплят в Гонконге в 1997 и 2001 годах. Для сравнения можно включить вирус H2N3 выделенный из уток кряквы в Альберте в 1977 и 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, NA, M1 и NS1, необходимо удалить любые варианты сплайсинга.
Удалите варианты соединения с 1997 года H5N1
ntSeq97{7}(1) = [];% M2
ntSeq97{8}(1) = [];% NS2
Удалить варианты сращивания с 1977 года H2N3
ntSeq77{2}(2) = [];% PB1-F2
ntSeq77{7}(1) = [];% M2
ntSeq77{8}(1) = [];% NS2
Удалить варианты соединения с 1985 года H2N3
ntSeq85{2}(2) = [];% PB1-F2
ntSeq85{7}(1) = [];% M2
ntSeq85{8}(1) = [];% NS2
Необходимо выровнять нуклеотидные последовательности для вычисления отношения Ка/Кс. Выровнять белковые последовательности для каждого гена (доступно в поле «трансляция»), используя nwalign функция, затем вставить промежутки в нуклеотидную последовательность, используя seqinsertgaps. Используйте функцию dnds для расчета несинонимических и синонимичных показателей замещения для каждого из восьми генов в геномах вируса. Если вы заинтересованы в просмотре выравниваний последовательности, установите для опции «verbose» значение true при использовании 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
Примечание: Результаты отношения Ка/Кс могут отличаться от показанных на [1] из-за вариантов сращивания последовательностей.
Визуализация отношений Ка/Кс в гистограмме 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});
Создайте новый структурный массив, содержащий поля, соответствующие аминокислотной последовательности (Sequence) и исходной информации (Header). Можно извлечь исходную информацию из 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 с помощью метода Джукса-Кантора. Использовать seqneighjoin реконструировать филогенетическое дерево методом соединения соседей. Seqneighjoin возвращает phytree объект.
distHA = seqpdist(alignHA,'method','Jukes-Cantor','alpha','aa'); HA_NJtree = seqneighjoin(distHA,'equivar',alignHA);
Используйте view метод, связанный с phytree объекты для открытия дерева в инструменте филогенетического дерева.
view(HA_NJtree);
Другим способом визуализации взаимосвязи между последовательностями является использование многомерного масштабирования (MDS) с расстояниями, рассчитанными для филогенетического дерева. Эта функциональность предоставляется cmdscale функция в Toolbox™ статистики и машинного обучения.
[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);
Более подробную информацию об использовании отношений Ка/Кс, филогенетики и МДС для анализа последовательностей см. Cristianini и Hahn [5].
ПРИМЕЧАНИЕ: Для создания следующего рисунка требуется панель инструментов сопоставления.
С помощью инструментов из Mapping Toolbox можно построить график расположения каждого вируса на карте Африки и Азии. Для этого необходимы широта и долгота для каждого местоположения. Сведения о поиске геопространственных данных в Интернете см. в разделе Поиск геопространственных данных в Интернете. Широта и долгота для столицы каждого географического региона, где были выделены вирусы, приведены в этом примере.
создание структуры геострукта; regionHA, которая содержит географическую информацию для каждого отображаемого элемента или последовательности. Геострукт должен иметь поля Геометрия (Geometry), Лат (Lat) и Лон (Lon), которые определяют тип элемента, широту и долготу. Эта информация используется функциями сопоставления в панели инструментов сопоставления для отображения геопространственных данных.
[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);
Геострукт также может иметь поля атрибутов, содержащие дополнительную информацию о каждом элементе. Добавить поля атрибутов Имя и Кластер в regionHA структура. Поле Кластер содержит номер кластера последовательности, который будет использоваться для идентификации членства последовательностей в кластере.
[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')
ПРИМЕЧАНИЕ: Для экспорта данных в файл в формате KML требуется панель инструментов сопоставления.
Использование kmlwrite из панели инструментов сопоставления можно записать информацию о местоположении и аннотациях для каждой последовательности в файл в формате KML. Google Earth отображает географические данные из KML-файлов в браузере Earth. Панель инструментов сопоставления kmlwrite функция переводит геострукт, такой как regionHA, в файл в формате KML, который будет использоваться Google Earth. Для получения дополнительной информации о kmlwrite см. раздел Экспорт данных векторных точек в KML.
Можно далее аннотировать каждую последовательность информацией из раздела «Функции» файла GenBank с помощью featureparse функция. Затем эту информацию можно использовать для заполнения геострукта. regionHAи отображать его в табличной форме в виде тега описания для каждой метки в браузере Google Earth. В геострукте обязательными полями являются поля Геометрия (Geometry), Лат (Lat) и Лон (Lon). Все остальные поля считаются атрибутами 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 , которая используется для форматирования таблицы описания для каждого маркера. Структура атрибута определяет порядок и форматирование каждого атрибута. Его также можно использовать, чтобы не отображать один из атрибутов в геострукте. 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>';
Напишите 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 сохранялось с помощью меню Google Earth «Файл - > Сохранить изображение». Так выглядят данные в файле KML при загрузке в Google Earth. Чтобы этот вид переместился и увеличил масштаб региона над Азией.
Для просмотра информации о последовательности щелкните по метке. Регистрационный номер в каждой таблице данных является гиперссылкой на файл последовательности GenPept в базе данных белков NCBI.
При необходимости удалите новый KML-файл из папки вывода KML.
delete(filename)
close all
[1] https://computationalgenomics.blogs.bristol.ac.uk/case_studies/birdflu_demo
[2] Лейвер, В.Г., Бишофбергер, Н. и Вебстер, Р.Г., «Обезоруживающие вирусы гриппа», Scientific American, 280 (1): 78-87, 1999.
[3] Suzuki, Y. and Masatoshi, N., «Происхождение и эволюция генов вируса гриппа гемагглютинина», Molecular Biology and Evolution, 19 (4): 501-9, 2002.
[4] Гамбарян, А., и др., «Эволюция рецептор-связывающего фенотипа вирусов гриппа А (H5)», Вирусология, 344 (2): 432-8, 2006.
[5] Криштианини, Н. и Хан, М.В., «Введение в вычислительную геномику: подход к тематическим исследованиям», Cambridge University Press, 2007.
[6] Изображения Google Earth были получены с помощью Google Earth Pro. Дополнительную информацию о Google Earth и Google Earth Pro можно найти на сайте http://earth.google.com/