Для позиций -3..4 от нуклеотида старта трансляции (для последовательности Козак) белок-кодирующих генов генома курицы Gallus gallus с помощью написанного скрипта была построена позиционно весовая матрица (PWM):

Комманда:

python3 pwm_build.py counts gallus_kozak.tsv 0.01 0.5

pseudoscore per nucl:   0.01
genome GC content:      0.5
PWM:
      -3     -2     -1      1      2      3      4
A  0.565  0.077 -0.511  1.386 -7.824 -7.824 -0.223
C -1.272  0.565  0.631 -7.824 -7.824 -7.824 -0.654
G  0.565 -0.223  0.215 -7.824 -7.824  1.386  0.751
T -1.608 -1.021 -1.272 -7.824  1.386 -7.824 -0.580

Для построения PWM была использована таблица счетов нуклеотидов для -3..4 позиций от старта трансляции. Псевдокаунты взяты равные 0.01 в силу наличия нулевых счетов только в позициях ATG кодона, где встреча других нуклеотидов (и реальное отклонение от 0) предполагается редким. GC состав генома – 50%.

Введение

Инициация трансляции у прокариот может происходить по 3 основным механизмам в зависимости от устройства 5'-UTR mRNA^[1]:

• SD-зависимая инициация – для mRNA, 5'-UTR которых несёт последовательность Shine-Dalgarno 10-15 nt upstream от старт-кодона; 30S рибосома связывается с SD последовательностью вблизи старт-кодона, 70S рибосома собирается уже на старт кодоне при наличии факторов инициации и инициаторной метионил-tRNA; характерна для бактерий и архей
• SD-независимая инициация – для mRNA c 5'-UTR без SD, но каким-либо другим сигналом для связывания 30S рибосомы; характерна для бактерий, где осуществляется за счёт связывания 5'-UTR mRNA рибосомальным белком RPB1, но видимо, существует с другим механизмом и у архей
• Безлидерная инициация – для mRNA с 5'-UTR короче 5 nt (вплоть до 0); mRNA 5'-фосфорилированным 5'-нуклеотидным остатком либо сразу заходит в mRNA-канал 70S рибосомы, либо связывается старт-кодоном при участии факторов инициации с 30S рибосомой, после чего добавляется и 50S рибосома; характерна для бактерий и особенно для архей, видимо, возможна и у эукариот.

Для дистальных цистронов полицистронных mRNA характерно 5'-UTR зависимая (SD-зависимая у архей) инициация трансляции, для проксимальных цистронов, так же как и для моноцистронных mRNA, возможны и лидерные, и безлидерные варианты инициации трансляции^[1,2].

Для бактерий доля CDS, проходящих инициацию трансляции по тому или иному типу, сильно различается в разных в филумах, тогда как у неэукариотических архей, кроме Euryarchaeota, преобладает безлидерный вариант. Наличие SD-зависимых CDS, однако, показано для представитей DPANN и TACK групп. Среди Euryarchaeota преобладают SD-зависимые CDS, за исключением Haloarchaea, использующих в основном безлидерную инициацию^[1,3].

Материалы и методы

Исследовалась сборка хромосомы Euryarchaeota археи Pyrococcus abyssi GE5 (NCBI AC AL096836.1).

Вся работа была проведена в Jupyter Notebook (загрузить файл проекта).

Построение PWM для SD-последовательностей P. abyssi

Получение первичной выборки для построения PWM:

1. С помощью скрипта features2CDSs.py из загруженного с NCBI файла feature table исследуемой сборки был получен файл со спиcком координат всех CDS сборки;
2. Были отобраны все CDS имеющие названия продуктов, и из них – те, в названии которых нет слов "hypothetical" и "putative", и при этом имеющие длину >= 300 nt – "хорошие" CDS. Такая процедура предполагает минимализацию вероятности столкнуться с неверно аннотированным старт-кодоном;
3. Из отобранных CDS случайной выборкой было получено 50 CDS;
4. Для полученной выборки из 50 CDS программой seqret из fasta файла сборки выделены (с учётом ориентации цепи) upstream-последовательности длиной 40 nt;
5. К upstream-последовательностям CDS была добавлена обратно комплементарная последовательность 3'-конца 16S rRNA длиной 40 nt в соответствии с координатами указанными в gbk файле сборки.

Собранные последовательности были переданы программе MEME:

meme upstreams.fasta -dna -oc . -nostatus -time 14400 -mod zoops\
-nmotifs 1 -minw 5 -maxw 5 -objfun classic -markov_order 0

