Батч-загрузка данных из GENBANK

By | 5 June 2011

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

Наша цель сегодня – загрузить несколько различных участков ДНК кода для всех видов Ctenotus, доступных в Genbank. Информация о том, какие участки генома были сиквенированы и с каким успехом, мы получили из статьи D. Rabosky (2007) в PRSofB .

Для батч-загрузки данных из БД, мы воспользуемся скриптом Perl, который осуществляет поиск и загрузку всех записей, удовлятворяющих следующим поисковым условиям.

  1. Имена видов
  2. все записи должны принадлежать роду Ctenotus

  3. Название участков генома
    Доступны следующие участки генома ктенотусов (информация эта, подчеркиваю, может быть получена не только из статьи Дэна Рабоски, но и с помощью простого изучения базы GB) 

    • nuclear ATP synthetase b-subunit intron
    • mitochondrial ND4
    • tRNAs: tRNA-ser, tRNA-his, tRNA-leu
    • mitochondrial 12S rRNA
    • mitochondrial 16S rRNA
  4. Обратите внимание, что митохондриальная ND4 (NADH dehydrogenase subunit 4) идет в базе вместе с tRNA’азми, поэтому не имеет смысла обрабатывать такие последовательности отдельно друг от друга. В результате, у нас получаются 4 участка для загрузки: (1) ATP (2) ND4 + tRNA’ы (3) 12S и (4)  16S.

  5. Длина последовательности участка
  6. В случае с ящерицами, это не так важно, но в целом рекомендуется указывать границы протяженности участка ДНК (например от 300 до 5000 нуклеотидов). Но это условие подразумевает, что вы заранее знаете длину последовательности и можете поставить ограничение. Например, известно заранее, что авторы статьи сиквенировали 671 bp для ATP  – таким образом мы сразу исключаем слишком короткие последовательности, явным образом неудовлетворяющие нашей задаче. Однако если длина последовательности неизвестна – лучше это условие пропустить.

В результате, у нас образуются три поисковых правила (прочитать о структуре БД Genbank и номенклатуре и свойствах полей можно в практикуме 2):

“Ctenotus”[Organism]  AND ATP[All Fields] AND 300:20000[SLEN]

В таком запросе мы указываем, что:

  • поле Organism, обязательно должно содержать слово Ctenotus
  • другие поля должны содержать слово ATP (название интрона) хотя бы единожды
  • поле SLEN должно быть в пределах от 300 до 20000 нуклеотидов

Важно понимать, что условие на название интрона очень гибкое –  это сделано мной умышленно, потому что к сожалению, не все записи в БД Genbank имеют название сиквенированной последовательности в правильном поле.

Скрипт, которым мы воспользуемся, решает следующую задачу:

На входе, скрипт получает файл, где каждая строка представляет собой набор условий для поиска. На данный момент, закодированы только три поля: название организма, название гена, и длина гена. Кроме того, в скрипте реализован негибкий вариант поиска – когда по умолчанию запрос идет с логическим оператором AND между условиями. Например, если наш входной файл с условиями запроса выглядит следующим образом:

ctenotus,ATP,300:2000
ctenotus,NADH,300:20000
ctenotus,12S,300:20000
ctenotus,16S,300:20000

Это значит, что первая поисковая строка будет равнозначна следующему ручному запросу в Genbank:

ctenotus[ORGN] AND ATP[All Fields] AND 300:2000[SLEN]
Вторая:
ctenotus[ORGN] AND NADH[All Fields] AND 300:20000[SLEN]
и т.д.

Скрипт выглядит следующим образом:

