Monthly Archives: October 2011

Выпуск новостей 22.10.2011 (#3)

Публикуем очередной выпуск статей, в которых мы продолжаем рассказывать о методах и программах для реконструкции филогенетических деревьев.  Кроме того, возращаясь к теме выравнивания генетических последовательностей, мы рассказываем о программе Gblocks.

В разделе скрипты два обновления: наиболее интересным, с нашей точки зрения является скрипт на (bio)python для парсинга данных через интерфейсы BLAST/ENTREZ с конвертацией в формат фаста. Надеемся, что он пригодится не только новичкам в области филогенетики и эволюционной биологии, но и опытным пользователям – в среднем, он позволяет сократить время нужное для поиска и компоновки данных из Genbank на 95%.  Разумеется, мы будем очень рады услышать замечания и сообщения о багах.

Раздел “Практические руководства”

#Реконструкция филогенетического дерева с использованием PhyML – часть 1
#Реконструкция филогенетического дерева с использованием PhyML – часть 2
#Визуализация филогенетических деревьев и оценка качества их реконструкции
#Использование Gblocks для корректировки генетических данных после выравнивания – часть 1

Раздел “Скрипты”

#Скрипт для парсинга данных через BLAST/Entrez с автоматической конвертацией в формат .FASTA
#Biopython для парсинга данных из GenBank – интерфейсы BLAST/Entrez 

Использование Gblocks для корректировки генетических данных после выравнивания – часть 1

Как уже писалось ранее (здесь и здесь) выравнивание генетических последовательностей является нетривиальной задачей, требующей серьезных усилий особенно при наличии большого числа таксонов. В связи с этим, любая автоматизация и объективность при корректировке выровненных последовательностей позволяет существенно сократить время обработки данных. В этом упражнении мы рассмотрим основные параметры программы Gblock, позволяющей автоматизировать процесс корректировки генетических данных после выравнивания. 

Gblocks – программа, с помощью которой можно проводить выборку консервативных (conserved) блоков генетического кода из уже выровненных последовательностей (нуклеотидных, аминокислотных). Алгоритм программы работает следующим образом: на первом шаге уровень консервативности каждой позиции генетической последовательности оценивается и классифицируется на три группы: неконсервативные, консервативные и сильно консерватиные (nonconserved, conserved, highly conserved). Затем, блоки последовательных неконсервативных позиций длиннее определенного значения исключаются (иным языком, неконсервативные блоки удаляются из рассмотрения). В оставшихся блоках боковые позиции оцениваются на предмет консевартивности и удаляются до тех пор, пока с обеих сторон блок не будет окружен сильно консервативными позициями. Таким образом, оставшиеся блоки гарантированно являются корректно выровненными и представляют собой (с большой степенью вероятности) гомологи.

Результат обработки аминокислотных последовательностей программой Gblocks

Затем, удаляются все гэпы (пробелы, пустые позиции, gaps). При этом в данной программе гэпы могут быть определены тремя различными способами. После этого, некосервативные позиции примыкающие к пробелам также удаляются до тех пор, пока не достигнута консервативная позиция. В завершении, все блоки меньше определенной длины также удаляются из рассмотрения.

Важными настройками в главном меню программы являются следующие:

t. тип входных данных (протеины, кодоны или ДНК – Proteins, Codos, DNA).

o. позволяет указать входной файл (NBRF/PIR или FASTA). По умолчанию, при выборе файла в формате .fasta тип входных данных  автоматически становится Protein (протеин).

b. открывает меню настроек алгоритма Gblocks.

s. открывет меню настроек выходных файлов.

g. запускает программу.

q. выходит из программы.

Важные настройки в меню b (алгоритм Gblocks):

  1. Устанавливает минимальное число последовательностей для определения консервативных позиций. Это значение должно быть больше половины числа последовательностей во входном файле (например, если во входном файле присутствуют 11 последовательностей, то значение будет равно 6). Большие значения этого параметра уменьшают число выбранных позиций
  2. Устанавливает минимальное число последовательностей для определения боковых позиций. Это значение также должно быть равно или больше параметру 1 данного меню (минимальное число последовательностей для определения консервативных позиций). Большие значения этого параметра уменьшают число выбранных позиций
  3. Устанавливает максимальное число последовательных неконсервативных позиций. Все сегменты со значением большим чем заданное число последовательных неконсервативных позиций будут удалены. Большие значения этого параметра увеливают число выбранных позиций
  4. Устанавливает минимальную длину блока после удаления гэпов (пробелов).  Блоки меньше чем указанное значения будут удалены. Большие значения этого параметра уменьшают число выбранных позиций
  5. Переключение между тремя опциями для поиска гэпов:
    • None: в результирующем файле не должно быть пробелов. Все позиции, содержащие хотя бы один пробел рассматриваются как пробел и будут удалены
    • With Half: только позиции, где 50% и более процентов последовательностей содержит пробел рассматриваются как пробел.
    • All: все позиции с пробелами сохраняются в результирующем файле.
  6. r – восстановление исходных параметров меню.

В следующем упражнении мы рассмотрим возможности использования данной программы для коррекции выровненных генетических данных на примере, как обычно, ящериц рода Ctenotus.

Визуализация филогенетических деревьев и оценка качества их реконструкции

В одном из последних упражнений мы реконструировали филогенетическое дерево для рода ящериц Ctenotus с помощью программы PhyML.  Чтобы ускорить расчеты, реконструкция была продублирована на внешнем сервере. Обработка на домашнем компьютере заняла 45 часов, в сравнении с 7,5 часам на внешнем сервере. Таким образом, еще раз обращаем внимание пользователей на то, что для реконструкции филогенетических деревьев (современными методам и с современными объемами данных) требуются мощные сервера. 

Итак, проанализируем полученное филогенетическое дерево. Прежде всего нас интересуют следующие параметры:

  1. какую поддержку имеют внутренние узлы дерева?
  2. как именно распределены значения поддержки по узлам дерева?
  3. монофилетична ли исследуемая группа (клада) по отношению к внешней группе (outgroup)?
  4. дополнительно: сохранены ли таксономически монофилетичные группы?

Для первой (быстрой) визуалиции воспользуемся встроенным плагином на сайте сервера. Полученное дерево выглядит следующим образом:

Филогенетическое дерево рода Ctenotus

Нас интересуют значения обозначенные красным цветом около узлов филогенетического дерева. Нами было проведено 100 бутстрэпов, соответственно максимально возможное значение поддержки узла составляет 100. В среднем, процент поддерки узла от 80% и выше от числа бутстрэпов считается хорошим, от 70%-80% удовлетворительным, ниже 70% – слабым.

В нашем случае мы видим, что всего лишь 7 узлов имеют поддержку более 80%, и 1 узел – удовлетворительную, остальные узлы имеют слабое разрешение. Слабое разрешение в филогенетическом смысле обозначает то, сколько раз при бутстрэппинге была встречена данная конфигурация таксонов на узле. Например, 100% поддержка узла, ведущего к нашей внешней группе (виды Prasinohaema virens и Sphenomorphus muelleri), означает, что во всех бустрэпах, эти два вида были монофилетичны по отношению к остальным видам. Таким образом, мы находим ответ на третий вопрос – как расположена внешняя группа по отношению к нашей кладе.

Вернемся теперь к вопросу о том, как именно распределены значения поддержки по дереву. Наиболее “опасными” являются т.н. “глубокие” узлые (то есть расположенные ближе к корню дерева). Малые значения поддержки глубокого узла показывают, что большая часть реконструированных деревьев содержит разный набор видов в кладе, берущей начало от данного узла. Это значит, что число замен в генетических последовательностях использованных для реконструкции дерева недостаточно, чтобы разрешить (resolve) филогенетические отношения в группе. К сожалению, единственный метод улучшения такой ситуации – добавление генетических данных с большим числом полиморфных позиций (polymorphic site). 

Положительным моментом является то, что по отношению к внешней группе вся наша клада оказалась монофилетичной. 

Выводы:

  1. менее 20% узлов имеют высокую и удовлетворительную поддержку
  2. часть глубоких узлов имеет поддержку менее 50%, таким образом необходимо добавление генетических последовательностей: если их нет в Genbank (как можно проверить наличие данных в базе Genbank можно прочесть в этой статье), то сиквенирование
  3. клада Ctenotus монофилетична по отношению к внешней группе.

 В следующем упражнении мы разберемся с тем, какие файлы получаются в результате обработки данных программой PhyML и как их визуализировать с помощью программы FigTree.

Скрипт для парсинга данных через BLAST/Entrez с автоматической конвертацией в формат .FASTA

Одна из наиболее трудоемких задач встречающихся в биоинформатике, филогенетике и эволюционной биологии, как неоднократно писалось на этом сайте, это сбор и подготовка генетических последовательностей из специализированных баз данных. В связи с тем, что для некоторых задач поиск данных, оценка их качества, конвертация из одного формата в другой, получение определенным образом отформатированных заголовков, и многое другое становится рутинной и часто повторяющейся задачей, мной был написан скрипт, который позволяет минимизировать временные затраты.

Путем многочисленных проб и ошибок разработан следующий алгоритм для обработки данных:

  1. Поиск данных через BLASTn по референсной последовательности. Данный вариант является намного более точным, чем поиск через интерфейс NCBI по ключевым словам (например, название гена). К сожалению, база данных Genbank каждый день растет в размерах, и точность вносимой в описание данных информации остается на совести пользователей. Поэтому если выбирать между поиском через NCBI или BLASTn по последовательности, второй вариант позволяет найти на самом деле гомологичные данные. Условным недостатком такого метода является необходимость указать референсную последовательность, относительно которой будет проводится поиск в БЛАСТ. Чем более полно референсная последовательность покрывает искомый участок генома (ген, часть гена, несколько интронов-экзонов, и т.д.), тем выше шансы найти все гомологичные записи присутствующие на данный момент в базе.
  2. Проверка цепи ДНК (plus/minus strand). Одна из основных проблем ,возникающих при необходимости провести выравнивание последовательностей, это определение последовательностей сиквенированных с разных цепей ДНК (в английском варианте – plus/minus strand). В связи с этим, данный скрипт автоматически определеяет расположение последовательности и при необходимости конвертирует в нужную форму (делает reverse compliment).
  3. Проверка на наличие повторяющихся таксонов и автоматический выбор наиболее длинных последовательностей для каждого уникального таксона. Для построение филогенетических деревьев на межвидовом уровне пользователям часто приходится сталкиваться с необходимостью удаления повторяющихся таксонов. Наиболее логичным при этом является выбор наиболее длинных последовательностей из имеющихся
  4. Конвертация в формат FASTA. Здесь основным ньюансом является неоходимость выделения в заголовок только определенной информации (например, название таксона или порядковый номер записи в базе), и замена пробелов на знак подчеркивания (так как львиная доля программ для выравнивания последовательностей негативно относится к пробелам в заголовке)

Таким образом, разработанный алгоритм принимает на входе три настройки:

  • номер референсной записи из базы (уникальный ID из Genbank)
  • название гена (эта информация условна, и используется только для того, чтобы просвоить уникальный идентификатор результатам поиска)
  • название группы в пределах которой будет идти поиска гомологичных последовательностей (можно оставить NULL, если интересует любой гомолог, но эта опция не рекомендуется к использованию, за исключением случаев когда человек четко представляет возможные результаты своего шага -> размер выходного файла и блокировку аккаунта. Для любителей экспериментов было установлено максимальное число последовательностей в выходном файле – 1000)

Установленный настройки:  поиск через BLASTn/nr, e-value = 0.00001; максимальная длина последовательности 5000 bp; в fasta файле будут уникальные записи по названию таксона. Для запуска программы требуется установленный python/biopyhon. Пример записка скрипта:

python blast_2_fasta.py AB093081 matK Fagopyrum

Поиск будет осуществлен для гена MaturaseK (по референской записи AB093081, а не по названию гена) в пределах рода Fagopyrum

В следующей статье, будет выложен аналогичный код, но без фильтра по названию таксона (то есть сохранение всех уникальных записетей из базы).

Одна из наиболее трудоемких задач встречающихся в биоинформатике, филогенетике и эволюционной биологии, как неоднократно писалось на этом сайте, это сбор и подготовка генетических последовательностей из специализированных баз данных. В связи с тем, что для некоторых задач поиск данных, оценка их качества, конвертация из одного формата в другой, получение определенным образом отформатированных заголовков, и многое другое становится рутинной и часто повторяющейся задачей, мной был написан скрипт, который позволяет минимизировать временные затраты.

Путем многочисленных проб и ошибок разработан следующий алгоритм для обработки данных:

  1. Поиск данных через BLASTn по референсной последовательности. Данный вариант является намного более точным, чем поиск через интерфейс NCBI по ключевым словам (например, название гена). К сожалению, база данных Genbank каждый день растет в размерах, и точность вносимой в описание данных информации остается на совести пользователей. Поэтому если выбирать между поиском через NCBI или BLASTn по последовательности, второй вариант позволяет найти на самом деле гомологичные данные. Условным недостатком такого метода является необходимость указать референсную последовательность, относительно которой будет проводится поиск в БЛАСТ. Чем более полно референсная последовательность покрывает искомый участок генома (ген, часть гена, несколько интронов-экзонов, и т.д.), тем выше шансы найти все гомологичные записи присутствующие на данный момент в базе.
  2. Проверка цепи ДНК (plus/minus strand). Одна из основных проблем ,возникающих при необходимости провести выравнивание последовательностей, это определение последовательностей сиквенированных с разных цепей ДНК (в английском варианте – plus/minus strand). В связи с этим, данный скрипт автоматически определеяет расположение последовательности и при необходимости конвертирует в нужную форму (делает reverse compliment).
  3. Проверка на наличие повторяющихся таксонов и автоматический выбор наиболее длинных последовательностей для каждого уникального таксона. Для построение филогенетических деревьев на межвидовом уровне пользователям часто приходится сталкиваться с необходимостью удаления повторяющихся таксонов. Наиболее логичным при этом является выбор наиболее длинных последовательностей из имеющихся
  4. Конвертация в формат FASTA. Здесь основным ньюансом является неоходимость выделения в заголовок только определенной информации (например, название таксона или порядковый номер записи в базе), и замена пробелов на знак подчеркивания (так как львиная доля программ для выравнивания последовательностей негативно относится к пробелам в заголовке)

Таким образом, разработанный алгоритм принимает на входе три настройки:

  • номер референсной записи из базы (уникальный ID из Genbank)
  • название гена (эта информация условна, и используется только для того, чтобы просвоить уникальный идентификатор результатам поиска)
  • название группы в пределах которой будет идти поиска гомологичных последовательностей (можно оставить NULL, если интересует любой гомолог, но эта опция не рекомендуется к использованию, за исключением случаев когда человек четко представляет возможные результаты своего шага -> размер выходного файла и блокировку аккаунта. Для любителей экспериментов было установлено максимальное число последовательностей в выходном файле – 1000)

Установленный настройки:  поиск через BLASTn/nr, e-value = 0.00001; максимальная длина последовательности 5000 bp; в fasta файле будут уникальные записи по названию таксона. Для запуска программы требуется установленный python/biopyhon. Пример записка скрипта:

python blast_2_fasta.py AB093081 matK Fagopyrum

Поиск будет осуществлен для гена MaturaseK (по референской записи AB093081, а не по названию гена) в пределах рода Fagopyrum

В следующей статье, будет выложен аналогичный код, но без фильтра по названию таксона (то есть сохранение всех уникальных записетей из базы).
Код:

# ---------------------------------------------------------------------------
# blast_2_fasta.py
# Searching genbank data with BLAST/ENTREZ interfaces and converting it to a fasta files.
# Features:
# * check strands and do reverse compliment automatically
# * records the longest sequences for a given gene/taxa name
# * allow flexible naming in resulting fasta file (only taxa name, taxa + id, only id, etc)
# * automatically detects sequences that are too long (e.g. full genome scans) and excludes them (by default 5000 bp)
# * has a flexible way to change BLAST searching settings (e.g. e-value, maximum number of returned sequences)
# * write additionally information files with exlcuded sequences
# Author: Anna Kostikova (anna.kostikova[]gmail.com)
# Date: 19.10.2011
# ---------------------------------------------------------------------------

#parse BLAST results and save as genbank file, then convert into fasta file, keeping only non repeating,longest records (non-repeating by organism name), max sequence threshold = 5000 by default
import os
import sys
import re
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML
from Bio import Entrez
from Bio import SeqIO
Entrez.email = ‘youremail@domain.com’

#Do blastn search and search results in xml file
def Qblast_search_2_xml(seq, gene_name, organism, e_val=0.00001, max_aligns=1000, format = “XML”):
print “Starting QBlast search…”
if organism == ‘NULL’:#this setting allow to NOT specify an organism. Search will be restricted to a max of 1000 records
result_handle = NCBIWWW.qblast(“blastn”, “nr”, sequence=seq, expect=e_val, hitlist_size=1000, alignments=max_aligns, format_type=format) #setting for the blastn search
else:
org = organism + “[orgn]” #set restiction to a particular taxon (e.g. genus, family, order, etc.)
result_handle = NCBIWWW.qblast(“blastn”, “nr”, sequence=seq, expect=e_val, entrez_query=org, hitlist_size=1000, alignments=max_aligns, format_type=format) #setting for the blastn search
blast_results = result_handle.read() # handle results
save_file = open(organism+”_”+gene_name+”.xml”, “w”)#name for an xml file
print “Writing QBlast search results…”
save_file.write(blast_results) #write to xml file and save
save_file.close() #close xml
print “\nFinished QBlast.”

#Parse XML from BLAST to creat list of all downloaded accessions
def Blast_xml_accessions(organism, gene_name):
result_handle = open(organism+”_”+gene_name+”.xml”) #Open blast xml and load the blast record
blast_records = NCBIXML.parse(result_handle) #do parsing of xml file
blast_record = blast_records.next() # for each record
output = {} #Store all accession ids
for x in blast_record.alignments:
output[x.accession] = [x.length]
print “Total ” + str(len(output)) + ” records retrieved.”
return output

#Download .gb file based on accession list provided
def Download_gb_files(acc_list,organism, gene_name):
local_file=open(organism+”_”+gene_name+”.gb”, ‘w’)#Do entrez.efetch to retrieve all records in genbank file, store locally
counter = 1
for x in acc_list:#for each record in accessions list, search NCBI
print “Downloading genbank records ” + str(counter) + “…”
try:
handle = Entrez.efetch(db=”nucleotide”,id=x,rettype=”gb”) #download record with Entrez.efetch
local_file.write(handle.read()) #write to a file
handle.close()
except:
print “Accsession id is not found”
counter = counter + 1
local_file.close()
print “\nFinished genbank processing”

#Read-in genbank file and make a disctionary of taxa names (hence, it will be unique names)
def Genbank_find_longest(gb_filename, max_len_threshold=5000):
input_handle = open(gb_filename, “r”)
mylist = []
for seq_record in SeqIO.parse(input_handle, “genbank”) :#parse genbank file and create 2 dictionaries
#print “Dealing with GenBank record %s” % seq_record.id
data_len = len(seq_record.seq)
mylist.append(seq_record.annotations[“organism”])

keys=dict.fromkeys(mylist).keys() #dictionary species(key):sequence_length(value) pair. For now value is empty
values=[0] * len(keys)
dic_length = dict(zip(keys, values))
print “\nYour dictionary is %s length (%s species in *.gb file)\n” % (len(dic_length),len(mylist))
dic_access = dict(zip(keys, values)) #dictionary species(key):accession_id(value) pair. For now value is empty
input_handle.close()

input_handle = open(gb_filename, “r”)
for seq_record in SeqIO.parse(input_handle, “genbank”) :#fill-in dictionaries
#print “Dealing with GenBank record %s” % seq_record.id
sp = seq_record.annotations[“organism”]
if (not dic_length.get(sp) is None and len(seq_record.seq) < max_len_threshold): #if species from file is present in dictionary and its length is < than max threshold
if dic_length.get(sp) dic_length[sp] = len(seq_record.seq)
dic_access[sp] = seq_record.id
else:
if (dic_length.get(sp) < len(seq_record.seq)): #if longer than the seq already in a dictionary – replace dic_length[sp] = len(seq_record.seq) dic_access[sp] = seq_record.id #else: # print “Record %s for species %s is longer than %s, skipping…” % (dic_access[sp],sp,seq_record.id) else: print “%s (uid %s) is %s bp longer than sequence threshold (%s bp), skipped. To increase threshold, change default in Genbank_find_longest function” % (sp,seq_record.id,(len(seq_record.seq)-max_len_threshold),max_len_threshold) input_handle.close() print “\nFinshed processing dictionaries. Start retrieval of the longest sequences…” return(dic_access) #Convert genbank file to fasta file. In addition make 2 files: 1) species+accessions 2) skipped species #It also checks for strand position and makes a reverse compliment if necessary def Genbank_2_fasta(dictionary, gb_filename): fas_filename = gb_filename[:-3] + “.fasta” dic_filename = gb_filename[:-3] + “.txt” skip_filename = gb_filename[:-3] + “_skip.txt” input_handle = open(gb_filename, “r”) output_handle = open(fas_filename, “w”) output_handle_dic = open(dic_filename, “w”) output_handle_skip = open(skip_filename, “w”) for seq_record in SeqIO.parse(input_handle, “genbank”) : #print “Dealing with GenBank record %s” % seq_record.id uid = seq_record.id if uid in dictionary.values(): #print seq_record.features[0].strand if seq_record.features[0].strand == 1: output_handle.write(“>%s\n%s\n” % (
re.sub(“\s”, “_”, seq_record.annotations[“organism”]), #this way we replace space with an underscore
seq_record.seq))
output_handle_dic.write(“%s,%s\n” % (
re.sub(“\s”, “_”, seq_record.annotations[“organism”]),
seq_record.id))
elif seq_record.features[0].strand == 2: #check for reverse compliment
output_handle.write(“>%s_%s\n%s\n” % (
re.sub(“\s”, “_”, seq_record.annotations[“organism”]),
seq_record.seq.reverse_complement()))
output_handle_dic.write(“%s,%s\n” % (
re.sub(“\s”, “_”, seq_record.annotations[“organism”]),
seq_record.id))
else:
print “Sequences %r has an unknown strand, skipping” % uid
else:
#print “This %s sequences is not the longest for this gene. Skipping.” % uid
output_handle_skip.write(“%s,%s,%s,%s\n” % (seq_record.id,seq_record.annotations[“organism”],len(seq_record.seq),uid))
print “\nGenbank to fasta converstion is successfully done.\nCheck your folder for %s.gb, %s.fasta, %s.txt, %s_skip.txt files” % (gb_filename[:-3],gb_filename[:-3],gb_filename[:-3],gb_filename[:-3])
input_handle.close()
output_handle.close()
output_handle_dic.close()
output_handle_skip.close()