Для многих CDS были найдены последовательности, характерные для SD мотива архей (GGTG[GA]), но не для 3'-конца 16S rRNA.

Поскольку без aSD-последоватеьности на 3'-конце 16S rRNA SD-последовательности CDS становятся нефункциональными, а для геномов Thermococci, и в частности P. abyssi, показан большой процент SD-несущих CDS и предполагается наличе активного механизма SD-зависимой инициации трансляции, было решено проверить запись о гене 16S rRNA исследуемой археи на неверную аннотацию.
Поиском blastn аннотированного участка 16S rRNA среди нуклеотидных последовательностей (NCBI Nucleotide collection) Thermococci была выбрано 12 последовательностей 16S rRNA, аннотированных как complete sequences. Программой MAFFT было построено выравнивание этих 12 последовательностей последовательностью 16S rRNA, с аннотированной в исследуемой сборке. Из полученного выравнивания видно, что во всех 12 выбранных последовательностях 3'-концевой участок несколько длиннее, чем в последовательности сборки. Было решено продлить 3' конец на 6 nt по сравнению с координатой в аннотации.

Upstream-последовательности CDS и обратно комплементарная последовательность 3'-конца 16S rRNA длиной 40 nt в соответствии с новыми координатами (теперь: 205039..206547) были переданы программе MEME:

meme upstreams.fasta -dna -oc . -nostatus -time 14400 -mod zoops\
 -nmotifs 1 -minw 5 -maxw 5 -objfun classic -markov_order 0       

В результате поиска мотив с паттерном GG[TA]G[GA] был обнаружен в 21 из 50 upstream-последовательностей CDS и в 3'-концевом участке 16S rRNA:

Полученная PWM мотива (транспонированная)

Найденный мотив подходит под консенсус последовательности SD архей GGTGA. E-value мотива оставляет желать лучшего: E-value = 1.9e+1, – то есть такой мотив мог вполне случайно найтись в наборе последовательностей данного размера и состава. Однако при случайном перемешивании нуклеотидов внутри последовательностей этот мотив не находится, позиции в пределах 10-15 nt upstream от старт-кодона тоже указывают на то, что находки не случайны. Здесь выделяется находка в upstream_32 необычной позицией в 24 nt upstream от старт-кодона; её было решено сохранить, ибо её последовательность GGAGA, в целом, подходящая под возможные варианты SD архей, может повысить чувствительность полученной PWM.

Из 50 отобранных upstream последовательностей SD был найден лишь в 21 последовательности (в 42%), что значительно меньше предыдущих оценок для P. abyssi^[3]. Это может быть связано с небольшим размером выборки, неправильной аннотацией старт-кодонов или недостаточностью алгоритма поиска. Поиск мотива одного типа длины 5 nt, встречающегося ровно 1 раз в каждой последовательности (параметр -mod oops), не даёт разумных результатов: находится похожий мотив, но в большом количестве ложноположительных сайтов.

Было замечено, что результат несколько улучшается при поиске мотива длиной от 5 до 9 nt, встречающегося не чаще 1 раза в каждой последовательности:

meme upstreams.fasta -dna -oc . -nostatus -time 14400 -mod zoops\
 -nmotifs 1 -minw 5 -maxw 9 -objfun classic -markov_order 0

При этом мотив длиной 9 nt и содержащий в себе участок GGTG[AG] находится в 25 из 50 upstream-последовательностей CDS, хотя у двух из них сайт необычно отдалён от определённого старт-кодона; также находка с хорошим p-value обнаруживается для 16S rRNA:

Полученная PWM мотива (транспонированная)

Мотив имеет хороший E-value = 5.4e-4 и находится в 25 из 50 отбранных upstream-последовательностей (в 50%), что всё равно меньше предыдущей оценки числа генов с SD-зависимой инициацией трансляции для P. abyssi^[3]. Более того, несмотря на хорошее совпадение мотива с 3'-концевой последовательностью 16S rRNA, для P. abyssi был показан более короткий SD мотив^[1].

В результате, было решено использовать обе построенные PWM для поиска генов с SD-зависимой инициацией трансляции в исследуемой хромосоме.

Поиск SD-последовательностей в геноме P. abyssi

I. Сначала был проведён тестовый поиск обеими простроенными PWM в 500 upstream-последовательностях длиной 20 nt из CDS, отобранных случайным образом из набора "хороших" белок-кодирующих генов (см. выше):

