Введение

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

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

Расчет отношения Ка/Кс для HIV-1 генов

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

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

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

Расчет отношения Ка/Кс с помощью скользящих окон

Часто различные области одного гена могут подвергаться различным селективным давлениям. В этих случаях вычисление 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, мы наблюдаем несколько пиков над порогом 1. Эти области, по-видимому, подвергаются положительному отбору, который способствует аминокислотному разнообразию, поскольку он обеспечивает некоторое преимущество в пригодности.

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

Можно выполнить аналогичный анализ для других генов, которые отображают глобальное отношение Ка/Кс менее 1. Вычислите глобальное отношение Ка/Кс для генов 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 кодирует вирусные ферменты, такие как обратная транскриптаза, интеграза и протеаза. Эти ферменты необходимы для выживания вируса, и, следовательно, селективное давление для сохранения их функции и структурной целостности является достаточно высоким. Следовательно, этот ген, по-видимому, подвергается очищающему отбору, и мы наблюдаем значения отношения Ка/Кс менее 1 для большей части длины гена.

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

Анализ отношения Ка/Кс и эпитопов в GP120

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

Из генов ENV извлекают последовательности, кодирующие GP120. Вычислите Ka/Ks над скользящим окном размером, равным 45 кодонам. Построить график и перекрыть тренд Ka/Ks с расположением четырех эпитопов Т-клеток для 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)');

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

Ссылки

[1] Криштианини, Н. и Хан, М.В., «Введение в вычислительную геномику: подход к тематическим исследованиям», Cambridge University Press, 2007.

Siebert, S.A., et al., «Естественный отбор на генах gag, pol и env вируса иммунодефицита человека 1 (HIV-1)», Molecular Biology and Evolution, 12 (5): 803-813, 1995.