Этот пример показывает анализ источника и диффузию эпидемии SARS. Это основано на обсуждении вирусной филогении, представленной в Главе 7 "Введения в Вычислительную Геномику. Подход Тематических исследований" [1].
SARS (Тяжелый Острый Респираторный синдром) является недавно появившейся болезнью, вызванной новым типом коронавируса (SARS-CoV). Это состоит из 29 571 основного длинного, одноцепочечного RNA и отображает характеристический остроконечный белок конверта, который напоминает корону.
Первые случаи SARS появились в последних 2002 в китайской провинции Гуандун и превратились в основную вспышку за следующие несколько месяцев (январь в течение февраля 2003). Большинство зараженных индивидуумов получило болезнь в Больнице Гуанчжоу. Доктор, который работал в этой больнице, переместился в Гонконг в феврале 2003 и остановился в отеле Metropole. Доктор и много других гостей отеля все стали зараженными вирусом и переместились в различные места назначения (Вьетнам, Канада, Сингапур, Тайвань) явление носителем болезни и вируса.
Путем анализа филогенетических отношений между выборками SARS-CoV, которые были собраны в последних 2002 и в 2 003, мы можем восстановить историю эпидемии SARS и изучить, как это распространилось во всем мире за такой короткий период времени.
Мы рассматриваем последовательности нуклеотида 13 деформаций человеческих коронавирусов SARS, которыми известны местоположение и дата набора. Последовательности соответствуют скачку S белок, который ответственен за привязку с определенными приемниками и рассматривается главным аллергенным определителем. Поскольку гималайская пальмовая куница, как полагают, является источником человеческого SARS-CoV, мы также считаем выборку выведенной из пальмовой куницы. Ради удобства данные о последовательности хранятся в структуре MATLAB® под названием spike
состоя из Header
и Sequence
поля для каждой вирусной деформации. Данные могут также быть загружены с базы данных GenBank® с помощью инвентарных номеров, сохраненных в таблице accNum
.
% Load the genomic data for the human and palm civet SARS strains load sarsdata.mat % Display the accession numbers and collection dates of the sequence % dataset. accNum
accNum = 14x3 table GenbankAccession CollectionDate Location ________________ _______________ ______________________ {'AY278489'} {'DEC-16-2002'} {'GZ 12/16/02' } {'AY394997'} {'DEC-26-2002'} {'ZS 12/26/02' } {'AY395004'} {'JAN-04-2003'} {'ZS 01/04/03' } {'AY394978'} {'JAN-24-2003'} {'GZ 01/24/03' } {'AY394983'} {'JAN-31-2003'} {'GZ Hospital' } {'AY304495'} {'FEB-18-2002'} {'GZ 02/18/03' } {'AY278554'} {'FEB-21-2003'} {'Metropole 02/21/03'} {'AY278741'} {'FEB-26-2003'} {'Hanoi 02/26/03' } {'AY274119'} {'FEB-27-2003'} {'Toronto 02/27/03' } {'AY283794'} {'MAR-01-2003'} {'Singapore 03/01/03'} {'AY291451'} {'MAR-08-2003'} {'Taiwan 03/08/03' } {'AY345986'} {'MAR-19-2003'} {'Hong Kong 03/19/03'} {'AY394999'} {'MAY-15-2003'} {'Hong Kong 05/15/03'} {'AY627048'} {' '} {'Palm civet' }
Получите матрицу расстояния, должен был создать филогенетическое дерево путем вычисления симметрической матрицы, которая содержит попарные генетические расстояния с коррекциями Jukes-Cantor. Проигнорируйте сайты последовательности, представляющие разрывы.
JC_distances = seqpdist(spike,'method','jukes-cantor','alphabet','NT', ... 'indels','pairwise-delete','squareform',true); numSeq = size(JC_distances,1);
Путем графического вывода матрицы расстояния мы можем ценить присутствие подмножества последовательностей, которые более тесно связаны друг с другом (центральный кластер, представленный более темными тонами). Последняя последовательность, которая сопоставлена гималайской пальмовой кунице, является самой удаленной большинству членов набора. Это ожидается, потому что это - нечеловеческий коронавирус.
figure;
imagesc(JC_distances);
colormap(bone);
colorbar;
title('Pair-wise distances (spike protein nt sequences)');
Используя расстояния, вычисленные выше, создайте филогенетическое дерево с помощью соединяющего соседа метода. В этом случае мы принимаем равное отклонение и независимость эволюционных оценок расстояния.
tree1 = seqneighjoin(JC_distances,'equivar',spike); plot(tree1,'orient','left'); title('Neighbor-joining tree using Jukes-Cantor model');
Дерево изображает историю эпидемии. Ранние заражения все произошли в Гуанчжоу и Чжуншане, помеченном как GZ и ZS соответственно. Международные случаи (Гонконг, Сингапур, Ханой, Тайвань, Торонто) все связаны друг с другом и, кажется, переходят от случая, прослеженного до отеля Metropole в Гонконге.
Поскольку дата набора каждой деформации SARS известна, мы можем наблюдать прогрессию вирусных мутаций в зависимости от времени. Рассмотрите попарные расстояния согласно модели Kimura, которая различает переходные и tranversional уровни мутации. Затем ограничьте свой анализ расстоянием каждой человеческой деформации от деформации гималайской пальмовой куницы. Наконец, постройте генетическое расстояние по сравнению с датой набора.
K_distances = seqpdist(spike,'method','Kimura','squareform',true, ... 'alphabet','NT','indels','pairwise-delete'); % sequence of the palm civet civet = find(~cellfun(@isempty,strfind({spike.Header},'civet'))); % compute the genetic distance with respect to the palm civet's strain scores = zeros(numSeq-1,1); dates = zeros(numSeq-1,1); d = regexp({spike.Header},'\d+/\d+/\d+','match','once'); for i = 1:numSeq-1 % genetic distances with respect to the palm civet's strain scores(i,1) = K_distances(civet,i); % convert the collection dates into numbers dates(i,1) = datenum(d{i}); end refDate = datenum('01/01/03','mm/dd/yy'); % reference date figure(); plot(dates-refDate,scores,'k*'); ylabel('Genetic distance (relative to the palm civet)'); xlabel('Time distance from 01/01/03 (days)'); hold on;
Относительно последовательности пальмовой куницы мы замечаем, что генетическое расстояние увеличивается приблизительно линейным способом со временем. Выполните полиномиальный подбор кривой и интерполяцию наименьшего квадрата, чтобы обрисовать в общих чертах прогрессию вирусных мутаций в зависимости от времени и оценить аппроксимированную дату источника эпидемии. Запуск заражения соответствует более или менее корню аппроксимации полиномом, i.e., любая дата, которая является на нулевом генетическом расстоянии от последовательности пальмовой куницы.
[P,S] = polyfit(dates-refDate,scores,1); x = -max(dates-refDate):.1:max(dates-refDate); [y,delta] = polyconf(P,x,S); % estimate 95% prediction interval plot(x,y,'b-'); plot(x,y+delta,'r-',x,y-delta,'r-'); % confidence interval line([-max(dates-refDate) max(dates-refDate)],[0 0],'LineStyle',':'); title('Estimate of origin of SARS epidemic'); originDist = roots(P); % estimated distance between origin and reference date estimated_origin = datestr(floor(originDist+refDate)) plot(originDist,0,'*b'); annotation(gcf,'textarrow',[0.245 0.245],[0.30 0.35], ... 'String',{'estimated origin'},'color',[0 0 1]);
estimated_origin = '17-Sep-2002'
Поскольку болезнь, вызванная новой деформацией человеческого SARS-CoV, кажется, произошла у пальмовой куницы, мы можем принять, что местоположение корня для филогенетического дерева человеческих деформаций рядом с узлом, сопоставленным с гималайской пальмовой куницей.
civetNode = getbyname(tree1,'civet'); tree2 = reroot(tree1,civetNode,0); plot(tree2,'orient','left'); title('Rerooted Neighbor-joining tree using Jukes-Cantor model');
Перевнедренное дерево лучше иллюстрирует прогрессию эпидемии SARS. Начиная с ранних заражений в провинции Гуандун в 2 002 (GZ 12/16/02 и ZS 12/26/02), вирусное распространение в Больнице Гуанчжоу в ранних 2003 (Больница GZ 01/31/03) и достигнутый Гонконг через доктора, который работал в упомянутой больнице и кто остановился в отеле Metropole (Metropole 02/21/03). Вирус затем несли через границы через зараженных гостей отеля Metropole.
Предположение, что выборки представляют коронавирус SARS в различных моментах времени, мы можем наблюдать вирусную эволюцию, когда филогенетическое дерево (созданный на основе генетических расстояний) создается. Мы можем симулировать различные шаги в древовидной реконструкции. movie
функция анимирует процесс создания дерева.
d = regexp({spike.Header},'\d+/\d+/\d+','match','once'); d{end} = datestr(estimated_origin,'mm/dd/yy'); allDates = datenum(d); [dummy,order] = sort(allDates); % sort according to collection date for i = 2:numSeq sp = order(1:i); tr1 = seqneighjoin(JC_distances(sp,sp),'equivar',spike(sp)); tr2 = reroot(tr1,getbyname(tr1,'civet'),0); h = plot(tr2,'leaflabels',true,'terminallabels',false); h1 = findobj(h.leafNodeLabels,'string',spike(sp(i)).Header); h1.Color = 'r'; axis([-.0002 .0045 0 15]) fs(i-1) = h.axes.Parent; M(i-1) = getframe(fs(i-1)); end close(fs) % close figures
% movie(figure,M,1,1) % <== uncomment this line to play the animation
Мы можем также визуализировать диффузию вируса с помощью ориентированного графа, где каждый узел представляет зараженного индивидуума, и веса ребер сопоставлены к генетическому расстоянию между последовательностями. Во-первых, создайте матрицу смежности на основе даты набора выборок, таких, что возможные пути пробегают узлы, которые совместимы в терминах дат набора. Затем используйте ранее вычисленные расстояния Jukes-Cantor, чтобы присвоить веса ребрам между узлами. И наконец, определите кратчайший путь из узла, сопоставленного с пальмовой куницей и любым узлом.
% adjacency matrix based on collection dates gValid = bsxfun(@lt,allDates,allDates'); % weight matrix for the graph g1 = sparse((gValid .* JC_distances)); % find shortest paths from civet node to all nodes [dist,paths,pred_tree] = graphshortestpath(g1,civet); % create adjacency matrix for the winning shortest path g2 = sparse(pred_tree(1:13),1:13,1,14,14).*g1; % plot the graph spikeGraph = view(biograph(g2,{spike.Header})); % nodes relative to Guangdong province (GZ and ZS) guangdong = find((~cellfun(@isempty,strfind({spike.Header},'GZ'))) |... (~cellfun(@isempty,strfind({spike.Header},'ZS')))); % node relative to the Metropole Hotel metropole = find(~cellfun(@isempty,strfind({spike.Header},'Metropole'))); % node relative to the Guangzhou Hospital hospital = find(~cellfun(@isempty,strfind({spike.Header},'Hospital'))); % highlight some of the important nodes set(spikeGraph.Nodes(civet),'Color',[1 1 1]) % white (palm civet) set(spikeGraph.Nodes(guangdong),'Color',[1 .7 .7]) % pink (Guangdong) set(spikeGraph.Nodes(metropole),'Color',[0.8 0.8 1]) % lavander (Metropole Hotel) set(spikeGraph.Nodes(hospital),'Color',[1 0.3 0.3]) % red (GZ Hospital)
Этот график подсвечивает важную роль, которую играют некоторые события заражения:
Гималайская пальмовая куница является источником заражения
Отель Metropole является корнем ответвления для международной эпидемии
Больница Гуанчжоу представляет мост, соединяющий провинцию Гуандун (GZ и ZS) и отель Metropole в Гонконге.
[1] Cristianini, M. и Хан, M.W. "Введение в вычислительную геномику: подход тематических исследований", издательство Кембриджского университета, 2007.