Analyse der GenBank-Datei:holen Sie sich locus Tag vs Produkt
Frage
Grundsätzlich besteht eine GenBank-Datei aus Geneinträgen (angekündigt durch 'gene' gefolgt von dem entsprechenden 'CDS' -Eintrag (nur einer pro Gen), wie die beiden, die ich hier unten zeige.Ich möchte locus_tag vs product in einer tabulatorgetrennten zweispaltigen Datei erhalten.'gene' und 'CDS' werden immer mit Leerzeichen vorangestellt und gefolgt.
Eine vorherige Frage schlug ein Skript vor.
Das Problem ist, dass es den Anschein hat, dass, weil 'Produkt' manchmal '/' Zeichen in seinem Namen hat, Konflikte mit diesem Skript auftreten, das, soweit ich verstehen kann, '/' als Feldtrennzeichen zum Speichern von Informationen verwendet in einem Array?
Ich möchte das lösen, indem ich entweder dieses Skript ändere oder ein anderes erstelle.
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"
Lösung
Da Ihre Beispiel-GenBank-Datei unvollständig war, ging ich online, um eine Beispieldatei zu finden, die in einem Beispiel verwendet werden könnte, und ich fand diese Datei.
Mit diesem Code und dem Bio::GenBankParser
modul, es wurde analysiert, um zu erraten, nach welchen Teilen der Struktur Sie suchten.In diesem Fall "Features", die sowohl eine locus_tag
feld und ein product
Feld.
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;
}
}
}
Nutzung:
perl script.pl input.txt > output.txt
Ausgabe:
MG_001 DNA polymerase III, beta subunit
MG_470 CobQ/CobB/MinD/ParA nucleotide binding domain-containing protein
Die Ausgabe von Ihrem Einzeiler für die gleiche Eingabe wäre:
MG_001 DNA polymerase III, beta subunit
MG_470 CobQ/CobB/MinD/ParA nucleotide binding
domain-containing protein
Vorausgesetzt natürlich, dass Sie das hinzufügen /s
Modifikator für den regulären Ausdruck, um mehrzeilige Einträge zu berücksichtigen (die leeduhem in den Kommentaren darauf hingewiesen):
m!/(?:locus_tag|product)="(.+?)"!sg
# ^---- this
Andere Tipps
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.