################## end of functions ######################

#INPUT: reference accession id, gene name, organism name

def usage():
print “Usage: blast_2_genbank.py sequence_id gene_name organism”
print “For instance (windows OS): python blast_to_genbank_def.py FJ872097 matK Rheum”
print “For any taxa, set 3rd argument to NULL (default): python blast_to_genbank_def.py FJ872097 matK NULL”
sys.exit( 0 )

if __name__ == ‘__main__’:
args = sys.argv[ 1: ]

if len( args ) != 3:
usage()

#reading arguments from the file
seq1 = os.path.normpath( args[ 0 ] )
gene_name1 = os.path.normpath( args[ 1 ] )
organism1 = args[ 2 ]

#call to functions
#example: python blast_to_genbank_def.py AB093081 matK Fagopyrum

Qblast_search_2_xml(seq1,gene_name1,organism1)
list_of_accessions = Blast_xml_accessions(organism1,gene_name1)
Download_gb_files(list_of_accessions,organism1,gene_name1)
outputfile = organism1+”_”+gene_name1+”.gb”
my_dic = Genbank_find_longest(outputfile)
Genbank_2_fasta(my_dic, outputfile)

Объявление: загрузка данных в архивах с сайта

Небольшое объявление для читателей сайта, которые хотят не просто прочитать все упражнения, но еще и проверить работают ли они:) (то есть проверить на тех данных, которые используются в упражнениях).

