Практикум 8. UniProt Proteomes, EMBOSS

1. Поиск протеома, соответствующего геномной сборке

Идентификатор геномной сборки моей бактерии Leptospira interrogans serovar Lai str. 56601 в RefSeq: GCF_000092565.1

Ссылка на страницу из базы NCBI Datasets Genome
https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000092565.1/

Идентификатор последней версии сборки INSDC: GCA_000092565.1

Поисковый запрос по UniProt Proteomes: (genome_assembly:GCA_000092565.1)
https://www.uniprot.org/proteomes?query=%28genome_assembly%3AGCA_000092565.1%29

Идентификатор протеома: UP000001408
https://www.uniprot.org/proteomes/UP000001408

Статус протеома: Reference proteome

2. Поиск и скачивание референсного протеома

Задание сформулировано следующим образом: вам нужно найти референсный протеом, наиболее близкий к протеому вашей бактерии/археи, и скачать белковые записи из UniProtKB, которые этому протеому принадлежат.

Найденный протеом уже референсный, но поскольку в задании сказано найти наиболее близкий, то я буду искать на уровне серотипа, и если ничего не найдется, то на уровне вида.

Поисковой запрос (proteome_type:1) AND (taxonomy_id:57678) дает мне тот же самый протеом UP000001408.

Поисковый запрос (proteome_type:1) AND (taxonomy_id:173) приводит к аналогичному результату.
В целом, здравый смысл подсказывает, что на нем и стоит остановиться. Для выполнения всех заданий будет использоваться именно этот протеом.

Краткая информация о протеоме: 3676 белков; BUSCO C:100% (S:100% D:0%) F:0% M:0%

Для скачивания была применена команда:

wget 'https://rest.uniprot.org/uniprotkb/stream?compressed=true&format=txt&query=(proteome: UP000001408)' -O ~/term2/pr8/UP000001408.swiss.gz

3. Оценка количества ферментов в протеоме

Моей задачей было попытаться оценить, сколько в протеоме белков, обладающих какой-либо ферментативной активностью, используя поисковые запросы на сайте UniProt и команды bash на kodomo.

3.1 Оценка количества ферментов в протеоме с помощью поисковых запросов

Поисковый запрос (proteome:UP000001408) AND (ec:*) позволяет найти все записи, имеющие какой-либо идентификатор ЕС.
(add field: Function → Enzyme classification [EC] → *)
Всего обнаружилось 585 записей.

Поисковый запрос (proteome:UP000001408) AND (ec:2) позволяет обнаружить все трансферазы моей бактерии. В результате было найдено 238 записей. Таким образом, можно подумать, что трансферазы составляют 40% от всех ферментов бактерии Leptospira interrogans serovar Lai str. 56601. (но на самом деле это не точно)

Поисковый запрос (proteome:UP000001408) AND ((ec:*) OR (keyword:KW-0560) OR (keyword:KW-0808) OR (keyword:KW-0378) OR (keyword:KW-0456) OR (keyword:KW-0413) OR (keyword:KW-0436) OR (keyword:KW-1278)) дал 1059 находок. Как и первый он позволяет найти все записи, имеющие какой-либо идентификатор ЕС, а также записи, в ключевых словах которых указаны 7 основных классов ферментов: оксидоредуктазы (EC 1), трансферазы (EC 2), гидролазы (EC 3), лиазы (EC 4), изомеразы (EC 5), лигазы (EC 6) или транслоказы (EC 7).

К последнему поисковому запросу можно еще что-нибудь добавить, например, protein_name, но я не буду этого делать. Приведенных поисковых запросов и их результатов вполне достаточно, чтобы увидеть разницу (585 vs 1059) по количеству записей почти в 2 раза и снова прийти к мысли о том, что с фильтрами нужно работать осознанно.

Если бы меня спросили, какой результат кажется мне более достоверным, я бы сказала что на общее количество 3676 белков 1059 обладающих ферментативной активностью мне кажутся более вероятным вариантом, нежели всего 585.

3.2 Оценка количества ферментов с помощью конвейера bash

Я планирую написать команду, которая подсчитывает количество строк в моем сжатом файле (по условию его распаковывать нельзя) базы данных белковых последовательностей (UP000001408.swiss.gz). Будут считаться строки, содержащие любой код классификации ферментов (EC) от 0 до 7 (по сути, это то же самое, что я делала ранее, но при помощи «сложного» поискового запроса).