1. Из набора "хороших" CDS были случайно отобраны 500 записей, для них были определены координаты upstream-последовательнотей длиной 20 nt с учётом ориентации кодирующей цепи;
2. Из набора определённых на предыдущем шаге координат и fasta файла сборки хромосомы скриптом fragments2fasta.py был получен fasta файл с отобранными 500 upstream-последовательностями;
3. Поиск производился с последовательным применением обеих построенных PWM программой FIMO только на данных цепях сначала с p-value cutoff 0.001, затем – с p-value cutoff 0.01.

При поиске мотива GGTGA (далее: короткого мотива) в 500 upstream-последовательностях программой FIMO с p-value cutoff 0.001 было найдено всего 132 сайта данного мотива, 26.4% от числа входных последовательностей, что намного меньше ожидаемого числа CDS с SD-зависмой инициацией трансляции у взятой археи^[3]. Все найденные последовательности совпадают с консенсусом мотива.
При поиске FIMO с p-value cutoff 0.01 нашлось 450 сайтов данного мотива, последовательности сайтов стали более разнообразными. Дальнейший анализ проводился для этого набора сайтов:

Несмотря на разумный порог p-value, не все находки в наборе действительно значимы (multiplicity problem), поэтому решено было отобрать часть из них на основе q-values. При FDR = 0.05 было слишком много потерь, поэтому был принят FDR = 0.1.
Среди найденных сайтов были те, которые попадали дважды на одну upstream-последовательность, что биологически бессмысленно. Поэтому для каждой последовательности с множеством найденных сайтов была отобрана первая в таблице находка с наименьшим p-value (та же, что и с наибольшим score).

Была составлена таблица с названиями и координатами генов, в которых был найден сайт с мотивом, и указанным расстоянием от сайта до старт-кодона.

Таким образом, было найдено 344 сайта мотива, то есть 344 из 500 взятых CDS (68,8%) несут SD последовательность. Из распределения последовательностей найденных сайтов видно, что большинство соответствует консенсусу или близко к нему:

GGTGA    131
GGAGG     79
GGTGG     50
GGGGA     29
GGAGA     24
GGTGT     20
GGTGC      7
GGCGA      4

Для отобранных сайтов мотива построена гистограмма расстояний от III нуклеотида поледовательности мотива до старт-кодона (или точнее, бар-плот координат этого нуклеотида):

Для отобранных последовательностей сайтов мотива построено LOGO мотива (сервис WebLogo3, параметр no adjustment for composition):

При поиске мотива GGAGGTGAK (далее: длинного мотива) в 500 upstream-последовательностях программой FIMO с p-value cutoff 0.001 было найдено всего 169 сайтов данного мотива, 32.8% от числа входных последовательностей, что снова мало. Находки имели различные последовательности.
При поиске FIMO с p-value cutoff 0.01 нашлось 555 сайтов данного мотива. Дальнейший анализ проводился для этого набора сайтов:

При отборе наход по q-value, как и в прошлый раз был взят FDR = 0.1. Для избавления от множественных сайтов на upstream-последовательностях, из каждой группы таких сайтов снова выбирался первый в таблице сайт с наименьшим p-value.

Была составлена таблица с названиями и координатами генов, в которых был найден сайт с мотивом, и указанным расстоянием от сайта до старт-кодона.

Таким образом, было найдено 329 сайтов мотива, то есть 329 из 500 взятых CDS (65.8%) несут SD последовательность. Разнообразие найденных последоватеьностей велико (не более 8 сайтов один тип последовательности мотива). Было оценено количество характерных для SD мотива последовательностей в ожидаемо наиболее консервативном участке найденных сайтов (4..8 nt), а также количество сайтов, где в этом участке присутствует нехарктерный для SD мотива цитозин:

include GGTGA    				105
include GGTGA or GGAGG or GGTGG or GGGGA	174
include C					43 

Для отобранных сайтов мотива построена гистограмма расстояний от VI нуклеотида последовательности мотива до старт-кодона:

Для отобранных последовательностей сайтов мотива построено LOGO мотива (сервис WebLogo3, параметр no adjustment for composition):

Несмотря на хорошее сходство с последовательностью 3'-концевого участка 16S rRNA, длинная PWM находит в upstream-последовательностях 500 отобранных "хороших" CDS несколько меньше достоверных сайтов мотива, и значительно большая часть найденных сайтов имеет нехарактерные для искомого мотива последовательности по сравнению с поиском короткой PWM.
Для дальнейшего поиска мотива SD во всех белок-кодирующих генах хромосомы археи будет использоваться короткая PWM.

