Введение

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

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

Анализируя филогенетические связи между выборками SARS-CoV, которые были собраны в конце 2002 года и в 2003 году, мы можем восстановить историю эпидемии 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'        }

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

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

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

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

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'

Перезагрузка 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');

Перезагруженное дерево лучше иллюстрирует прогрессирование эпидемии атипичной пневмонии. Начиная с ранних инфекций в провинции Гуандун в 2002 году (GZ 12/16/02 и ZS 12/26/02), вирус распространился в больнице Гуанчжоу в начале 2003 года (GZ Hospital 01/31/03) и достиг Гонконга через врача, который работал в упомянутой больнице и который остановился в отеле Metropole (Metropole 02/21/03). Вирус затем переносился через границы через инфицированных гостей отеля Metropole.

Наблюдение за Phylogenetic Tree при его построении

Принимая, что выборки представляют коронавирус SARS в различные точки времени, мы можем наблюдать эволюцию вируса, когда создается филогенетическое дерево (построенное на основе генетических расстояний). Мы можем симулировать различные шаги в реконструкции дерева. The 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

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

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

 % 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. and Hahn, M.W. «Introduction to Computational Genomics: A Case Studies Approach», Cambridge University Press, 2007.