По (пока) неизвестным для разработчиков сайта причинам, скачка файлов с нашего фтп-сервера требует обязательного логина и пароля. Так как нам пока не удалось выяснить, каким образом сделать открытый доступ всем желающим, мы предлагаем всем заинтересованным оставить сообщение к данному посту с адресом электронной почты (это сообщение будет видно только администрации сайта). Мы, в ответ, с удовольствием предоставим полный доступ к фтп-сереверу (к папкам с данными).

В то же время, будем надеяться, что в скором времени нам удасться разрешить неприятную деталь с фтп-сервером и скачивать данные к упражнениям можно будет простым нажатием на ссылку в статье.

Реконструкция филогенетического дерева с использованием PhyML – часть 2

В предыдущем упражнении мы разобрались с основными настройками доступными в программе PhyML. В этом упражнении мы проведем реконструкцию филогенетического дерева для рода ящериц Ctenotus. Наш набор данных состоит из 4 участков генома (ATP, NADH, 16S и 12S) общей длинной 2955 позиций и 45 видов. Для получения более детальной информации о том, как именно создать такой набор данных в формате phylip можно обратиться к упражнениям здесь и здесь.

Для оценки топологии филогенетического дерева необходимо загрузить дистрибутив программы phуml в скомпилированном или нескомпилированном виде (для *nix систем).  Далее, нужно скачать файл в формате PHYLIP с выровненными нуклеотидными последовательностями с нашего ftp сервера (потребуется пароль и логин).

