Введение

SARS (Тяжелый Острый Респираторный синдром) является недавно появившейся болезнью, вызванной новым типом коронавируса (SARS-CoV). Это состоит из 29 571 основного длинного, одноцепочечного RNA и отображает характеристический остроконечный белок конверта, который напоминает корону.

Первые случаи SARS появились в последних 2002 в китайской провинции Гуандун и превратились в основную вспышку за следующие несколько месяцев (январь в течение февраля 2003). Большинство зараженных индивидуумов получило болезнь в Больнице Гуанчжоу. Доктор, который работал в этой больнице, переместился в Гонконг в феврале 2003 и остановился в отеле Metropole. Доктор и много других гостей отеля все стали зараженными вирусом и переместились в различные места назначения (Вьетнам, Канада, Сингапур, Тайвань) явление носителем болезни и вируса.

Путем анализа филогенетических отношений между выборками SARS-CoV, которые были собраны в последних 2002 и в 2 003, мы можем восстановить историю эпидемии SARS и изучить, как это распространилось во всем мире за такой короткий период времени.

Загрузка данных о последовательности деформаций 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)');

Построение присоединения соседа Phylogenetic Tree

Используя расстояния, вычисленные выше, создайте филогенетическое дерево с помощью соединяющего соседа метода. В этом случае мы принимаем равное отклонение и независимость эволюционных оценок расстояния.

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'

Повторно базируясь Phylogenetic Tree

Поскольку болезнь, вызванная новой деформацией человеческого 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.

При наблюдении Phylogenetic Tree, когда это создает

Предположение, что выборки представляют коронавирус 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)

Этот график подсвечивает важную роль, которую играют некоторые события заражения:

Ссылки

[1] Cristianini, M. и Хан, M.W. "Введение в вычислительную геномику: подход тематических исследований", издательство Кембриджского университета, 2007.