Введение

Вирус иммунодефицита человека 1 (HIV-1) является более географически распространенным из двух вирусных штаммов, которые вызывают синдром приобретенного иммунодефицита (СПИД) у людей. Поскольку вирус быстро и постоянно развивается, в настоящее время нет ни лекарства, ни вакцины против ВИЧ-инфекции. Вирус ВИЧ представляет очень высокую скорость мутации, которая позволяет ему избегать реакции нашей иммунной системы, а также действия конкретных лекарств. В то же время, однако, быстрая эволюция ВИЧ обеспечивает мощный механизм, который показывает важное понимание его функции и устойчивости к лекарствам. Путем оценки силы селективного давления (положительный и очищающий отбор) в различных областях вирусного генома, мы можем получить общее понимание того, как вирус развивается. В частности, мы можем определить, какие гены развиваются в ответ на действие целевой иммунной системы и какие гены сохраняются, потому что они участвуют в некоторых важных функциях вируса.

Несинонимичные мутации в последовательности ДНК вызывают изменение в транслируемой аминокислотной последовательности, в то время как синонимические мутации не имеют. Сравнение количества несинонимических мутаций (dn или Ka), и количество синонимических мутаций (ds или Ks), могут предположить, действует ли на молекулярном уровне естественный отбор для содействия фиксации выгодных мутаций (положительный отбор) или для удаления вредных мутаций (очищающий отбор). В целом, когда доминирует положительный выбор, отношение Ka/Ks больше 1; в этом случае предпочтение отдается разнообразию на уровне аминокислот, вероятно, из-за преимущества пригодности, обеспечиваемого мутациями. И наоборот, когда преобладает отрицательный выбор, отношение Ka/Ks меньше 1; в этом случае большинство изменений аминокислот являются вредными и, следовательно, отбираются против. Когда положительные и отрицательные силы выбора уравновешивают друг друга, отношение Ka/Ks близко к 1.

Извлечение информации о последовательности для двух геномов HIV-1

Загрузите две геномные последовательности HIV-1 (номера доступа GenBank ® AF033819 и M27323). Для каждого кодированного гена мы извлекаем соответствующую информацию, такую как нуклеотидная последовательность, переведенная последовательность и имя генного продукта.

hiv1(1) = getgenbank('AF033819');
hiv1(2) = getgenbank('M27323');

Для вашего удобства ранее загруженные последовательности включаются в MAT-файл. Обратите внимание, что данные в общедоступных хранилищах часто кураторятся и обновляются; поэтому результаты этого примера могут несколько отличаться при использовании современных наборов данных.

load hiv1.mat

Извлеките информацию о последовательности генов с помощью featureparse функция.

genes1 = featureparse(hiv1(1),'feature','CDS','Sequence','true');
genes2 = featureparse(hiv1(2),'feature','CDS','Sequence','true');

Вычисление отношения Ka/Ks для генов HIV-1

Выровнять соответствующие белковые последовательности в двух геномах HIV-1 и использовать полученное выравнивание в качестве руководства для вставки соответствующих промежутков в нуклеотидные последовательности. Затем вычислите отношение Ka/Ks для каждого отдельного гена и постройте график результатов.

KaKs = zeros(1,numel(genes1));
for iCDS = 1:numel(genes1)
        % align aa sequences of corresponding genes
        [score,alignment] = nwalign(genes1(iCDS).translation,genes2(iCDS).translation);
        seq1 = seqinsertgaps(genes1(iCDS).Sequence,alignment(1,:));
        seq2 = seqinsertgaps(genes2(iCDS).Sequence,alignment(3,:));

        % Calculate synonymous and nonsynonymous substitution rates
        [dn,ds] = dnds(seq1,seq2);
        KaKs(iCDS) = dn/ds;
end

% plot Ka/Ks ratio for each gene
bar(KaKs);
ylabel('Ka / Ks')
xlabel('genes')
ax = gca;
ax.XTickLabel = {genes1.product};
% plot dotted line at threshold 1
hold on
line([0 numel(KaKs)+1],[1 1],'LineStyle', ':');
KaKs

KaKs =

  Columns 1 through 7

    0.2560    0.1359    0.3013    0.1128    1.1686    0.4179    0.5150

  Columns 8 through 9

    0.5115    0.3338

Все рассматриваемые гены, за исключением TAT, имеют общую Ka/Ks менее 1. Это соответствует тому, что большинство белкокодирующих генов считаются находящимися под эффектом очистки селекции. Действительно, большинство наблюдаемых мутаций являются синонимами и не влияют на целостность кодированных белков. В результате количество синонимических мутаций обычно превышает количество несинонимических мутаций. Случай TAT представляет собой хорошо известное исключение; на уровне аминокислот белок TAT является одним из наименее консервативных среди вирусных белков.

Вычисление коэффициента Ka/Ks с помощью раздвижных окон

