# Вступление Мой выбор пал на ген простагландинсинтазы PTGS2 и его окружение. В таблице я привёл описательную информацию о выбранно фрагменте. | Предмет интереса| Величина | | ------ | :------: | |Идентификатор соотв. записи | NC_060925.1| |Номер хромосомы |1 | | 5'-coord| 186025753 | |3'-coord |186036111 | | Длина | 10358| | `.fasta`-файл| [near_PTGS2.fasta](https://drive.google.com/file/d/1gfRm5BOdmqks6wBUCAiVyz7c9RZSYJlV/view?usp=drive_link)| # Аннотация выбранного фрагмента В выбранный фрагмент целиком попадает PTGS2, кроме того, с 3'-конца мы видим первый экзон гена PACEER. Посмотрев на визуализацию гена PTGS2 в геномном браузере, понимаем: - ген кодируется цепью ДНК, комплементарной для данной - ген содержит 10-11 экзонов Но не кажется ли нам странным, что левые (на картинке) экзоны окрашены иначе, чем остальные? Мы сможем объяснить это наблюдение, если посмотрим на расширенную визуализацию: Теперь нам всё понятно - с гена могут транскрибироваться две мРНК: одна заканчивается на 3'-конце раньше, чем другая - видимо, здесь как-то хитро регулируется транскрипция. Таким образом, зелёная строка обобщает эти два транскрипта: тёмно-зелёным выделены участки, которые у них совпадают, салатовым - зоны различий. Таким образом, если учитывать PACEER, в выбранный фрагмент попадают 2 гена и 3 CDS (только что обсудили, для PTGS2 мы предполагаем два транскрипта). # Выбор таксона для поиска похожих последовательностей BLAST'ом Целевым для поиска я выбрал таксон *Euarchontoglires* - ветвь млекопитающих, состоящая из приматов и некоторых организмов, подобных грузунам. [Источник пикчи](https://genome.cshlp.org/content/17/6/760/F1.medium.gif) Теперь, когда мы станем запускать вариации BLAST'а, будем всякий раз искать по таксону *Euarchontoglires* (taxid:314146) и исключать из рассмотрения *Homo sapiens* (taxid:9606). # Megablast Итак, запустим самый "грубый" алгоритм, который найдёт нам самые явные совпадения, поскольку он использует якорь очень большого размера, в связи с чем им удастся вытащить только очень близкие последовательности. | Параметр| Значение| | ------ | :------: | |`max target sequences` | 5000 | |`word_size` | 28 | | `gap cost`| Linear | **Примечание про параметры:** я указал выбранный `gap_cost`, поскольку Linear уникален для Megablast'а. **Анализ выдачи:** Довольно много находок с покрытием поданного контига порядка 40-43%б идентичностью под 95% и evalue около машинного нуля, все они относятся к приматам. Такой query coverage, по-видимому, можно объяснить неконсервативностью интронов. При этом примечательно, что `max_score` у большинства этих находок больше `total_score`, что говорит о наличии нескольких участков совпадений. # Blastn Ищем нуклеотидный образец в нуклеотидной БД. Хороший пример задачи, для которой он отлично подходит, я привёл ниже в пункте про `blast+`. | Параметр| Значение| | ------ | :------: | |`max target sequences` | 1000 | |`word_size` | 15 | Увеличил размер слово и количество находок ради того, чтобы blastn вообще смог запуститься (при меньшей длине слова и большем количестве находок происходит ошибка на сервере). **Анализ выдачи:** На контрасте с megablast'ом, хотя здесь несколько больше находок с покрытием образца выше 50%, процент схожести для них находится в районе 70. Не все из этих последовательностей принадлежат приматам: две с наибольшим покрытием образца относятся к голому землекопу. Портрет этого красавчика можно взять [здесь](https://www.science.org/do/10.1126/science.aat1320/abs/rat_16x9.jpg)! В целом, результаты выходят более разнообразными. # Blastx Транслируем нуклеотидный запрос и ищем в белковой БД. Такой алгоритм логично применять в ситуации, когда нам интересно найти гомологи белка, который, вероятно, закодирован в исследуемой нуклеотидной последовательности. | Параметр| Значение| | ------ | :------: | |`max target sequences` | 1000 | | `word_size`| 3| **Анализ:** Поскольку интроны бывают довльно длинными, смотреть на `query_coverage` здесь бессмысленно: интроны не дают ему превзойти даже 25%. Процент схожести для большого кол-ва находок из топа (отсортированного по схожести) здесь порядка 80%. Ещё одна очень характерная для экзон-интронной особенность состоит во множественных находках в каждой из приведённых последовательностей: `total_score` тут больше `max_score` примерно в 5-6 раз, что даёт нижнюю оценку количества экзонов в гене (уж точно не менее 5). Наконец, для этого типа бласта у нас сравнительно большие `e-value`: если раньше они все были около 0, теперь они в районе 10^(-70). # tbastx Самый требовательный алгоритм: будет 6 раз бегать сквозь по-разному транслированные нуклеотидные БД. Этот вариант blast'а хорош для той же задачи, которую мы привели для blastx, с той лишь разницей, что предполагаемые находки находятся в неразмеченных (неаннотированных) геномах. | Параметр| Значение| | ------ | :------: | |`max target sequences` | 10 | | `word_size`| 3| Опять же, параметры поставил максимально простые, чтобы хотя бы запустилось. Так и не запустился:''\) # Blast+ В этом сюжете мы попробуем найти в T2T-геноме человека гомологи 16S- и 23S-рибосомальной РНК - первая входит в состав малой субъединицы рибосомы и содержит последовательность, комплементарную сайту Шайна-Дальгарно мРНК, вторая включена в больную субъединицу. Раз мы имеем дело с функциональными РНК, для нас не имеют смысла белковые и транслированные вариации BLAST'а, а ввиду большой филогенетической дистанции мы выбираем blastn, а не megablast, из возможных нуклеотидных алгоритмов. Итак, `sudo apt install ncbi-blast+` Проиндексируем файл `db.fna` командой: ``` makeblastdb -in db.fasta -dbtype nucl ``` Далее "бластанём" каждую рРНК отдельно. ``` blastn -task blastn -query 16s_ecoli.fasta -db db.fna -out 16s_blast.out -max_target_seqs 5000 -word_size 8 blastn -task blastn -query 23s_ecoli.fasta -db db.fna -out 23s_blast.out -max_target_seqs 5000 -word_size 8 ``` `word_size` 8 я взял, чтобы поиск не был мучительно долгим. А теперь посмотрим на результаты! **[Предполагаемые гомологи](https://drive.google.com/file/d/1DYTEd9XJSXJdIwrJdOeXohEdMEK-4EC6/view?usp=drive_link) 16S-рРНК E. Coli** Последовательности из топа по e-value (у них он порядка 10^(-7)) находятся на 13, 14, 15, 21, 22 хромосомах. **[Предполагаемые гомологи](https://drive.google.com/file/d/1rFJikoDj-e7Y8sDNQP6LkVQz5vmWgVWj/view?usp=drive_link) 23S-рРНК E. Coli** Лучшие находки имеют e-value порядка 10^(-17), лежат на 12, 13, 14, 15, 21, 22 хромосомах. Получилась довольно ожидаемая картина: гомологи нашлись на хромосомах, кодирующих рибосомальные РНК. Наверное, в другой сборке человеческого генома (в референсной, например), мы бы нашли ещё один очень характерный тип находок - рибосомальные РНК митохондрий, но в нашей T2T-сборке их пока нет.