Kodomo

Пользователь

Учебная страница курса биоинформатики,
год поступления 2015

Ещё немного о BioPython

Сегодня мы разберем, как работать с файлами аннотаций и как изнутри программы скачивать данные из ncbi. Более подробно см. документацию BioPython:

http://biopython.org/DIST/docs/tutorial/Tutorial.html

Работа с файлами GenBank

Из предыдущего занятия мы знаем, что последовательности обрабатываются как объекты класса Seq. Для обработки аннотированных последовательностей есть класс SeqRecord, который кроме собственно последовательности умеет хранить информацию о разных биологических объектах, содержащихся в этой последовательности: генах, сайтах и пр. Аннотации часто хранятся в формате GenBank, и BioPython умеет с этим форматом работать. Допустим, у нас есть аннотированный геном бактерии Bacillus subtilis. Прочтем его. 

   1 from Bio import SeqIO
   2 record = SeqIO.read("sequence.gb", "genbank")

Теперь пройдемся по всем аннотированным участкам (feature) и посмотрим, каких типов они бывают. 

   1 features = {}
   2 for f in record.features:
   3     type = f.type
   4     features[type] = features.get(type,0)+1
   5 print features

У каждого объекта feature класса Feature есть его положение – feature.location, – в котором записаны его начало, конец и цепь, на которой он лежит. Кроме того, у каждой feature в feature.qualifiers лежит словарь, где ключами являются все те поля, которые были в файле. Например, у некоторого гена были поля: /gene="dnaN" /locus_tag="BSU00020" /db_xref="GeneID:939970". Тогда словарь feature.qualifiers будет выглядеть так: {"gene":"dnaN", "locus_tag":"BSU00020", "db_xref":"GeneID:939970"}

Доступ к NCBI из Biopython

Бывает удобно скачивать информацию из баз NCBI на лету, внутри кода программы. Плохая новость: много так не накачаешь – они трепетно отслеживают частоту и объемы скачиваний и банят за превышения. Хорошая новость: биопитон эти ограничения знает и частично сам отслеживает. Нам понадобится подпакет Entrez. При подключении необходимо представиться: на этот адрес ncbi будет слать гневные письма, прежде чем забанить. 

   1 from Bio import Entrez
   2 Entrez.email = "something@something.com"

1. Доступные базы данных.

Получим список доступных баз данных: 

   1 from Bio import Entrez
   2 Entrez.email = "something@something.com"
   3 handle = Entrez.einfo()
   4 record = Entrez.read(handle)
   5 print record["DbList"]

Здесь record – словарь с единственным параметром.

2. Подключение к БД и поиск в ней

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

   1 handle = Entrez.einfo(db="pubmed")
   2 record = Entrez.read(handle)
   3 print record["DbInfo"].keys()
   4 print record["DbInfo"]["Count"]
   5 print record["DbInfo"]["LastUpdate"]

Теперь найдем число статей про Last Universal Common Ancestor (LUCA) – нам потребуется функция esearch() 

   1 handle = Entrez.esearch(db="pubmed", term="LUCA")
   2 record = Entrez.read(handle)
   3 print record["IdList"]
   4 print record["Count"]

Параметр term – что, собственно, искать. По умолчанию Entrez записывает в IdList только 20 находок, это регулируется параметром RetMax. Достанем идентификаторы всех статей Кунина за 2014 год. 

   1 handle = Entrez.esearch(db="pubmed", term="Koonin[AUTH] AND 2014[MHDA]", RetMax=40) 
   2 record = Entrez.read(handle)
   3 print record["IdList"]
   4 print record["Count"]

Все возможные поля поиска можно посмотреть так: 

   1 handle = Entrez.einfo(db="pubmed")
   2 record = Entrez.read(handle)                    
   3 for field in record["DbInfo"]["FieldList"]:
   4     print("%(Name)s, %(FullName)s, %(Description)s" % field)

3. Скачивание из БД

Чтобы что-то скачать, понадобится функция efetch(). Список параметров можно посмотреть по ссылке. Давайте скачаем одну из статей про LUCA по её идентификатору. 

   1 handle = Entrez.efetch(db="pubmed", id="27886196", rettype="gb", retmode="text")
   2 record = handle.read()
   3 print record

Здесь скачается не вся статья, и даже не всегда абстракт – только базовая информация о публикации. Обратите внимание на разницу между handle.read() и Entrez.read(handle). Первое читает обычный текст, второе парсит XML. По умолчанию все ответы приходят в формате XML.

Теперь скачаем все последовательности гена rpl16 из опунций. 

   1 handle = Entrez.esearch(db="nuccore", term="Opuntia AND rpl16")
   2 record = Entrez.read(handle)
   3 print(record["Count"])
   4 gi_list = record["IdList"]
   5 handle = Entrez.efetch(db="nuccore", id=gi_list, rettype="gb", retmode="text")
   6 text = handle.read()
   7 print text

Полученный файл можно сразу скормить модулю Seq: 

   1 records = SeqIO.parse(handle, "gb")
   2 handle = Entrez.efetch(db="nuccore", id=gi_list, rettype="gb", retmode="text")
   3 records = SeqIO.parse(handle, "gb")
   4 for record in records:
   5     print("%s, length %i, with %i features" % (record.name, len(record), len(record.features)))

Запрос можно уточнять и прописывать отдельные поля: term="Opuntia[ORGN] AND rpl16[GENE]". Ну а как посмотреть все существующие поля с из описаниями, вы уже знаете.