#!/usr/bin/perl
use Bio::DB::GenBank;
use Bio::DB::Query::GenBank;
open (DRR, "conditions.txt");
my @alldir=<DRR>;
$num=@alldir;
print "number of searches: $num\n\n";
for ($inc = 0; $inc < $num; $inc++) {
#read line from file and split
$item = @alldir[$inc];
chomp $item; #clean line
($species,$gene,$length) = split(",", $item); #split by comma as separator
print "searching for: $species and $gene\n\n"; #print species and gene for search
#make up a query
$query_string=$species."[ORGN] AND ".$gene."[All Fields] AND ".$length."[SLEN]";
$file=$species."_".$gene.".gb"; #create output file name
print "Query: $query_string\n"; #print out query string
$query = Bio::DB::Query::GenBank->new(-db=>'nucleotide', -query=>$query_string); #form new genbank object
$count = $query->count; #get number of sequences for a query
print "Number of sequences found: $count\n"; #print it out
# get a genbank database handle
$gb = Bio::DB::GenBank->new; #(-retrievaltype=>'tempfile',-format=>'Genbank');
$stream = $gb->get_Stream_by_query($query); #returns a Bio::SeqIO stream object
$seqio=Bio::SeqIO->new(-format=>'Genbank',-file=>">$file");
$i=0;
#make a stream to write out sequences
 while ($seq = $stream->next_seq) {
$seqio->write_seq($seq);
$i++;
print "Writing entry $i out of $count\r";
}
$seqio->close();
if($count>0) {       print "\nDone.\n\n";    } #final check if something is found
else {       print "Done.\n";    } #if no sequences found
}
close (DRR); #close file
exit(); #exit

NB: Этот скрипт заслуживает внимания и аккуратной чистки (написано не в режиме strict, без exceptions, без обращения к пользователю за именем файла, и т.д.), но на данный момент мы выспользуемся этим неаккуратным вариантом.

Чтобы запустить скрипт нужно сделать несколько вещей (внимание – скрипт требует установленнго ActivePerl и библиотек Bioperl):

  1. Создаем рабочую директорию Work в удобном для вас месте
  2. В ней создаем новый текстовый файл, который называем conditions.txt и копируем туда наши поисковые условия:
  3. ctenotus,ATP,300:2000
    ctenotus,NADH,300:20000
    ctenotus,12S,300:20000
    ctenotus,16S,300:20000
  4. Также копируем туда скрипт, представленный выше и называем его genbank.pl
  5. Далее запускаем командую строку (Start > Run > cmd.exe) и перемещаемся (с помощью команды cd) в рабочую директорию Work

    Сохранение скрипта и файла с условиями в общей папке (рабочей директории)

  6. В рабочей директории Work запускаем наш скрипт:
    D:\work> perl genbank.pl 

    Запуск скрипта на языке Perl

  7. Если все идет корректно, то в окне cmd начнут появлятся наши поисковые запросы:

    Так выглядит запущенный скрипт Perl при получении данных из GENBANK

  8. И в результате в нашей рабочей директории окажутся 4 файла с расширением gb.

    Результирующая директория с загруженными файлами в формате Genbank

Первая часть нашей задачи на сегодня выполнена. Мы загрузили все нужные данные из генбанка. Теперь перед нами стоит следующая задача – конвертация записей из формата genbank в формат fasta (прочитать подробнее о форматах в БД Genbank). Эта операция сопровождается своеобразной чисткой исходного файла – поэтому назвается парсинг. Подробности – во второй части этой практической работы

Скачать готовые файлы в формате .gb можно отсюда