Конвейер: zgrep 'EC=[0-7]' UP000001408.swiss.gz | wc -l

Результат: файл содержит 1198 таких строк, что более или менее согласуется с результатом последнего «сложного» поискового запроса.

Конвейер был составлен по аналогии с командами из подсказок:
zgrep 'что-нибудь' some-file.gz | sort -u
zcat some-file.gz | wc -l

4. Анализ протеома консольными средствами

Можно начать с подсчета ферментов каждого класса:

Конвейер: zgrep -o 'EC=[1-7]\.' UP000001408.swiss.gz | sort | uniq -с 

Результат:

Подсчет трансмембранных белков

Можно посчитать количество каких-то особых белков, например трансмембранных:

Конвейер: zcat UP000001408.swiss.gz | awk '/Transmembrane/{if(!skip) count++; skip=10} skip>0{skip--} END{print count}'

awk читает файл построчно и выполняет действия, когда строка соответствует шаблону
awk 'шаблон1{действие1} шаблон2{действие2} ...'

/Transmembrane/ - это шаблон (ищем строки, содержащие слово "Transmembrane")

Когда строка подошла под заданный шаблон, выполняется блок в {}:
if(!skip) - если переменная skip равна 0 (т.е. мы не в режиме пропуска)
count++ увеличиваем счётчик найденных уникальных вхождений на 1 (если строка соответствует шаблону, т.е. содержит слово "Transmembrane")
skip=10 переустанавливаем переменную skip в 10 (режим пропуска следующих 10 строк)

Перебираем строки дальше: блок skip>0{skip--}
skip>0 шаблон: «режим пропуска активен» (skip больше нуля: 10, 9, 8, 7...) то
{skip--} выполняется действие: skip уменьшается на 1 (так мы отсчитываем строки, которые нам осталось пропустить)

Переходим дальше (и так пока skip не обнулится)
блок END{print count}
END - специальный шаблон, который выполняется после обработки всех строк файла
print count выводит итоговое значение счётчика

Результат: 825 находок.

Конечно, такая команда - это своеобразный костыль, но просто считать, сколько раз встретилось слово "Transmembrane" не вариант, так как на один трансмембранный белок это слово может встретиться и 2, и 3 раза. Тогда и в итоговом результате количество трансмембранных белков будет сильно превышено (не будет соответствовать действительности). Лично я придумала (предварительно изучив структуру файла) просто пропускать следующие 10 строк после того, как попадается слово "Transmembrane", так что если оно и будет попадаться 2-3 раза на одну запись (белок), то счетчик его проигнорирует, и посчитает только первое соответствие шаблону.

Белки, связанные с патогенностью

Раз уж я когда-то выбрала для обзора патогенную лептоспиру, почему бы не посмотреть, сколько белков связаны с патогенностью (по ключевому слову Virulence):

Конвейер: zcat UP000001408.swiss.gz | grep -c '^KW.*Virulence'

Результат: 5 находок.

Можно попробовать узнать ID этих белков.

Конвейер: zcat UP000001408.swiss.gz | awk '/^ID/{id=$0; flag=0} /^KW/,/^$/{if(/Virulence /) flag=1} /^\//{if(flag) print id}'

Блок /^ID/{id=$0; flag=0}

/^ID/ шаблон: ищем строки, начинающиеся с "ID" (начало записи белка)

{id=$0; f=0} действие: id=$0 сохраняет всю текущую строку в переменную id (так мы запоминаем ID белка), flag=0 сбрасывает индикатор в 0 (так будет, пока не найдется заданное ключевое слово)

Дальше блок /^KW/,/^$/ {if(/Virulence /) flag=1}

Шаблон /^KW/,/^$/ это диапазон строк (действие начинается, когда мы встречаем строку, начинающуюся с ^KW, продолжается для всех следующих строк и прекращается, когда встречаем пустую строку (^$)

Действие {if(/Virulence /) flag=1} с условием: если в текущей строке есть слово Virulence с пробелом после, то ставим индикатор flag =1 (нашли нужное слово)

Зачем пробел после Virulence? Чтобы быть уверенным в том, что это самостоятельное слово.

Дальше блок /^\/\// {if(flag) print id}

Шаблон /^\/\// значит, что строка начинается с // (конец записи белка)

Действие {if(flag) print id} с условием: если индикатор не равен 0, выводится сохранённый ID белка

Результат: