تحليل ملف GenBank:الحصول على علامة الموقع مقابل المنتج
سؤال
في الأساس، يتكون ملف GenBank من إدخالات جينية (يتم الإعلان عنها بواسطة "gene" متبوعة بإدخال "CDS" المقابل لها (واحد فقط لكل جين) مثل المدخلين اللذين أعرضهما هنا أدناه.أرغب في الحصول على locus_tag vs Product في ملف مكون من عمودين مفصول بعلامات جدولة.دائمًا ما يسبق "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"
المحلول
نظرًا لأن نموذج ملف 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 لحساب الإدخالات متعددة الأسطر (والتي leeduhem تمت الإشارة إليه في التعليقات):
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.