Для анализа генов бактерии Shigella flexneri 2a, которой принадлежит уже рассмотренный мною в прошлом семестре протеом, нужно сначала удостовериться, что протеом корректно собран и все белки принадлежат одному виду.
cat UP000001006.swiss | grep "TaxID" | cut -c5-20 | sort -u > IDs.txt
Теперь получим список сборок, принадлежащих этому таксону.
esearch -db 'assembly' -query '623 [TXID]' | efetch -format 'docsum' | grep 'AssemblyAccession>GCF' | cut -c24-36 | sort -u > ases
Ввиду того, что протеом собирался на нескольких сборках, будем искать сборку по ссылкам с нуклеотидных записей INCDS, на которые чаще всего ссылаются в протеоме.
cat UP000001006.swiss | grep 'DR EMBL;' | cut -d ';' -f2 | tr -d [:blank:] | sort | uniq -c | sort -n -r -k1 > out
head out -n1
elink -db 'nuccore' -id 'AE005674' -target 'assembly' | efetch -format 'docsum' | grep 'AssemblyAccession'
Теперь скачаем сборку и таблицу локальных особенностей при помощи полученного AC.
datasets download genome accession GCF_000006925.2 --include 'genome,gff3'
Распакуем скаченный файл.
unzip ncbi_dataset.zip
Для того чтобы корректно запустить blastp, необходимо узнать подробности генетичекого кода S. flexneri.
esearch -db 'taxonomy' -query '623 [UID]' | efetch -format 'xml' > gen_code.xml
Теперь найдём открытые рамки считывания, кодирующие не менее 50 аминокислот, и получим их трансляции при помощи программы getorf.
getorf -sequence ncbi_dataset/data/GCF_000006925.2/GCF_000006925.2_ASM692v2_genomic.fna -outseq orfs -find 0 -table 11 -minsize 150
Удостоверимся, что все получившиеся трансляты имеют длину не менее 50.
infoseq orfs -filter -only -length | sort -n > lengths.txt
Теперь можно на основе этих последовательностей создать базу данных.
makeblastdb -dbtype prot -in orfs -out ORFs
По данным кодам из SwissProt создадим файл-список для команд Emboss.
echo 'sw:P0AED9^sw:P0AEE8^sw:P23941' | tr '^' '\n' > SW.seq
Теперь получим последовательности из SwissProt.
seqret @SW.seq -filter > query_MTases.fasta
Запустим алгоритм blastp по созданной нами базе данных.
blastp -db ORFs -query query_MTases.fasta -out blastp_MT.fasta -evalue 0.001 -outfmt 7 -word_size 4
Алгоритм нашёл гомологи для всех белков. Две находки имели Evalue 0, и я выбрал ту, у которой был больше процент идентичности - гомолог m6A-МТазы E. coli.
Параметры:Вес-578,Evalue-0.0,Находка-NC_004337.2_24986
Координаты назодки получены следующим образом:
grep 'NC_004337.2_24986' orfs
Координаты:[3485129 - 3484287]
Теперь создадим файл CDS.tsv, содержащий строчки из таблици локальных особенностей, соответстыующие локусу NC_004337.2. И из этих строчек возьмём только некоторые строчки: 4,5-координаты, 7-направление цепи, 9-доп.инф.
grep 'NC_004337.2' ncbi*t/data/GCF*/*.gff | grep 'CDS' | cut -f4,5,7,9 > CDS.tsv
Теперь нужно в стандартном выводе добавить строчку с кординатами нашей находки, отсортировать стандартный вывод по первой координате и взять окрестность нашей строчки:
echo -e '3484287\t3485129\t-\tMY-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'MY-ORF' > neighbours.tsv
Нашлось пересечение с геном ДНК-аденин-метилазы, по видимому это и есть искомый ген гомолога. Параметры перекрывания приведены в таблице:
| Frame | Start | End | Strand |
|---|---|---|---|
| MyOrf | 3484287 | 3485129 | - |
| NP_709160.1 | 3484284 | 3485120 | - |
Как видно по таблице, кординаты не совсем совпадают. Открытая рамка найденная пограммой сдвинута немного вправо относительно гена.
Попробуем найти CDS через конвейер EDirect.
elink -db 'nuccore' -id 'NC_004337.2' -target 'protein' | efilter -query '2.1.1.72' | efetch -format 'fasta'
Однако найти ничего не удалось. Должно быть в описаниях белков моей бактерии не используется такая классификация ферментов.