Análise GenBank arquivo:obter locus tag vs produto
Pergunta
Basicamente, um GenBank arquivo consiste no gene entradas (anunciado pela 'gene', seguido por seu correspondente 'CDS' entrada (somente uma por gene), como os dois que eu mostrar aqui abaixo.Eu gostaria de obter locus_tag vs produto em uma guia, delimitado por duas colunas do arquivo.'gene' e 'CDS' são sempre precedido e seguido por espaços.
Uma pergunta anterior sugeriu um script.
O problema é que parece que porque o "produto" tem, por vezes, o carácter '/' dentro de seu nome, seus conflitos com este script, que, até onde posso entender, é usando '/' como separador de campo para armazenar informações em uma matriz?
Eu gostaria de resolver este problema, modificar esse script ou a construção de outro.
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"
Solução
Como o exemplo do GenBank arquivo incompleto, fui on-line para encontrar um arquivo de exemplo que pode ser usado em um exemplo, e eu encontrei este arquivo.
Usando este código e o Bio::GenBankParser
módulo, foi analisada a adivinhar o que as partes da estrutura que foram depois.Neste caso, "recursos" que continha uma locus_tag
campo e um product
de campo.
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;
}
}
}
Uso:
perl script.pl input.txt > output.txt
Saída:
MG_001 DNA polymerase III, beta subunit
MG_470 CobQ/CobB/MinD/ParA nucleotide binding domain-containing protein
A saída do seu one-liner para a mesma entrada seria:
MG_001 DNA polymerase III, beta subunit
MG_470 CobQ/CobB/MinD/ParA nucleotide binding
domain-containing protein
Supondo, claro, que você adicionar o /s
modificador para a regex para dar conta de várias linhas de entradas (o que leeduhem apontado nos comentários):
m!/(?:locus_tag|product)="(.+?)"!sg
# ^---- this
Outras dicas
Depois de ter lido o seu duplicado pergunta http://www.biostars.org/p/94164/ (por favor, não fazer double post como este), aqui está um mínimo de Biopython resposta:
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))
Com pequenas mudanças você pode escrever isso para um ficheiro em vez disso.