Разбор файла GenBank:получить тег локуса по сравнению с продуктом
Вопрос
По сути, файл GenBank состоит из записей генов (объявленных "gene", за которыми следует соответствующая запись "CDS" (только по одной на ген), подобных двум, которые я показываю здесь ниже.Я хотел бы получить locus_tag vs product в файле с двумя столбцами, разделенными табуляцией."gene" и "CDS" всегда предваряются пробелами и следуют за ними.
В предыдущем вопросе предлагался сценарий.
Проблема в том, что кажется, что из-за того, что 'product' иногда имеет символ '/' внутри своего имени, у него возникают конфликты с этим скриптом, который, насколько я могу понять, использует '/' в качестве разделителя полей для хранения информации в массиве?
Я хотел бы решить эту проблему, либо изменив этот скрипт, либо создав другой.
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"
Решение
Поскольку ваш образец файла 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
модификатор регулярного выражения для учета многострочных записей (которые лидухем указано в комментариях):
m!/(?:locus_tag|product)="(.+?)"!sg
# ^---- this
Другие советы
Having read your duplicated question http://www.biostars.org/p/94164/ (please don't double post like this), here's a minimal Biopython answer:
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))
With minor changes you can write this out to a file instead.