Практикум 9. EMBOSS, Entrez Direct, NCBI Datasets
Протеом бактерии Leptospira borgpetersenii serovar Ceylonica
Был скачен в прошлом семестре, но для напоминания вот код, с помощью которого это было сделано: wget 'https://rest.uniprot.org/uniprotkb/stream?compressed=true&format=txt&query=(proteome:UP000011783)' -O UP000011783.swiss.gz
Этап 1: получение AC геномной сборки и TaxID организма Leptospira borgpetersenii serovar Ceylonica
Идентификатор протеома: UP000263483
Ссылка на страницу в UniProt Proteomes: UniProt Proteomes: UP000263483
Ссылку на страницу сборки в NCBI Assembly (Datasets Genome): Genome assembly ASM351614v1
Данная ссылка была получена так: на странице в UniProt Proteomes взяла код геномной сборки из GenBank (GCA_003516145.1), затем в Datasets Genome можно найти соответсвующую геномную сборку и вставить ссылку на нее.
Код доступа геномной сборки Leptospira borgpetersenii serovar Ceylonica из GenBank: GCA_003516145.1
Код доступа геномной сборки Leptospira borgpetersenii serovar Ceylonica в RefSeq: GCF_003516145.1
Этап 2: скачивание последовательности генома и таблицы локальных особенностей Leptospira borgpetersenii serovar Ceylonica
Код для скачивания последовательности генома и feature table:
datasets download genome accession GCF_003516145.1 --include genome,gff3
Затем с помощью следующего кода был распакован архив:
unzip ncbi_dataset.zip
Этап 3: поиск и трансляция открытых рамок считывания
Для начала необходимо определить какой вариант генетического кода использует Leptospira borgpetersenii serovar Ceylonica.
Для нахождения NCBI Taxonomy ID использовали следующий код:
esearch -db taxonomy -query "Leptospira borgpetersenii serovar Ceylonica" | efetch -format xml
Получили NCBI Taxonomy ID: 508536
Информацию о варианте генетического кода можно получить в базе NCBI Тaxonomy, запустив следующий код:
efetch -db 'taxonomy' -id '508536' -format 'xml'
Таким образом, данная бактерия использует таблицу генетического кода №11
Следующим шагом было нахождение открытых рамок считывания (ORF) и получение их трансляции. Для этого был использован код приведенный ниже. Важно отметить, что рамки считывания были найдены только между стоп-кодонами (для исключения возможности трансляции с неправильнгого старт-кодона), также одним из условий была длина аминокислотной последовательности не меньше 50 аминокислот.
getorf -sequence GCF_003516145.1_ASM351614v1_genomic.fna -outseq protein.orf -table 11 -find 0 -minsize 150
Пояснение: -find 0 - искать рамки считывания между стоп-кодонами
Создание локальной базы данных proteome с использованием полученного файла:
makeblastdb -in protein.orf -dbtype prot -out proteome
Команда infoseq для проверки того, что все последовательности белков будут не короче 50 аминокислот:
infoseq -sequence protein.orf -only -length | awk '$1 < 50' #awk '$1 < 50' выведет значение в каждой строке, если первый элемент меньше 50 (у нас в каждой строке 1 значение), в файле protein.orf - аминокислоты; весь код не выдаст ничего, так как таких значений нет
Таким образом, была дополнительно проверена корректность выполнения задания.
Этап 4: получение последовательностей гомологичных метилтрансфераз
Для выполнения этого пункта была использована команда seqret и создан файл query.fasta, который содержит последовательности трех ДНК-метилтрансфераз: P0AED9 (Dcm, m5C-МТаза, E.coli), P0AEE8 (Dam, m6A-МТаза, E.coli), и P23941 (m4C-МТаза, Bacillus amyloliquefaciens)
echo "sw:P0AED9 sw:P0AEE8 sw:P23941" | tr ' ' '\n' | seqret -filter '@stdin' - outseq query.fasta
Этап 5: поиск по сходству последовательностей
С помощью метода blastp, с использованием раннее созданной локальной базы данных proteome был проведен поиск последовательностей, обладающих максимальным сходством
blastp -query query.fasta -outfmt 7 -db proteome -out blast_out
Пояснение: формат 7 = Tabular with comment lines
Название рамки находки с наибольшим весом: NZ_CP026671.1_19973
Координаты в геноме: 2009022 - 2009960
Вес находки: 62.0
Обнаружен гомолог: Dam, m6A-МТаза, E.coli (P0AEE8)
Для поиска координат в геноме использовался код:
grep 'NZ_CP026671.1_19973' protein.orf
Ссылка на табличную выдачу в формате 7: Выдача blastp
Из таблицы локальных особенностей были отобраны строки, соответствующие CDS, а также столбцы – с координатами (4 и 5), цепью (7) и дополнительной информацией (9).
grep 'CDS' genomic.gff | cut -f1,3,4,5,7,9 > CDS.tsv
Для проверки того, что были отобраны только CDS:
cut -f2 CDS.tsv | sort -u
Затем была добавлена строчка, которая соответствует моей находке, а также были взяты 3 сверху и 3 снизу CDS
echo -e 'NZ_CP026671.1\tCDS\t2009022\t2009960\t-\tFOUND-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'FOUND-ORF' > neighbors.tsv
Ссылка на таблицу: neighbors.tsv
Было найдено, что WP_016762643.1 (2009019 2009900) пересекается с исследуемой CDS (2009022 2009960), но координаты совпадают не полностью. Также стоит заметить, что WP_016762643.1 - это Dam family site-specific DNA-(adenine-N6)-methyltransferase, что полностью подтверждает находку.
Этап 6: поиск по аннотациям кодирующих участков
Поиск был осуществлен с помощью следующего конвейера edirec:
elink -db nuccore -id 'NZ_CP026671.1' -target protein| efilter -query '2.1.1.72[ECNO]' | efetch -format acc
Пояснение: elink получает id белков (целевая база данных), которые кодируются нуклеотидной последовательностью (исходная база данных, для последовательности указываем id).
efilter фильтрует, используя EC код
efetch скачивает результат в заданном формате (Accession Number)
Был найден белок: WP_016762643.1, то есть результат совпал с предыдущим.