Двойным щелчком запускаем программу PhyML и в первом подменю вводим название файла ctenotus_4genes.phy и нажимаем enter. Перед нами появится первая часть настроек. Переход к следующей странице настроек осуществляется с помощью кнопок + и –. Нажатие кнопки Y и затем enter приводит к запуску анализа. Для выбора настроек каждого подменю (задается заглавными буквами – D, I, M и т.д.) необходимо несколько раз нажать эту букву на клавиатуре – таким образом вы сможете прокручивать возможные варианты подменю.

Выбор файла с нуклеотидными или аминокислотными последовательностями

Во втором подменю выбраны следующие настройки:

D – тип данных DNA – оставляем без изменений, так как наши последовательности – это нуклеотидные последовательности.

I – формат данных – формат наших данных, полученных с помощью программы SequenceMatrix (смотреть упражнение) при объединении четырех участков генома (ATP, NADH, 16S и 12S), interleaved, поэтому эту опцию мы также оставляем без изменений

M и R – мы будем проводить реконструкцию только для одного набора данных, поэтому эти две опции мы также не меняем

Выбор типа и формата данных

Для перехода в следующее подменю нажимаем  + и enter.  В третьем подменю указываем следующие настройки:

M – выбор модели нуклеотидных замен. Для этих данных, мы выбираем модель  GTR (для этого необходимо 3 раза нажать на букву M на клавиатуре)

