Этот пример показывает, как вычислить отношения Ka/Ks для восьми генов в геномах вируса H5N1 и H2N3, и выполнить филогенетический анализ гена HA из вируса H5N1 выделенного из цыплят по всей Африке и Азии. Для филогенетического анализа вы восстановите соседнее дерево и создадите 3-D график расстояний последовательности с помощью многомерного масштабирования. Наконец, вы будете отображать географические местоположения, где каждый последовательность HA был найден на региональной карте. Последовательности, используемые в этом примере, были выбраны из тематического исследования птичьего гриппа на веб-сайте вычислительной геномики [1]. Примечание.В заключительном разделе этого примера требуется Mapping Toolbox™.
Есть три типа вируса гриппа: Тип A, B и C. Все геномы гриппа состоят из восьми сегментов или генов, которые кодируют для полимеразы B2 (PB2), полимеразы B1 (PB1), полимераза A (PA), гемагглютинин (ХА), нуклеобелок (NP), нейраминидаза (NA), матрица (M1) и неструктурные белки (NS1). Примечание: Вирус типа C имеет гемагглютининэстеразу (HE), гомолог HA.
Из трех типов гриппа тип А может быть наиболее разрушительным. Он поражает птиц (его естественный резервуар), людей и других млекопитающих и был основной причиной глобальных эпидемий гриппа. Тип B поражает только людей, вызывающих локальные эпидемии, и тип C не имеет тенденцию вызывать серьезные заболевания.
Инфлюэнзы типа А далее классифицируются на различные подтипы согласно изменений в аминокислотных последовательностях белков HA (H1-16) и NA (N1-9). Оба белка расположены снаружи вируса. HA присоединяет вирус к хозяину камеры затем помогает в процессе слияния вируса с камерой. NA клипирует недавно созданный вирус из камеры, чтобы он мог перейти к здоровой новой камере. Различие в аминокислотной композиции в белке и рекомбинация различных белков HA и NA способствуют способности инфлюэнза типа A скачкообразным видам хозяина (то есть птице к человеку) и широкой области значений тяжести. Многие новые лекарства разрабатываются для нацеливания на HA и NA белки [2,3,4].
В 1997 году H5N1 подтип вируса птичьего гриппа, вируса гриппа типа А, совершил неожиданный переход к людям в Гонконге, вызвав смерть шести человек. Чтобы контролировать быстро распространяющуюся болезнь, вся птица в Гонконге была уничтожена. Последовательный анализ вируса H5N1 показан здесь [2,4].
Исследование коэффициентов Ka/Ks для каждого генного сегмента H5N1 вируса даст некоторое представление о том, как каждый изменяется с течением времени. Ka/Ks - отношение несинонимических изменений к синонимам в последовательности. Для более подробного объяснения коэффициентов Ka/Ks смотрите Анализ Синонимических и несинонимических частот замещения. Чтобы вычислить 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 функция. The 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 генах, необходимо удалить любые варианты сращивания.
Удалите варианты сращивания из H5N1 1997 года
ntSeq97{7}(1) = [];% M2
ntSeq97{8}(1) = [];% NS2
Удалите варианты сращивания из H2N3 1977 года
ntSeq77{2}(2) = [];% PB1-F2
ntSeq77{7}(1) = [];% M2
ntSeq77{8}(1) = [];% NS2
Удалите варианты сращивания из H2N3 1985 года
ntSeq85{2}(2) = [];% PB1-F2
ntSeq85{7}(1) = [];% M2
ntSeq85{8}(1) = [];% NS2
Вам нужно выровнять нуклеотидные последовательности, чтобы вычислить отношение Ka/Ks. Выравнивание белковых последовательностей для каждого гена (доступно в поле 'translation') с помощью 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
Примечание: Результаты отношения 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});
Создайте новый массив структур, содержащий поля, соответствующие аминокислотной последовательности (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 объекты для открытия дерева в Phylogenetic Tree Tool.
view(HA_NJtree);
Другой способ визуализации связи между последовательностями - это использование многомерного масштабирования (MDS) с расстояниями, рассчитанными для филогенетического дерева. Эта функциональность обеспечивается cmdscale функция в Statistics and Machine Learning 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, эквивалентен поддереву, корневым узлом которого является Branch 14, который является Node 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, филогенетики и MDS для анализа последовательности см. Cristianini и Hahn [5].
ПРИМЕЧАНИЕ. Вам нужен Mapping Toolbox, чтобы создать следующий рисунок.
Используя инструменты из Mapping Toolbox, можно построить график местоположения, где был выделен каждый вирус, на карте Африки и Азии. Для этого нужны широта и долгота для каждого места. Для получения информации о поиске геопространственных данных в Интернете смотрите Find Geospatial Data Online. Для этого примера предусмотрены широта и долгота для столицы каждой географической области, где были выделены вирусы.
Создайте гео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 в своем браузере Earth. Mapping Toolbox функция переводит геоstruct, такой как regionHA, в файл в формате KML, который будет использоваться Google Earth. Для получения дополнительной информации о kmlwrite см. «Экспорт данных векторной точки в KML».
Можно далее аннотировать каждую последовательность с информацией из раздела Features файла GenBank с помощью featureparse функция. Затем можно использовать эту информацию для заполнения геоstruct, regionHA, и отобразить его в таблицу форме как тег-описание для каждого плацемента в браузере Google Earth. В геоstruct обязательными полями являются поля Геометрия, Лат и Лон. Все другие поля считаются атрибутами знака заполнения.
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>';
Напишите regionHA геоstruct в файл в формате 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] Laver, W.G., Bischofberger, N. and Webster, R.G., «Disarming Flu Viruses», Scientific American, 280 (1): 78-87, 1999.
[3] Suzuki, Y. and Masatoshi, N. «, Источник и эволюция генов гемагглютинина вируса гриппа», Молекулярная биология и эволюция, 19 (4): 501-9, 2002.
[4] Gambaryan, A., et al. «, Эволюция фенотипа связывания рецепторов вирусов гриппа A (H5)», Вирусология, 344 (2): 432-8, 2006.
[5] Cristianini, N. and Hahn, M.W., «Introduction to Computational Genomics: A Case Studies Approach», Cambridge University Press, 2007.
[6] Изображения Google Earth были получены с помощью Google Earth Pro. Для получения дополнительной информации о Google Earth и Google Earth Pro посетите http://earth.google.com/