Этот пример показывает, как анализ синонимичных и несинонимичных мутаций на уровне нуклеотида может предложить шаблоны молекулярной адаптации в геноме HIV-1. Этот пример основан на обсуждении естественного отбора на молекулярном уровне, представленном в Главе 6 "Введения в Вычислительную Геномику. Подход Тематических исследований" [1].
Вирус иммунодефицита человека 1 (HIV-1) более географически широко распространен из двух вирусных деформаций, которые вызывают Синдром приобретенного иммунодефицита (AIDS) в людях. Поскольку вирус быстро и постоянно развивается, в данный момент нет никакого средства исправления, ни вакцины против ВИЧ-инфекции. Вирус иммунодефицита человека представляет очень высокий уровень мутации, который позволяет ему уклоняться от ответа нашей иммунной системы, а также действия определенных наркотиков. В то же время, однако, быстрая эволюция HIV обеспечивает мощный механизм, который показывает важное понимание его функции и сопротивления наркотикам. Путем оценки силы селективных давлений (положительный и очищения выборов) через различные области вирусного генома, мы можем получить общее понимание того, как вирус развивается. В частности, мы можем определить, какие гены развиваются в ответ на действие целенаправленной иммунной системы и какие гены сохраняются, потому что они вовлечены в часть вируса существенные функции.
Несинонимичные мутации к последовательности DNA вызывают изменение в переведенной последовательности аминокислот, тогда как синонимичные мутации не делают. Сравнение между количеством несинонимичных мутаций (dn или Ka) и количеством синонимичных мутаций (ds или Ks), может предложить, действует ли на молекулярном уровне естественный отбор, чтобы способствовать фиксации выгодных мутаций (положительный выбор) или удалить вредные мутации (очищающий выбор). В целом, когда положительный выбор доминирует, отношение Ka/Ks больше, чем 1; в этом случае разнообразие на уровне аминокислоты одобрено, вероятно, из-за преимущества фитнеса, обеспеченного мутациями. С другой стороны, когда отрицательный выбор доминирует, отношение Ka/Ks - меньше чем 1; в этом случае большинство изменений аминокислоты вредно и, поэтому, выбрано против. Когда положительные и отрицательные силы выбора балансируют друг друга, отношение Ka/Ks близко к 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 геном и используйте получившееся выравнивание в качестве руководства, чтобы вставить соответствующие разрывы в последовательности нуклеотида. Затем вычислите отношение 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, сопоставленный с геном 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. Эти области, кажется, подвергаются положительному выбору, который способствует разнообразию аминокислоты, когда это обеспечивает некоторое преимущество фитнеса.
Можно выполнить подобные исследования других генов, которые отображают глобальное отношение 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 (Специфичный для группы Антиген) ген обеспечивает основную физическую инфраструктуру вируса. Это кодирует для 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.
Гликопротеин GP120 связывает с приемником CD4 любой целевой ячейки, особенно T-ячейка помощника. Это представляет первый шаг ВИЧ-инфекции и, поэтому, GP120 был среди первых белков, изученных с намерением нахождения вакцины HIV. Интересно определить, какие области 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, существуют некоторые области, где отношение больше, чем одно, указывая, что эти области, вероятно, будут являться объектом положительного выбора. Интересно, местоположение некоторых из этих областей соответствует присутствию антигенных детерминант ячейки T, идентифицированных сотовыми методами. Эти сегменты отображают высокую изменчивость аминокислоты, потому что разнообразие аминокислоты в этих областях позволяет вирусу уклоняться от распознавания иммунной системы хоста. Таким образом мы можем прийти к заключению, что источник изменчивости в этом области, вероятно, будет иммунной реакцией хоста.
[1] Cristianini, N. и Хан, M.W., "Введение в вычислительную геномику: подход тематических исследований", издательство Кембриджского университета, 2007.
[2] Siebert, S.A., и др., "Естественный отбор на затычке, политике и огибающих Генах Вируса иммунодефицита человека 1 (HIV-1)", Молекулярная биология и Эволюция, 12 (5):803-813, 1995.