F – мы включаем опцию оценки частот нуклеотидов нажав один раз на F

V – мы не включаем опцию (то есть оставляем fixed, как задано по умолчанию)

С – число категорий можно оставить равным 4 или увеличить до 16 (в зависимости от мощности вашего компьютера)

A – мы оставляем эту опцию включенной, таким образом наша результирующая модель замены нуклеотидов выглядит как GTR+Gamma

G – мы меняем эту опцию с mean на median (cо среднего на медиану) один раз нажав на G

Выбор модели нуклеотидных замен

Нажимаем +, enter и переходим в четвертое подменю:

В данном подменю единственный параметр который мы меняем это S – мы задаем поиск топологии через Best of NNI и SPR (дважды нажать на S). Все остальные параметры оставляем неизменными

Нажимаем +, enter и переходим в пятое подменю:

В пятом подменю мы включаем непараметрический бутстрэппинг – один раз нажав на B, после чего вводим число повторов 100, и затем указываем Y при ответе на вопрос Print Bootstrap trees. Эта последняя настройка, которую нужно поменять. Для запуска программы нажимаем Y. В зависимости от мощности компьтера оценка топологии дерева может занять до 10 часов.

Бутстрэппинг для оценки поддержки узлов филогенетического дерева

В следующем упражнении мы визуализируем полученное дерево и разберемся с тем, как оценивать качество реконструкции.

