解析GenBank文件:获取Locus Tag VS产品
题
基本上,GenBank文件在基因条目中组成(由“基因”宣布,后跟其相应的'CDS'条目(每个基因只有一个),如下我在此显示。我想在a中获取locus_tag与产品选项卡分隔两个列文件。始终先和后跟空间。
问题是,它似乎是因为“产品”有时在其名称中有时'/'字符,它与此脚本发生冲突,就我所能理解,正在使用'/'作为现场分隔符来存储数组中的信息?
我想解决这个问题,要么修改这个脚本或建立其他脚本。
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
修饰程序添加到Regex以解释MultiLine条目(哪个
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.
不隶属于 StackOverflow