Введение

ОРВИ (тяжелый острый респираторный синдром) - недавно возникшее заболевание, вызванное новым типом коронавируса (SARS-CoV). Он состоит из 29 571 одноцепочечной РНК длиной в основание и демонстрирует характерный белок остроконечной оболочки, который напоминает корону.

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

Анализируя филогенетические взаимосвязи между образцами SARS-CoV, которые были собраны в конце 2002 года и в 2003 году, мы можем реконструировать историю эпидемии атипичной пневмонии и понять, как она распространилась по всему миру за столь короткий период времени.

Загрузка данных последовательности штаммов SARS

Рассмотрим нуклеотидные последовательности 13 штаммов коронавирусов SARS человека, для которых известно местоположение и дата сбора. Последовательности соответствуют белку шипа S, который отвечает за связывание со специфическими рецепторами и считается основной антигенной детерминантой. Поскольку считается, что гималайская пальмовая цивета является источником человеческой атипичной пневмонии, мы также рассматриваем образец, полученный из пальмовой циветы. Для удобства данные последовательности хранятся в структуре 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'        }

Вычисление парных расстояний последовательности

Получите матрицу расстояний, необходимую для построения филогенетического дерева, вычисляя симметричную матрицу, которая содержит парные генетические расстояния с коррекциями Джукса-Кантора. Игнорировать сайты последовательности, представляющие промежутки.

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 известна, мы можем наблюдать прогрессирование мутаций вируса с течением времени. Рассмотрим парные расстояния по модели Кимуры, которая различает переходные и трансверсциональные скорости мутаций. Затем ограничьте свой анализ расстоянием каждого штамма человека от штамма циветы гималайской пальмы. Наконец, постройте график генетического расстояния по сравнению с датой сбора.

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;

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

[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');

Перезагруженное дерево лучше иллюстрирует прогрессирование эпидемии ОРВИ. Начиная с ранних инфекций в провинции Гуандун в 2002 году (GZ 12/16/02 и ZS 12/26/02), вирус распространился в больнице Гуанчжоу в начале 2003 года (GZ Hospital 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

Визуализация диффузии вируса через направленный график

Мы также можем визуализировать диффузию вируса с помощью направленного графа, где каждый узел представляет инфицированного человека, а веса рёбер связаны с генетическим расстоянием между последовательностями. Во-первых, создайте матрицу смежности на основе даты сбора выборок, так что возможные пути проходят через узлы, которые совместимы с точки зрения дат сбора. Затем используйте ранее вычисленные расстояния Джукса-Кантора для назначения весов ребрам между узлами. И, наконец, определите кратчайший путь от узла, связанного с пальмовой циветой, и каждого другого узла.

 % 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] Криштианини, М. и Хан, М.В. «Введение в вычислительную геномику: подход к тематическим исследованиям», Cambridge University Press, 2007.