Реконструкция филогенетического дерева с использованием PhyML – часть 1

PhyML – одна из наиболее популярных программ для реконструкции филогенетических деревьев методом максимального правдоподобия. Программа работает как с нуклеотидными, так и с аминокислотными последовательностями. Ключевым преимуществом программы, по сравнению с другими похожими приложениями, является большой спектр моделей нуклеотидных замен, а также широкие возможности исследования топологического пространства филогенетических деревьев (то есть вариантов реконструкций). Кроме того, в программе реализованы два метода по оценке поддержки узлов дерева – непараметрический бутстрэппинг (non-parametric bootstrap) и приблизительный тест отношения правдоподобий (approximate likelihoodbranch supports).

PhyML был разработан для работы с наборами данных среднего и большого размера: теоретически возможна работа с числом нуклеотидных последовательностей до 4000 и максимальной длинной в 2 000 000 позиций.Тем не менее, вычислительные емкости компьютера, на котором осуществляется поиск и оценка топологии дерева, должны быть достаточно большими. На практике, основное ограничение – это число нуклеотидных последовательностей и рекомендуемая цифра в данный момент (2011) лежит в пределах 500 последовательностей с длинной до 10000 позиций. Входной формат данных программы – Phylip (interleaved или sequential)

В данном упражнении мы кратко познакомимся с различными опциями доступными в программе PhyML, а в следующем упражнении проведем реконструкцию дерева для ящериц рода Ctenotus.

Подменю 1

