Разбор файла GenBank: получение тега локуса по сравнению с продуктом

По сути, файл GenBank состоит из записей генов (объявленных «геном», за которым следует соответствующая запись «CDS» (только одна для каждого гена), как две, которые я показываю здесь ниже. Я хотел бы получить locus_tag против продукта в разделенном табуляцией файл с двумя столбцами. "gene" и "CDS" всегда предшествуют и сопровождаются пробелами.

В предыдущем вопросе был предложен сценарий.

Проблема в том, что кажется, что, поскольку «продукт» иногда имеет символ «/» внутри своего имени, он конфликтует с этим сценарием, который, насколько я понимаю, использует «/» в качестве разделителя полей для хранения информации в множество?

Я хотел бы решить эту проблему, либо изменив этот скрипт, либо создав другой.

perl -nE'
  BEGIN{ ($/, $") = ("CDS", "\t") }
  say "@r[0,1]" if @r= m!/(?:locus_tag|product)="(.+?)"!g and @r>1
' file


 gene            complement(8972..9094)
                 /locus_tag="HAPS_0004"
                 /db_xref="GeneID:7278619"
 CDS             complement(8972..9094)
                 /locus_tag="HAPS_0004"
                 /codon_start=1
                 /transl_table=11
                 /product="hypothetical protein"
                 /protein_id="YP_002474657.1"
                 /db_xref="GI:219870282"
                 /db_xref="GeneID:7278619"
                 /translation="MYYKALAHFLPTLSTMQNILSKSPLSLDFRLLFLAFIDKR"
 gene            68..637
                 /locus_tag="HPNK_00040"
 CDS             68..637
                 /locus_tag="HPNK_00040"
                 /codon_start=1
                 /transl_table=11
                 /product="NinG recombination protein/bacteriophage lambda
                 NinG family protein"
                 /protein_id="CRESA:HPNK_00040"
                 /translation="MIKPKVKKRKCKCCGGEFKSADSFRKWCSAECGVKLAKIAQEKA
                 RQKAIEKRNREERAKIKATRERLKSRSEWLKDAQAIFNEYIRLRDKDEPCISCRRFHQ
                 GQYHAGHYRTVKAMPELRFNEDNVHKQCSACNNHLSGNITEYRINLVRKIGAERVEAL
                 ESYHPPVKWSVEDCKEIIKTYRAKIKELK"

python parsing perl bioperl genbank
person biotech    schedule 27.02.2014    source источник
comment
Вы понимаете неправильно. Единственная проблема с / заключается в том, что он конфликтует с разделителем по умолчанию оператора сопоставления m//, но это уже устранено путем замены разделителя на !, как в m!!. Это не имеет ничего общего с массивами или разделителями полей.   -  person TLP    schedule 27.02.2014
comment
Этот формат (GenBank) кажется своего рода стандартным форматом, поэтому я уверен, что есть модуль, который будет анализировать его для вас, что, скорее всего, будет проще и безопаснее, чем то быстрое исправление, которое вы пытаетесь здесь.   -  person TLP    schedule 27.02.2014
comment
Вот один из них: Bio::GenBankParser, например.   -  person TLP    schedule 27.02.2014
comment
Bio::GenBankParser просто разбирает записи GenBank в формат YAML, та же проблема, по крайней мере, для меня.   -  person biotech    schedule 27.02.2014
comment
Нет, это не правда. Я не знаю, откуда вы это взяли, в документации точно не упоминается, что он анализирует GenBank в YAML. Однако в другом результате поиска, genbank-parser.pl, указано parse GenBank records into YAML. Возможно, вы перепутали два результата поиска? Насколько я могу судить, это будет анализировать запись GenBank в хэш-ссылку, например. print $rec->{'ACCESSION'};   -  person TLP    schedule 27.02.2014
comment
Это нормально, что второй product имеет новую строку в своем значении?   -  person Lee Duhem    schedule 27.02.2014
comment
Да, там новая линия, лидухем   -  person biotech    schedule 27.02.2014
comment
metacpan.org/pod/Bio::Perl кажется полезным для некоторых людей, смехотворное количество заглавных букв в своих вопросах. И это, по-видимому, позволяет некоторые махинации с чем-то, что называется форматом генбанка.   -  person DeVadder    schedule 27.02.2014
comment
@popnard Я добавил простой скрипт Perl, который анализирует файл GenBank и распечатывает записи, которые вам нужны. Похоже, это тот же результат, что и ваша сломанная однострочная строка, если она будет исправлена, но извлечена более безопасным и простым способом. И это не YAML.   -  person TLP    schedule 27.02.2014
comment
Если вы хотите анализировать файлы GenBank (которые являются стандартным форматом данных последовательности) с помощью Perl, вам действительно следует использовать BioPerl, а именно Bio::SeqIO. См. этот HOWTO — bioperl.org/wiki/HOWTO:SeqIO и соответствующую аннотацию функций. HOWTO — bioperl.org/wiki/HOWTO:Feature-Annotation.   -  person neilfws    schedule 28.02.2014

Ответы (2)


arrow_upward
2
arrow_downward

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

Используя этот код и модуль Bio::GenBankParser, он был проанализирован, чтобы определить, какие части структуры вы были после. В данном случае «функции», содержащие как поле locus_tag, так и поле product.

use strict;
use warnings;
use feature 'say';
use Bio::GenBankParser;

my $file = shift;
my $parser = Bio::GenBankParser->new( file => $file );
while ( my $seq = $parser->next_seq ) {
    my $feat = $seq->{'FEATURES'};
    for my $f (@$feat) {
        my $tag = $f->{'feature'}{'locus_tag'};
        my $prod = $f->{'feature'}{'product'};
        if (defined $tag and defined $prod) {
            say join "\t", $tag, $prod;
        }
    }
}

Использование:

perl script.pl input.txt > output.txt

Вывод:

MG_001  DNA polymerase III, beta subunit
MG_470  CobQ/CobB/MinD/ParA nucleotide binding domain-containing protein

Результат вашего однострочного ввода для того же ввода будет:

MG_001  DNA polymerase III, beta subunit
MG_470  CobQ/CobB/MinD/ParA nucleotide binding
                     domain-containing protein

Предполагая, конечно, что вы добавляете модификатор /s в регулярное выражение для учета многострочных записей (что leeduhem указано в комментариях) :

m!/(?:locus_tag|product)="(.+?)"!sg
#                                ^---- this
person TLP    schedule 27.02.2014
comment
Модуль установлен, но у меня есть это сообщение об ошибке: «синтаксическая ошибка в /Users/bernardo/Documents/BioLinux/A0_scripts/parse_gbk_2.pl, строка 13, рядом с надписью join Execution of /Users/bernardo/Documents/BioLinux/A0_scripts/parse_gbk_2. pl прерван из-за ошибок компиляции.' - person biotech; 27.02.2014
comment
Да, вам нужно use feature 'say', чтобы использовать say(). - person TLP; 27.02.2014
comment
@popnard Всегда лучше использовать модуль. Большая часть этого кода была взята из документации к модулю, кроме распечатки. - person TLP; 27.02.2014
comment
ОК, сработало для одного файла, но не для других: ' ОШИБКА (строка 5): недопустимый раздел: ожидалась закомментированная строка, или заголовок, или локус, или источник базы данных, или определение, или строка доступа, или строка проекта, или строка версии. , или ключевые слова, или исходная строка, или организм, или ссылка, или признаки, или число оснований, или контиг, или происхождение, или комментарий, или разделитель записи. вместо 59273' - person biotech; 27.02.2014
comment
Неработающий файл находится здесь: drive.google.com/file/ d/0B8-ZAuZe8jldcTNZUzdXd0lCREU/ - person biotech; 27.02.2014
comment
Этот сработал, кажется, это строка 5, где скрипт ожидает чего-то, чего там нет. drive.google.com/file/d/0B8-ZAuZe8jldVV9ObHV6NnBMOWM/ - person biotech; 27.02.2014
comment
Удалена проблемная строка и, кажется, работает хорошо: Проект DBLINK: 59273 Биопроект: PRJNA59273 - person biotech; 27.02.2014
comment
Я ничего не знаю о формате, поэтому ничего не могу сказать об этом. Но похоже, что ваши файлы могут быть каким-то образом повреждены. - person TLP; 27.02.2014

arrow_upward
1
arrow_downward

Прочитав ваш продублированный вопрос http://www.biostars.org/p/94164/ ( пожалуйста, не дублируйте подобные посты), вот минимальный ответ Biopython:

import sys
from Bio import SeqIO
filename = sys.argv[1] # Takes first command line argument input filename
for record in SeqIO.parse(filename, "genbank"):
    for feature in record.features:
        if feature.type == "CDS":
            locus_tag = feature.qualifiers.get("locus_tag", ["???"])[0]
            product = feature.qualifiers.get("product", ["???"])[0]
            print("%s\t%s" % (locus_tag, product))

С небольшими изменениями вы можете вместо этого записать это в файл.

person peterjc    schedule 28.02.2014