Практикум 9
Резюме: В ходе работы над данным практикумом был исследован геном Streptococcus thermophilus, на основе которого был получен протеом, который анализировался в 8 практикуме прошлого семестра (ID: UP000001170, количество белков: 1577).
1. Получение AC геномной сборки
Для начала нужно было из swiss файла узнать TaxID организма, геном которого будет анализироваться. TaxID прописан в поле OX, поэтому искать будем по соответствующим строкам. В моем файле были разные строки этого поля:
OX NCBI_TaxID=264199 {ECO:0000313|EMBL:AAV60042.1, ECO:0000313|Proteomes:UP000001170};OX NCBI_TaxID=264199;поэтому пришлось с ними справлятся. Конвейер:
grep '^OX' UP000001170.swiss | cut -f4 -d ' ' | tr -d ';' | uniq -c
Вывод (вроде нормально получилось):
1577 NCBI_TaxID=264199
Далее был получен список геномных сборок в человекочитаемом виде для данного TaxID с помощью следующей команды:
datasets summary genome taxon 264199 --as-json-lines | dataformat tsv genome | less -S
Нашлась одна геномная сборка в двух версиях: GenBank (GCA_000011825.1) и RefSeq (GCF_000011825.1). Для выполнения следующих заданий использовалась сборка RefSeq.
2. Скачивание последовательности генома и таблицы локальных особенностей
Далее нужно было скачать файлы, необходимые для дальнейшого анализа.
datasets download genome accession GCF_000011825.1 --include gff3,genome
unzip ncbi_dataset.zip
Получили архив, затем его распаковали. В архиве были таблица локальных особенностей в формате GFF3 и последовательность генома в формате fasta (еще два json файла, но они нам не нужны).
3. Поиск и трансляция открытых рамок считывания
Так как дальше нужно будет транслировать последовательность генома, сначала надо определить какой вариант генетического кода использует Streptococcus thermophilus. Это было сделано с помощью программы esummary:
esearch -db assembly -query "GCF_000011825.1" | elink -target nuccore | esummary
Наверное, можно как-то проще, но мне очень захотелось использовать esummary, а как это сделать лучше я не придумал. Зато получается удобно читаемый файл с вот такой строчкой:
<GeneticCode>11</GeneticCode>Поиск и трансляция открытых рамок считывания были осуществлены с помощью программы getorf:
getorf GCF_000011825.1_ASM1182v1_genomic.fna translated.fna -minsize 150 -table 11 -circular
(Ну у бакетрий же кольцевой геном, поэтому -circular). ORF по умолчанию понимается, как последовательность, свободная от стоп-кодонов (как и надо). С помощью команды infoseq было проверно, что все прошло корректно, и среди транслированных ORFs нет короче 50:
infoseq translated.fna -only -length | sort -n | less
Длина всех последовательностей больше или равна 50. По ним была создана база данных для blastp:
makeblastdb -in translated.fna -dbtype prot -out ORFs
4. Получение последовательностей гомологичных метилтрансфераз
Далее из swissprot были скачены последовательности трех ДНК-метилтрансфераз с заданными кодами доступа c помощью следующего конвейера:
echo 'sw:P0AED9 sw:P0AEE8 sw:P23941' | tr ' ' '\n' | seqret -filter 'list::stdin' -outseq query_MTases.fasta
Получили fasta файл с последовательностями белков Dcm (m5C), E.coli; Dam (m6A), E.coli и m4C-МТазы из Bacillus amyloliquefaciens
5. Поиск по сходству последовательностей
С помощью blastp по запросу из последовательностей метилтрансфераз был осуществлен поиск гомологов по базе данных ORFs. Команда:
blastp -query query_MTases.fasta -db ORFs -out blastp.out -outfmt 7
Самой лучшей находкой оказалась NC_006448.1_7383, координаты: 808572 - 806629 (нашел grep в translated.fasta); по m4C-МТазе из Bacillus amyloliquefaciens. Вес находки: 46.2, E-value: 2.49e-06. Вес вроде бы не очень большой, но E-value достаточно маленькое, поэтому сомнений не вызывает (кстати логично, что лучшая находка именно с белком из Bacillus amyloliquefaciens, так как они с Streptococcus thermophilus довольно близки по систематическому положению). Выдача blastp.
Далее нужно было по координатам находки найти CDS из таблицы локальных особенностей генома находятся рядом c ней. Сначала из таблицы локальных особенностей для CDS на соответствующей последовательности были вырезаны столбцы с координатами, указанием +/- цепь и аннотацией. Конвейер:
grep '^NC_006448.1' genomic.gff | grep 'CDS' | cut -f 4,5,7,9 >CDS.tsv
Следующим шагом были найдены близкие CDS:
echo -e '806629\t808572\t-\tMY-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'MY-ORF' >neighbors.tsv
Выдача: neighbors.tsv. Ура, одна из найденных CDS - это DNA methyltransferase (WP_011225879.1), правда запись с parent_gene (STU_RS13745) прекращена в NCBI. Но белковая есть.
6. Поиск по аннотациям кодирующих участков
Наконец, последним этапом стал поиск по аннотации кодирующих участков в геноме. Поиск проводился по EC-коду, который соответствует моей находке (m4C - 2.1.1.113). Конвейер:
elink -target 'protein' -db 'nuccore' -id 'NC_006448.1' | efilter -query {'2.1.1.72','2.1.1.113','2.1.1.37'} | efetch -format 'fasta' | less
Нашелся один белок - WP_011225759.1, type I restriction-modification system subunit M. И это не тот же, что был обнаружен ранее (тот кстати тоже restriction-modification system, но type III). И EC код у него, как у m6A, а не m4C. Возможно причина в плохой аннотации (для ранее найденного белка просто не прописан EC-код, поэтому его не находит). По аннотации эту CDS мы бы не обнаружили.