В первом подменю необходимо указать название входного файла с нуклеотидными или аминокислотными, выровнеными последовательностями. Переход к следующей странице (подменю) настроек осуществляется с помощью кнопок + и –. Нажатие кнопки Y и затем enter приводит к запуску анализа. Для выбора настроек каждого подменю (опции подменю задаются заглавными буквами – D, I, M и т.д.) необходимо несколько раз нажать эту букву на клавиатуре – таким образом вы сможете прокручивать возможные варианты опций.

Подменю 2

D – тип данных – нуклеотидные (DNA) или аминокислотные (AA) последовательности

I – входной формат данных: interleaved – разбитый на блоки равной длины (Картинка 1) или sequential – где каждая запись представлена непрерывной строкой (Картинка 2)

Interleaved format

Sequential format

M – будет ли проводиться анализ несколькоих набров данных. При выборе этой опции необходимо указать число наборов данных во входном файле. Например, если в одном файле у вас присутствуют выровненные последовательности для двух родов растений, то можно последовательно провести реконструкции обоих наборов данных.

R – присвоение уникального индекса результирующим файлам. Например, если вы хотите провести анализ набора данных для трех моделей нуклеотидных замен, то эта опция позволит присвоить уникальное имя каждому анализу.

Подменю  3


M – модель замены нуклеотидов (доступные значения: JC60, K80, F81, F84, HKY85, TN93, GTR и пользовательская  – для нуклеотидных последовательностей; LG, WAG, Dayhoff, JTT, Blosum62, mtREV, rtREV, cpREV, DCMut, VT, mtMAM и пользовательская  – для аминокислотных последовательностей). Особенности различных моделей замены нуклеотидов и моделей аминокислотных последовательностей будут рассмотрены в отдельной статье.

F – возможность оценки частот нуклеотидов из данных в ходе реконструкции филогенетического дерева. Когда выбрана пользовательская модель замены нуклеотидов, данная опция позволяет  задать распределение частот нуклеотидов в эквилибриуме. Для аминокислотных последовательностей, стационарные частоты определяются либо моделью замен, либо оцениваются из исходного набора данных. Таким образом, значение опции F варьируется в зависимости от того, какие данные находятся в обработке.

T – возможность оценки отношения трансверсий/транзиций. Эта опция доступна только если выбраны модели K80, HKY85 или TN93. Отношение транзиций к трансверсиям это отношение числа транзиций к числу трансверсий для двух последовательностей ДНК.

V – пропорция позиций без нуклеотидных замен, то есть ожидаемая частота позиций без изменений. Параметр может либо иметь фиксированное значение, либо может быть оценен в ходе реконструкции филогенетического дерева. По умолчанию число зафиксировано на 0. Таким образом, считается, что даже если результирующие последовательности не различаются между собой, каждая позиция последовательности ДНК может, тем не менее, аккумулировать замены.

R – оценивать ли скорость замен отдельно для всех (некоторых) категория замещений

С – число категорий для скорости замен (substitution rates)

А – оценка парамтра Гамма-распределения

G – среднее или медиана для каждой категории скорости замен

Скорость мутационного процесса может варьироваться от одной нуклеотидной позиции в последовательности к другой. Эта неоднородность в скоростях мутационного процесса может быть аппроксимирована с помощью дискретного гамма-распределения. Выключение этой опции осуществляется нажатием на R. При включении опции R (то есть при оценке гетерогенности скорости мутаций между различными позициями нуклеотидной последовательности) по умолчанию число категорий дискретного распределения равно 4. Меньшее число категорий указывать не рекомендуется. Увеличие числа категорий ведет к более точной оценке скорости нуклеотидных замен, однако приводит к существенному замедлению процесса реконструкции дерева. Обычно задают число от 4 до 20. Изменение параметра возможно через подменю C. Подменю G позволяет выбрать, будет ли оцениваться среднее или медиана для каждой категории гамма-распределения.

Подменю 4


O – оптимизация топологии дерева. По умолчанию, топология дерева оптимизируется для максимизации значения правдоподобия. Однако, можно расчитать значение правдоподобия для фиксированной топологии дерева – это делается с помощью опции О.

U – когда опция оптимизаций топологии дерева включена (опция О), PhyML будет осуществлять поиск относительно некоторого исходного дерева. По умолчанию, таким деревом явялется дерево созданное по алгоризму BioNJ. Альтернативная опция –  топология методом “parsimony” или некоторое пользовательское дерево (в формате Newick). После этого, с помощью опции S, нужно задать то, каким образом будет осуществляться поиск более корректной топологии дерева.

S – в программе реализовано три метода для поиска топологии деревьев: NNI (nearest neighbor joining), SPR (subtree pruning and regrafting) и BEST (комбинация этих двух методов)

Подменю  5

B – поддержка внутренних узлов реконструированного дерева может быть оценена с помощью непараметрического бутстрэппинга. По умолчанию данная опция выключена, при включении необходимо указать число повторов для бутстрэппинга. Чем больше число повторов, тем более точным будет число поддержки узла. Однако, бутсрэппинг линейно увеличивает число нужного программе времени (так как каждый раз филогенетическое дерево оценивается по новой). В среднем, число бутстрэпов рекомендуется устанавливать равным 1000.  Минимальное число повторов – 100.