II. Затем был произведён поиск SD мотива, заданного короткой (более строгой) PWM, в upstream-последовательностях всех аннотированных CDS исследуемой сборки:

1. Для всех 1783 CDS сборки, координаты которых были получены ранее с помощью скрипта features2CDSs.py (см. выше), с учётом ориентации кодирующей цепи были определены координаты 20 nt upstream - 15 nt downstream последовательностей для поиска SD мотива, downstream участок взят для возможных случаев неверной аннотации старт-кодона;
2. Из набора определённых на предыдущем шаге координат и fasta файла сборки хромосомы скриптом fragments2fasta.py был получен fasta файл с отобранными 1783 upstream-downstream последовательностями;
3. Поиск производился применением короткой PWM программой FIMO только на данных цепях с p-value cutoff 0.001 и затем p-value cutoff 0.01.

При поиске короткого мотива с p-value cutoff 0.001 находилось всего 587 сайтов, что составляет 32.9% от числа входных последовательностей. Последовательности всех найденных сайтов совпадали с консенсусом искомого мотива.
При поиске с p-value cutoff 0.01 нашлось 2032 сайта данного мотива. Дальнейший анализ проводился для этого набора сайтов:

Для отбора по q-values был взят FDR = 0.1. Для избавления от множественных сайтов на upstream-последовательностях, из каждой группы таких сайтов выбирался первый в таблице сайт с наименьшим p-value.

Была составлена таблица с названиями и координатами генов, в которых был найден сайт с мотивом, и указанным расстояние от последнего нуклеотида сайта до аннотированного старт-кодона.

Таким образом, было найдено 1161 сайт мотива, то есть для 1161 из 1783 (65.1%) CDS хромосомы были обнаружены сайты SD мотива, что однако всё ещё меньше предыдущих оценок^[3].
Из распределения последовательностей найденных сайтов видно, что они все вполне похожи на последовательность SD:

GGTGA    560
GGAGG    247
GGTGG    240
GGAGA    114

Для отобранных сайтов мотива построена гистограмма расстояний от III нуклеотида поледовательности мотива до аннотированного старт-кодона (или точнее, бар-плот координат этого нуклеотида):

Как видно из гистограммы расстояний и таблицы найденных сайтов выше значительная часть сайтов была найдена downstream к аннотированному старт-кодону или "внутри" старт-кодона. Биологически это невозможно, объяснить такое состояние дел можно либо ложной находкой сайта SD, либо неверно аннотированным старт-кодоном.
Первое объяснение вполне вероятно: в силу нестрогого порога FDR до 10% наход могут быть ложноположительными. Однако "заскоков" за старт-кодон больше, чем 10% находок (23.4%) и распределение последовательностей этих находок вполне похоже на общее распределение (нет преобладания наименее характерных последовательностей):

GGTGA    130
GGTGG     57
GGAGA     51
GGAGG     34

Observed:

bad name	False	True
wrong start		
False		632	257
True		183	89

chi2-test p-value = 0.259

Таким образом, общую причину не удаётся выявить простыми методами.
В дальнейшем, при построении LOGO мотива считалось, что все найденные сайты определены верно, т.е. подразумевалось, что причина всегда в неверно аннотированном старт-кодоне (что, даже если неверно, вряд ли сильно скажется на LOGO из-за сходного распределения последовательностей).

Для отобранных последовательностей сайтов мотива построено LOGO мотива (сервис WebLogo3, параметр no adjustment for composition):

Выводы

В итоге, в хромосоме археи Pyrococcus abyssi GE5 последовательность SD была найдена для 65.1% белок-кодирующих генов, что значительно меньше определённого раньше значения 84.2% для этого же штамма, что, однако, всё равно отображает высокую распространённость SD-зависимой инициации трансляции у этой археи. Несовпадение с предыдущей оценкой связано с использовании в том исследовании более совершенного метода поиска SD с помощью рассчёта свободной энергии Гиббса связывания 5'-UTR mRNA с aSD 16S rRNA^[3]. Часть обнаруженных сайтов имела конфликт координат с аннотрированным старт-кодоном (расположены downstream от него), видимо, в некоторых случаях из-за ложности находки, а остальных случаях – из-за неверной аннотации старт-кодона для данного CDS.