Часто различные области одного гена могут подвергаться различным селективным давлениям. В этих случаях вычисление Ka/Ks по всей длине гена не дает подробной картины эволюционных ограничений, связанных с геном. Для примера общее значение Ka/Ks, сопоставленное с геном ENV, составляет 0,5155. Однако ген ENV кодирует оболочечную гликопротеиновую GP160, которая, в свою очередь, является предшественником двух белков: GP120 (остатки 31-511 в AF033819) и GP41 (остатки 512-856 в AF033819). GP120 подвергается воздействию на поверхности вирусной оболочки и выполняет первый этап ВИЧ-инфекции; GP41 нековалентно связан с GP120 и участвует во второй стадии ВИЧ-инфекции. Таким образом, мы можем ожидать, что эти два белка будут реагировать на различные селективные давления, и глобальный анализ всего гена ENV может скрыть диверсифицированное поведение. По этой причине мы проводим более мелкий анализ с помощью раздвижных окон разных размеров.

Выровнять гены ENV двух геномов и измерить отношение Ka/Ks по скользящим окнам размером 5, 45 и 200 кодонов.

env = 8; % ORF number corresponding to gene ENV

% align the two ENV genes
[score,alignment] = nwalign(genes1(env).translation,genes2(env).translation);
env_1 = seqinsertgaps(genes1(env).Sequence,alignment(1,:));
env_2 = seqinsertgaps(genes2(env).Sequence,alignment(3,:));

% compute Ka/Ks using sliding windows of different sizes
[dn1, ds1, vardn1, vards1] = dnds(env_1, env_2, 'window', 200);
[dn2, ds2, vardn2, vards2] = dnds(env_1, env_2, 'window', 45);
[dn3, ds3, vardn3, vards3] = dnds(env_1, env_2, 'window', 5);

% plot the Ka/Ks trends for the different window sizes
figure()
hold on
plot(dn1./ds1, 'r');
plot(dn2./ds2, 'b');
plot(dn3./ds3, 'g');
line([0 numel(dn3)],[1 1],'LineStyle',':');
legend('window size = 200', 'window size = 45', 'window size = 5');
ylim([0 10])
ylabel('Ka / Ks')
xlabel('sliding window (starting codon)')
title 'Env';

Выбор размера раздвижного окна может оказаться проблематичным: окна, которые являются слишком длинными (в этом примере 200 кодонов) средними значениями по длинным областям одного гена, таким образом скрывая сегменты, где Ka/Ks потенциально ведет себя своеобразно. Слишком короткие окна (в этом примере 5 кодонов), вероятно, приведут к результатам, которые являются очень шумными и поэтому не очень значимыми. В случае гена ENV, скользящее окно из 45 кодонов, по-видимому, является подходящим. На графике, хотя общий тренд ниже порога 1, мы наблюдаем несколько peaks над порогом 1. Эти области, по-видимому, подвергаются положительному отбору, который способствует разнообразию аминокислот, так как он обеспечивает некоторое преимущество соответствия.

Использование анализа скользящего окна для генов GAG, POL и ENV

Можно выполнить аналогичный анализ на других генах, которые показывают глобальное отношение Ka/Ks менее 1. Вычислите глобальное отношение Ka/Ks для генов GAG, POL и ENV. Затем повторите расчет с помощью скользящего окна.

gene_index = [1;2;8]; % ORF corresponding to the GAG, POL, ENV genes
windowSize = 45;

% display the global Ka/Ks for the GAG, POL and ENV genes
KaKs(gene_index)

for i = 1:numel(gene_index)
    ID = gene_index(i);
     [score,alignment] = nwalign(genes1(ID).translation,genes2(ID).translation);
    s1 = seqinsertgaps(genes1(ID).Sequence,alignment(1,:));
    s2 = seqinsertgaps(genes2(ID).Sequence,alignment(3,:));

    % plot Ka/Ks ratio obtained with the sliding window
    [dn, ds, vardn, vards] = dnds(s1, s2, 'window', windowSize);
    figure()
    plot(dn./ds, 'b')
    line([0 numel(dn)],[1 1], 'LineStyle', ':')
    ylabel('Ka / Ks')
    xlabel('sliding window (starting codon)')
    title(genes1(ID).product);
end

ans =

    0.2560    0.1359    0.5115

Ген GAG (Group-specific Antigen) обеспечивает базовую физическую инфраструктуру вируса. Он кодирует p24 (вирусный капсид), p6 и p7 (нуклеокапсидные белки) и p17 (матричный белок). Поскольку этот ген кодирует многие фундаментальные белки, которые структурно важны для выживания вируса, количество синонимических мутаций превышает количество несинонимических мутаций (то есть, Ka/Ks < 1). Таким образом, ожидается, что этот белок будет ограничен путем очистки отбора для поддержания вирусной инфекционности.