A – при выключенно бутстрэппинге, приблизительный тест отношения правдоподобий будет по умолчанию включен. Эта опция значительно быстрее бутсрэппинга, но предпочтительной является оценка качества реконструкции с помощью бутстрэппинга.

Biopython для парсинга данных из GenBank – интерфейсы BLAST/Entrez

Очень часто в работе с большими объемами генетических данных возникает необходимость загрузить и обработать нуклеотидные или аминокислотные последовательности в батч-режиме. В данном упражнении мы покажем как, с использованием небольшого скрипта на языке python,  можно осуществить поиск данных через интефейсы BLAST/Entrez.

Веб-интерфейс BLAST

Задача данного упражнения и скрипта: найти все записи в базе Genbank через интерфейс BLAST (нуклеотидный поиск) близкие к искомой, загрузить результат поиска в формате xml на локальный компьютер и затем, с помощью интерфейса Entrez загрузить записи в формате Genbank.

Для чего нам может потребоваться данная операция?

  • поиск большого числа различных генов через BLAST для определенных групп организмов
  • автоматизация отбора подходящих записей
  • загрузка файлов в формате генбанк для их последующей обработки

Для работы со скриптом требуется установленный Python, Biopython и NumPy.

Алгоритм скрипта устроен следующим образом:

  1. Пользователь должен указать два параметра:  номер записи из Генбанка, по которой будет осуществляться нуклеотидный поиск (то есть поиск подобных записей) и название организма (это может быть вид, род, семейство и т.д.)
  2. Скрипт с использованием класс NCBIWWW.qblast проведен поиск со следующими параметрами:
    blastn – обращение к интерфейсу blastn
    nr – использовать нуклеотидную коллекцию данных (nr/nt)
    acc_id – номер записи, по отношению к которой искать похожие последовательности
    entrez_query – дополнительные параметры поиска (ограничения), в нашем случае это будет только имя организма
    hitlist_size=1000 – максимальное число возвращаемых записей
    format_type=”XML” – формат результирующего файла
  3. В Интернет-версии интерфейса БЛАСТ это соответствовало бы:
    blastn – Database
    nr – Others
    acc_id – Enter accession number(s), gi(s), or FASTA sequence(s)
    entrez_query – Organism
    hitlist_size=1000 – Select the maximum number of aligned sequences to display
    format_type=”XML” – выбирается на странице с результатами поиска
  4. Результаты поиска будут сохранены в директории откуда запускается скрипт и будут иметь расширение XML и название по имени организма
  5. После загрузки файла в формате XML, мы должны обратить опять к базе, но теперь уже через интерфейс Entrez – это позволит нам скачать данные в формате GB. Это мы делаем с использованием строки
    Entrez.efetch(db=”nucleotide”,id=x,rettype=”gb”) – где прописываем, что нам нужны данные в формате генбанк, нуклеотидные последовательности.
  6. Результирующий файл также будет сохранен в директории откуда запущен файл и будет иметь разрешение gb

Скрипт:

# ---------------------------------------------------------------------------
# blast_to_genbank.py
# Loading genbank data with BLAST/ENTREZ interfaces
# Author: Anna Kostikova
# Date: 07.10.2011
# ---------------------------------------------------------------------------

#parse BLAST results and save as genbank file
import os
import sys
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML
from Bio import Entrez
from Bio import SeqIO
Entrez.email = 'youremail@server.com' #necessary for database admin to know that you are not a robot

#if len(sys.argv) != 2:
# stop the program and print an error message
#   sys.exit("Must provide accession id and genus/family name")

acc_id = "HM217109"
file = "rheum"
organism = file + "[orgn]"

#acc_id = str(sys.argv[1])
#file = str(sys.argv[2])
print 'Accesion_id:',acc_id
print 'File name:',file

#Do blastn search and search results in xml file
print "Starting QBlast search..."
result_handle = NCBIWWW.qblast("blastn", "nr", acc_id, entrez_query=organism, hitlist_size=1000, format_type="XML")
blast_results = result_handle.read()
save_file = open(file+".xml", "w")
print "Writing QBlast search results..."
save_file.write(blast_results)
save_file.close()
print "Finished QBlast."

#Open blast xml and load the blast record
result_handle = open(file+".xml")
blast_records = NCBIXML.parse(result_handle)
blast_record = blast_records.next()

#Store all accession ids
output = {}
for x in blast_record.alignments:
output[x.accession] = [x.length]
print "Total " + str(len(output)) + " records retrieved."

#Do entrez.efetch to retrieve all records in genbank file, store locally
local_file=open(file + ".gb", 'w')

counter = 1
for x in output:
print "Downloading genbank records " + str(counter) + "..."
try:
handle = Entrez.efetch(db="nucleotide",id=x,rettype="gb", )
local_file.write(handle.read())
handle.close()
except:
print "Accsession id is not found"
counter = counter + 1

local_file.close()
print "Finished genbank processing."