Выполняя предыдущее задание Вы научились делать филогенетический анализ белок-кодирующих генов. Это относительно простая задача, ведь такие последовательности как правило достаточно консервативны и медленно накапливают мутации. Более того, по аминокислотной последовательности возможно реконструировать набор нуклеотидных последовательностей, способных кодировать исследуемый белок. Это позволяет использовать простые алгоритмы выравнивания и поиска гомологичных последовательностей, такие как tBLASTn. Однако белок-кодирующие элементы составляют лишь несколько процентов нашего генома, подавляющая его часть занята некодирующими и регуляторными последовательностями. И их филогенетический анализ представляет нетривиальную и весьма сложную задачу. Кстати, вот в этой статье как раз так и пишут, что поиск гомологов некодирующих последовательностей is tedious task. Публикация старая, но рекомендую её почитать, там пишут интересное.
NCBI — база данных, содержащая литературу, геномы, транскриптомы, протеомы различных живых организмов, нуклеотидные и аминокислотные последовательности генов и их белков и огромное количество другой информации. Также тут можно найти множество биоинформатических программ, доступных как через web-интерфейс, так и для скачивания
GenBank — это база данных NIH, коллекция всех аннотированных общедоступных последовательностей ДНК и РНК.
BLAST — Basic Local Alignment Tool — пакет программ, используемых для выравнивания нуклеотидных и аминокислотных последовательностей. Наиболее часто используемая программа в молекулярной биологии.
tBLASTn — высокая чувствительность при локальном выравнивании аминокислотных последовательностей
MAFFT — программа для множественного выравнивания
Ensembl — база данных, содержащая геномы, транскриптомы, протеомы различных живых организмов, нуклеотидные и аминокислотные последовательности генов и их белков, а ещё там есть BioMart, который поможет вам вытянуть почти любую информацию почти про любой ген.
Конвертер форматов — web версия программы для конвертации .fasta в .nexus.
Seqmagick — программа для конвертации .fasta в .nexus.
Lifemap — филогенетическое дерево живых организмов. Позволяет посмотреть взаимное эволюционное расположение видов из всех трёх доменов жизни. В строку поиска можно вводить латинское название вида.
Итак, вы научились проводить филогенетический анализ белок-кодирующих генов, а теперь вам предстоит сделать анализ некодирующих последовательностей. На самом деле разницы почти никакой нет, но тут потребуются несколько иные программы, хотя полюбившиеся вам MAFFT и Mr.Bayes останутся. Из-за того, что весьма чувствительный tBLASTn заточен под анализ аминокислотных последовательностей, а BLASTn не годится для поиска эволюционно удалённых гомологов, предлагаю вам использовать программу HMMER. HMMER основан на скрытых марковских моделях и обладает хорошими показателями чувствительности при поиски эволюционно дальних последовательностей, даже если они накопили значительное число мутаций.
Вам будут выданы ensembl id некодирующих последовательностей, для которых нужно будет провести филогенетический анализ. Для этого потребуется скачать их нуклеотидную последовательность из знакомой вам по предыдущему заданию базы ensembl.org. Затем при помощи HMMER найти ортологи исследуемого гена в геномах 11 видов живых организмов, перечисленных ниже. Получить нуклеотидные последовательности однаруженных ортологов при помощи ресурса BioMart. Для обнаруженных ортологов нужно провести множественное выравнивание программой MAFFT по аглоритму E-INS-i. E-INS-i наиболее оптимален для множественного выравнивания небольшого количества последовательностей с неизвестной заранее структурой. Полученный .fasta файл с результатом множественного выравнивания потребуется перевести в формат .nexus при помощи web версии конвертера или standalone пакета seqmagick. Теперь все готово для реконструкции филогенетического дерева. Это лучше делать на локально установленной версии Mr.Bayes, хотя допустимо использовать и web-версию. Однако web-версии зачастую имеют ограниченные функции и лимит на объём анализируемых данных. Полученное филогенетическое дерево можно посмотреть через FigTree или Archeopteryx, обе этих программы не требуют установки и запускаются напрямую из бинарного файла. Дерево лучше приукрасить, выбрав адекватные шрифты и кегль, раскрасив клады в разные цвета для удобства визуального восприятия. Работа выполнена, вы превосходны!
Выбрать ген или смириться с назначенным.
Найти на ensembl.org нуклеотидную последовательность наиболее длинной сплайс-формы, кодируемой исследуемым геном.
Сохранить выбранную нуклеотидную последовательность в тексnовом файле с расширением .fasta.
Скачать геномы предложенных ниже организмов с ensembl.org. Их лучше скачивать скриптом через rsync. Адрес однотипен и содержит латинское название вида, которое можно использовать, как переменную. Например, для шимпанзе путь к файлу с геномной последовательностью выглядит так:
/pub/release-110/fasta/pan_troglodytes/dna/*dna.toplevel.fa.gz
Скрипт может выглядеть как-то так:
#!/bin/sh
s_list='species.txt'
while read species; do
rsync -av rsync://ftp.ensembl.org/pub/release-113/fasta/${species}/dna/*dna.toplevel.fa.gz
done < ${s_list}
Но обратите внимание, что версия базы постоянно обновляется, поэтому вместо release-113 может понадобиться подставить более свежий релиз.
Преобразовать .fasta файлы в базу .hmm. HMMER в любом случае будет генерировать такую базу на основе геномной последовательности, но при каждом новом выравнивании заново, поэтому для экономии времени лучше сгенерировать её заранее и один раз. Чтобы не вводить команду 11 раз, удобнее сделать скрипт, взяв путь к геному в качестве переменной.
Если вы ленивы или нетерпеливы, то готовые геномные базы можно скачать тут.
Подождать, выпить чай. Возможно, поспать. Генерация .hmm базы занимает весьма много времени, а если компьютер не очень мощный, то придётся подождать пару часов. Можете посмотреть, если ещё не видели, The Martian.
Запустить выравнивание нуклеотидной последовательности на геном в программе nHMMER. Удобнее это делать скриптом, взяв путь к базе в качестве переменной. Описание команды запуска и назначение ключей описано ниже. Если вы ленивы, то скрипт можно взять тут. Пути в скриптах отредактировать нужно.
Подождать, выпить чай. Возможно, поспать. HMMER суров и молчалив, как норвежские боги, и не говорит о прогрессе выравнивания.
Среди обнаруженных последовательностей нужно выбрать только соответствующие пороговому значению по e-value (e-10 для нуклеотидных выравниваний, а вообще, ну прочитайте вы ту статью).
Если вы работали с геномами, то нужно получить геномные координаты подходящих последовательностей из TSV файла, перевести их в формат, пригодный для использования в BioMart — chromosome_name:start:end:strand, например, chr1:235783:236789:1 или chr2:235783:236789:-1.
Скачать нуклеотидные последовательности по координатам при помощи ресурса BioMart.
С высокой долей вероятности BioMart вам ничего не выдаст, потому что обнаруженный Вами гомолог не представлен в текущей аннотации генома конкретного организма. В этом случае можно воспользоваться программой bedtools, но для этого нужно преобразовать ваш .tsv вывод в .bed файл:
Структура .bed файла выглядит вот так:
chrom chromStart chromEnd name score strand
Колонка chrom — номер хромосомы, скаффолда или контига и соответствует колонке target name в выводном .tsv файле
Колонка chromStart — координаты начала таргетной последовательности и соответствует колонке alifrom в выводном .tsv файле
Колонка chromEnd — координаты конца таргетной последовательности и соответствует колонке ali to в выводном .tsv файле
Колонка chromEnd — название последовательности в .bed файле и, скорее всего, в вашем случае ничему не соответствует. Можете поставить 0. Но может соответствовать колонкам accession или description of target в выводном .tsv файле
Колонка score — может быть от 0 до 1000. И обозначать он может все, что угодно. Например, p-value. Можете сюда поместить значения колонки E-value или score из выводного .tsv файла. Вообще, .bed формат создан для UCSC Genome Browser, поэтому нам пригодятся не все поля
Колонка strand — направление цепи и соответствует колонке strand в выводном .tsv файле
Рекомендую воспользоваться командой awk
awk -v OFS='\t' '{print $1 $7 $8 $2 $13 $12}' output.tsv > file.bed
Вытащить последовательность по координатам через bedtools:
bedtools getfasta -fi <genome.fa> -bed <file.bed> -fo <output.fasta>
Сохранить полученные последовательности в едином файлк формата .fasta.
Теперь из заголовков в скачанном файле нужно убрать все символы кроме букв "A-Z", цифр "0-9", знака ">" в начале заголовка и underscore (_ — подчёркивание). Неподходящие символы можно заменить на underscore. Это нужно потому, что при дальнейшей конвертации в .nexus и запуске файла в Mr.Bayes возможно появление ошибок из-за наличия несъедобных символов. Также длина заголовка не должна превышать 99 символов.
Добавить в файл исследуемую Вами последовательность гена человека, не забыв отформатировать заголовок и в ней.
Запустить множественное выравнивание в MAFFT с использованием, например, алгоритма E-INS-i. Пример запуска смотрите внизу этой страницы.
Результат множественного выравнивания можно посмотреть в программе Jalview. Она кросплатформенная и очень простая. Правда, посмотрите, чтобы понять, как оно хотя бы выглядит.
Сконвертировать полученный .fasta файл в формат .nexus в web-конвертере или в seqmagick. При использовании seqmagick используйте команду
seqmagick convert --output-format nexus --alphabet dna YOUR_FILE.fasta YOUR_FILE.nexus
Провести реконструкцию филогенетического дерева в Mr.Bayes по параметрам, описанным ниже.
Открыть сгенерированное дерево в редакторе FigTree или Archeopteryx на ваш выбор, выбрав файл с расширением .con.tre.
Проанализировать дерево, его биологическую корректность.
Сделать вывод о времени возникновения гена. Тут вам, возможно, потребуется понять, какой вид к какому таксону принадлежит. Для этого можно зайти вот сюда. А если хочется погрузиться в филогению, то можно позалипать в Time Tree.
Отредактировать дерево, выбрав шрифт и кегль, покрасив клады в разные цвета для лучшего восприятия. Этот шаг опционален.
Для генерации .hmm файла введите в комендную строку makehmmerdb GENOME.fa GENOME.hmm.
Если хотите побыстрее, то воспользуйтесь скриптом (+10 to speed, requires 150 of intelligence).
Для Linux (даже, если это консольная Ubuntu на WSL) или MacOS, команда выглядит так:
nhmmer -E 1e-9 --tblout YOUR_GENE.tsv YOUR_GENE.fasta GENOME.hmm > YOUR_GENE.txt
А теперь, что все это значит:
Ключ -E 1e-9 указывает программе, что все совпадения с e-value выше 10-9 не будут записываться в вывод. Это нужно, чтобы не загромождать файл с результатами лишней информацией, которая и так будет отброшена.
Ключ --tblout создаёт файл, обозначенный YOUR_GENE.tsv и содержащий результаты выравнивания в виде .tsv файла, из которого удобно, естественно, при помощи скрипта вытащить геномные координаты.
GENOME.hmm — база HMMER, сгенерированная на основе геномной последовательности в формате .fasta.
YOUR_GENE.fasta — выбранная вами нуклеотидная последовательность
YOUR_GENE.txt — файл, содержащий результаты выравнивания с самим выравниванием — выровненные участки, гэпы, замены, координаты.
NCBI, National Center for Biotechnological Information. NCBI предоставляет информацию о базах данных ДНК (GenBank), белковых доменов, РНК, базах данных статей научной литературы (PubMed) и таксономической информации (TaxBrowser), обеспечивает поиск данных о конкретном биологическом виде (Taxonomy). Также содержит различные стандартные программы биоинформатики (BLAST).
GenBank — открытая база данных NCBI, содержащая многие аннотированные последовательности ДНК, РНК и белков. Исследователи, если описывают новую последовательность, часто выгружают её именно в GenBank.
BLAST - Basic Local Alignment Search Tool (буквально, средство поиска основного локального выравнивания) - семейство компьютерных программ, служащих для поиска гомологов белков или нуклеиновых кислот, для которых известна первичная структура (последовательность) или её фрагмент. Семейство программ, реализованных в NCBI BLAST. Представленная еще в 1990 г, программа стала своего рода гуглом биологии.
Query sequence (search sequence, input sequence) - анализируемая Вами последовательность.
Subject sequences – database sequences. Последовательности, депонированные в базе и используемые для сравнения с интересующей последовательностью.
Local alignment – identification of regions of similarity (“hits”) within long sequences that are often widely divergent overall; a global alignment contains all letters from both the query and target sequences. «Локальное выравнивание» - алгоритм, основанный на поиске коротких совпадающих участков (“hits”, букв. «попадания») в сравниваемых последовательностях; альтернативный алгоритм global alignment сравнивает последовательности как целое). NCBI BLAST использует local alignment, потому что он быстрее.
Hit – a short high scoring similarity region homologous sequences are likely to contain. Короткие фрагменты, совпадающие между сравниваемыми последовательностями; то, что ищет local alignment для того, чтобы «правильно» ориентировать последовательности друг относительно друга; размеры фрагментов и допустимая степень различий между последовательностями регулируются настройками программы.
(biological) Homology - similarity attributed to descent from a common ancestor. Гомология, сходство между биологическими объектами по нуклеотидным последовательностям или другим признакам, обусловленное происхождением от общего предка.
Sequence match – буквально, совпадение. В контексте: best match, closest match – последовательность, наиболее сходная с исследуемой (query).
Sequence identity - the extent to which two (nucleotide or amino acid) sequences have the same residues at the same positions in an alignment, often expressed as a percentage. Сходство выравненных последовательностей белков или нуклеиновых кислот, выраженное в процентной доле совпадающих амнокислот или нуклеотидов.
E-value - the Expect value (E) is a parameter that describes the number of hits one can "expect" to see by chance when searching a database of a particular size. Статистика, характеризующая вероятность того, что выявленное при локальном выравнивании сходство между query и subject последовательностями случайно, с учетом объема базы, с которой сравнивается query. Чем меньше Е, тем надежнее результат выравнивания. Если последовательности идентичны, E=0.
Query cover (or coverage) – what percentage of the search sequence overlaps with the aligned segments «Охват, перекрытие» - на сколько процентов ваша исследуемая последовательность перекрывается по длине с той, с которой вы ее сравниваете.
Max(imum) score (=bit-score) and Total score - a score is a numerical value that describes the overall quality of an alignment. The higher the score, the better the alignment. Max score is calculated from the part of the subject sequence that aligns best to the query. Total score is calculated from all hits. Индексы качества выравнивания, чем выше значения индексов, тем больше похожи друг на друга сиквенсы. На практике, если Total score > Max score, значит несколько разных участков subject гомологичны одному и тому же участку query. В graphic summary blast цветом показано распределение Max score, красным – лучше всего совпадающие последовательности или их участки.
Sequence Accession number - the most general identifier used in the NCBI sequence databases. This is the identifier that should be used when citing a database record in a publication. The specific version of a record is also tracked by another identifier that is mainly for internal NCBI use called the GI number. Уникальный идентификатор последовательности, предназначенный для публичного использования. На этот идентификатор следует ссылаться в публикациях. Параллельный “GI” идентификатор, начинающийся на GI предназначен для внутреннего (NCBI) использования. Для нуклеотидных последовательностей, как правило, структура кода - 2+6 - две буквы и шесть цифр, например JX669269; когда цифры кончатся NCBI перейдет на формат 2+8.
FASTA – An early sequence similarity (local alignment) search tool. The term FASTA is also used to identify a text format for sequences that is widely used. Each sequence in the file is identified by a single line title preceded by the greater than sign (">"). «Фаста», на практике, простейший, общеупотребимый формат сохранения последовательностей в текстовом файле; индивидуальные последовательности идентифицируются по заглавным строчкам со значком > в начале, и чем хочешь дальше.
DNA barcoding (Баркодиирование ДНК, ДНК-штрихкодирование, генетический баркодинг) — метод молекулярной идентификации, который позволяет по коротким последовательностям ДНК определять принадлежность организма к определённому таксону. В отличие от методов молекулярной филогенетики, ДНК-баркодирование используется для определения места данного организма в уже существующей классификации, а не для построения филогенетических деревьев или их дополнения. Наиболее часто используемым «баркодинговым» признаком для животных является участок митохондриального гена цитохромоксидазы I — MT-COI — из примерно 600 пар нуклеотидов. Для баркодинговых локусов разработаны «универсальные праймеры», т.е. праймеры применимые для широкого круга объектов, например «метазойные» — для группы Metazoa — «Фолмеровские» праймеры LCO1490: 5'-ggtcaacaaatcataaagatattgg-3' и HC02198: 5'-taaacttcagggtgaccaaaaaatca-3' (Folmer et al. 1994), которые "налипнут" на митохондриальную ДНК любого животного.