Ген POL кодирует вирусные ферменты, такие как обратная транскриптаза, интеграза и протеаза. Эти ферменты необходимы для выживания вируса, и поэтому селективное давление для сохранения их функции и структурной целостности довольно высокое. Следовательно, этот ген, по-видимому, находится под очистительным отбором, и мы наблюдаем значения отношения Ka/Ks менее 1 для большей части длины гена.

Ген ENV кодирует предшественник GP120 и GP41 белков, встроенных в вирусную оболочку, которые позволяют вирусу присоединяться к клеткам-мишеням и сливаться с ними. GP120 заражает любую целевую камеру путем связывания с CD4 рецептором. Как следствие, GP120 должны поддерживать механизм распознавания хозяина камеры и в то же время избегать обнаружения иммунной системой. Эти две роли выполняются различными частями белка, что показано трендом в отношении Ka/Ks. Этот вирусный белок подвергается очистке (Ka/Ks < 1) и положительному отбору (Ka/Ks > 1) в различных областях. Аналогичный тренд наблюдается и у GP41.

Анализ коэффициента Ka/Ks и эпитопов в GP120

Гликопротеин GP120 связывается с CD4 рецептором любой целевой камеры, в частности с Т-клеткой-хелпером. Это представляет собой первый этап ВИЧ-инфекции, и, следовательно, GP120 был одним из первых белков, изученных с целью нахождения вакцины против ВИЧ. Интересно определить, какие области GP120, по-видимому, подвергаются очищающему отбору, как индикаторы белковых областей, которые функционально или структурно важны для выживания вируса и потенциально могут представлять собой лекарственные мишени.

Из генов ENV извлекают последовательности, кодирующие GP120. Вычислите Ka/Ks по скользящему окну размера, равному 45 кодонам. Постройте и перекройте тренд Ka/Ks с расположением четырех эпитопов T камеры для GP120.

% GP120 protein boundaries in genome1 and genome2 respectively
gp120_start = [31; 30]; % protein boundaries
gp120_stop = [511; 501];
gp120_startnt = gp120_start*3-2; % nt boundaries
gp120_stopnt = gp120_stop*3;

% align GP120 proteins and insert appropriate gaps in nt sequence
[score,alignment] = nwalign(genes1(env).translation(gp120_start(1):gp120_stop(1)), ...
                   genes2(env).translation(gp120_start(2):gp120_stop(2)));
gp120_1 = seqinsertgaps(genes1(env).Sequence(gp120_startnt(1):gp120_stopnt(1)),alignment(1,:));
gp120_2 = seqinsertgaps(genes2(env).Sequence(gp120_startnt(2):gp120_stopnt(2)),alignment(3,:));

% Compute and plot Ka/Ks ratio using the sliding window
[dn120, ds120, vardn120, vards120] = dnds(gp120_1, gp120_2, 'window', windowSize);

% Epitopes for GP120 identified by cellular methods (see reference [2])
epitopes = {'TVYYGVPVWK','HEDIISLWQSLKPCVKLTPL',...
            'EVVIRSANFTNDAKATIIVQLNQSVEINCT','QIASKLREQFGNNK',...
            'QSSGGDPEIVTHSFNCGGEFF','KQFINMWQEVGKAMYAPP',...
            'DMRDNWRSELYKYKVVKIEPLGVAP'};

% Find location of the epitopes in the aligned sequences:
epiLoc = zeros(numel(epitopes),2);
for i = 1:numel(epitopes)
    [sco,ali,ind] = swalign(alignment(1,:),epitopes{i});
    epiLoc(i,:) = ind(1) + [0 length(ali)-1];
end

figure
hold on
% plot Ka/Ks relatively to the middle codon of the sliding window
plot(windowSize/2+(1:numel(dn120)),dn120./ds120)
plot(epiLoc,[1 1],'linewidth',5)
line([0 numel(dn120)+windowSize/2],[1 1],'LineStyle',':')
title('GP120, Ka / Ks and epitopes');
ylabel('Ka / Ks');
xlabel('sliding window (middle codon)');

Хотя общий тренд отношения Ka/Ks меньше 1, существуют некоторые области, где это соотношение больше единицы, что указывает на то, что эти области, вероятно, будут подвергнуты положительному отбору. Интересно, что расположение некоторых из этих областей соответствует наличию Т камеры эпитопов, идентифицированных клеточными методами. Эти сегменты демонстрируют высокую аминокислотную изменчивость, потому что аминокислотное разнообразие в этих областях позволяет вирусу уклоняться от распознавания иммунной системы хозяина. Таким образом, мы можем сделать вывод, что источником изменчивости в этих областях, вероятно, является иммунный ответ хозяина.

Ссылки

[1] Cristianini, N. and Hahn, M.W., «Introduction to Computational Genomics: A Case Studies Approach», Cambridge University Press, 2007.

[2] Siebert, S.A., et al., «Natural Selection on the gag, pol, and env Genes of Human Immunodefict Virus 1 (HIV-1)», Молекулярная биология и эволюция 12 (5): 803-813, 1995.