解析GenBank文件：獲取基因座標籤與產物

-1

基本上，GenBank文件由基因條目組成（由'基因'及其對應的'CDS'條目（每個基因只有一個）發佈，就像我在下面展示的兩個基因條目。想獲得locus_tag VS產品在製表符分隔的兩列的文件。「基因」和「CDS」的空間總是前面和後面。解析GenBank文件：獲取基因座標籤與產物

A previous question suggested a script.

的問題是，它似乎是因爲「產品'有時'/'字符在它的名字裏面，它與這個腳本有衝突，據我所知，它使用'/'作爲字段分隔符來存儲信息到一個數組中？

我想解決這個問題，無論是修改這個腳本或建立其他的。

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"

來源

2014-02-27 biotech

你明白錯誤。 '/'唯一的問題就是它與匹配運算符'm //'的默認分隔符衝突，但是已經通過用'！'替換分隔符來處理，就像'm !!'中那樣。它與數組或字段分隔符無關。 – TLP

這種格式（GenBank）似乎是某種標準格式，所以我敢打賭有一個模塊可以爲你解析它，這很可能比你在這裏嘗試的這種快速修復更簡單也更安全。 – TLP

例如，下面是一個：['Bio :: GenBankParser']（http://search.cpan.org/perldoc?Bio%3A%3AGenBankParser）。 – TLP

當你的樣本基因庫文件是不完整的，我去網上找，可以在示例中使用的樣本文件，我發現this file。

使用此代碼和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改性劑添加到正則表達式佔多項（其中leeduhem在評論中指出）：

m!/(?:locus_tag|product)="(.+?)"!sg 
#        ^---- this

來源

2014-02-27 13:06:15 TLP

Module is installed，but I error error：'syntax error at /Users/bernardo/Documents/BioLinux/A0_scripts/parse_gbk_2.pl line 13，near「say join」 /Users/bernardo/Documents/BioLinux /由於編譯錯誤，A0_scripts/parse_gbk_2.pl中止。' – biotech

哦，對，你需要使用'say（'）來使用'say（）'。 – TLP

好人給出很好的答案，不斷學習！ – biotech

看了你的重複問題http://www.biostars.org/p/94164/（請不要雙擊後像這一點），這裏有一個最小的Biopython答案：

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))

有了細微的變化，你可以這樣寫出來的文件，而不是。

來源

2014-02-28 16:16:31 peterjc

解析GenBank文件：獲取基因座標籤與產物

回答

相關問題