18 thoughts on “Батч-загрузка данных из GENBANK

  1. Гость

    Не могу установить биоперл.
    Downloading BioPerl-Regular Releases packlist…not found

    И с остальными то же самое.
    Может ли это происходить из-за нашего прокси?

    1. Anna

      Вряд ли из-за прокси. А как именно ставите? На ActivePerl?
      Уже смотрели вот тут?

      1. Гость

        На ActivePerl, и через PPM и через командную строку – одно и то же.

        ppm> install PPM-Repositories
        Выдаёт
        Downloading ActiveState Package Repository packlist…not found
        No missing packages to install

        1. Гость

          на nmake15.exe ссылка видимо не рабочая, не удаётся скачать, поэтому этим способом не могу пока попробовать установить.

          1. Anna Post author

            хм, может и правда прокси…
            а если попробовать

            set http_proxy=http://user:passwd@proxy.company.com:8008

  2. Гость

    Через прокси вроде проходит.
    Теперь выдаёт следующее

    ppm> install PPM-Repositories
    No missing packages to install

    ppm> repo add http://bioperl.org/DIST
    ppm repo failed: Repo BioPerl-Regular Releases already set up with URL http://bioperl.org/DIST

    ppm> repo add uwinnipeg
    ppm repo failed: repo.packlist_uri may not be NULL

    ppm> install BioPerl
    ppm install failed: Can’t find any package that provides BioPerl

      1. Гость

        То же самое
        already set up
        и Can’t find any package

        Просто напасть какая-то.
        Если на домашнем компьютере всё установится, значит придётся биоперлом дома пользоваться.

        1. Anna Post author

          мм, а какая версия ActivePerl и заодно что за система?

          1. Гость

            ActivePerl сначала стоял 5.16 – не ставилось
            теперь 5.14.

            Win7 x64

  3. Гость

    Ура, наконец-то получилось поставить биоперл дома (на работе прокси так и не удалось обойти).
    Теперь выдаёт, что ошибка в скрипте
    syntax error at D: work\perl\gebank.pl line 5, near “&gt”
    syntax error at D: work\perl\gebank.pl line 8, near “&lt”
    syntax error at D: work\perl\gebank.pl line 18, near “>”
    syntax error at D: work\perl\gebank.pl line 18, near “$query_string)”
    syntax error at D: work\perl\gebank.pl line 24, near “-format”
    syntax error at D: work\perl\gebank.pl line 24, near “”>$file”)”
    Format not terminated at D: work\perl\gebank.pl line 39, at end of line
    syntax error at D: work\perl\gebank.pl line 39, at EOF
    Execution of D: work\perl\gebank.pl aborted due to compilation errors.

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

    1. Anna Post author

      Скачайте скрипт по ссылке – ftp://www.phylogenetics.ru/practicle7/script.rar (логин – public@phylogenetics.ru; пароль: cornus)
      а то из-за изменения версии вордпресса полетел модуль представления скриптов и туда (в скрипт) записалась куча лишних символов. На страничке я поправлю, но чтобы вам не ждать – загрузите скрипт с фтп.

      1. Гость

        Ура! Всё получилось! Свои группы удалось скачать!
        Шикарно, теперь буду активно использовать.
        Спасибо огромное за ответы и помощь.

  4. Sergey

    Анна, извините за оффтоп. Не подскажете, как сохранить определенный регион из последовательности в GenBank? Мне нужны ITS1ые, а они обычно идут со фланкирующими последовательностями. Как поступают в таком случае?

    1. Anna Post author

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

      Задача исключения фланкирующих кусков (да и любых других не нужных кусков), особенно если речь идет о большом их числе, сводится к тому, чтобы определить “консерватиный” участок и удалить все, что за его пределами. Соответственно, принципиально есть два решения:
      1. сделать multiple sequence alignment к референсной последовательности (которая содержит только интересующий участок) и в ручную “обрезать” все, что за пределами этого участка
      2. делать итеративное сравнение каждой нуклеотидной последовательности с референсной алгоритмом типа BLASTX (если реф. последовательность аминокислотная), записывать координаты кодирующего участка (если нужен он) и потом отрезать все за пределами координат. Для некодирующих участков референсная последовательность остается нуклеотидной и можно брать обычный BLAST

      Первый вариант проще технически, но менее автоматизируемый
      Первый можно сделать с помощью обычных программ типа Clustlw/muscle/mafft и обрезать в любом редакторе (bioedit, seaview,etc). Второй я бы лично делала скриптом на биопитоне-перле и использовала что-то типа exonerate, либо, если установлен blast – то прямо к бласту локально обращалась). Вместо exonerate любой другой софт делающий nucleotide to protein alignment подойдет

      1. Sergey

        Анна, спасибо за консультацию. Буду пробовать.

        1. Anna Post author

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

Comments are closed.