如何用Perl从GenBank文件中获取COMMENT之前的注释数据

问题说明 #

需要编写Perl脚本从GenBank文件中提取COMMENT行之前的注释信息，现有脚本会错误包含COMMENT行的部分内容，同时需利用GenBank文件首行以LOCUS开头的特性优化代码。

输入示例 #

LOCUS       mitochondrion_genome 19524 bp DNA HTG 17-DEC-2022
DEFINITION  Drosophila melanogaster primary_assembly mitochondrion_genome BDGP6.32 full
            sequence 1..19524 reannotated via EnsEMBL
ACCESSION   primary_assembly:BDGP6.32:mitochondrion_genome:1:19524:1
VERSION     mitochondrion_genomeBDGP6.32
KEYWORDS    .
SOURCE      fruit fly
  ORGANISM  Drosophila melanogaster
            Eukaryota; Opisthokonta; Metazoa; Eumetazoa; Bilateria;
            Protostomia; Ecdysozoa; Panarthropoda; Arthropoda; Mandibulata;
            Pancrustacea; Hexapoda; Insecta; Dicondylia; Pterygota; Neoptera;
            Endopterygota; Diptera; Brachycera; Muscomorpha; Eremoneura;
            Cyclorrhapha; Schizophora; Acalyptratae; Ephydroidea;
            Drosophilidae; Drosophilinae.
COMMENT     This sequence was annotated by FlyBase (https://www.flybase.org). Please visit the
            Ensembl or EnsemblGenomes web site, http://www.ensembl.org/ or
            http://www.ensemblgenomes.org/ for more information.

当前脚本 #

$genbank = <STDIN>;
chomp ($genbank);
open (READ, "<$genbank") or die;
@data = <READ>;
close READ;

$end= $#data;
for ($line= 0; $line<= $end; $line++){
    if ($data[$line] =~ /LOCUS/){
        @annotation = (@annotation, $data[$line]);
        until ($data[$line] =~ /COMMENT/){
            $line++;
            @annotation = (@annotation, $data[$line]);
}}
print @annotation;

错误输出 #

LOCUS       mitochondrion_genome 19524 bp DNA HTG 17-DEC-2022
DEFINITION  Drosophila melanogaster primary_assembly mitochondrion_genome BDGP6.32 full
            sequence 1..19524 reannotated via EnsEMBL
ACCESSION   primary_assembly:BDGP6.32:mitochondrion_genome:1:19524:1
VERSION     mitochondrion_genomeBDGP6.32
KEYWORDS    .
SOURCE      fruit fly
  ORGANISM  Drosophila melanogaster
            Eukaryota; Opisthokonta; Metazoa; Eumetazoa; Bilateria;
            Protostomia; Ecdysozoa; Panarthropoda; Arthropoda; Mandibulata;
            Pancrustacea; Hexapoda; Insecta; Dicondylia; Pterygota; Neoptera;
            Endopterygota; Diptera; Brachycera; Muscomorpha; Eremoneura;
            Cyclorrhapha; Schizophora; Acalyptratae; Ephydroidea;
            Drosophilidae; Drosophilinae.
COMMENT     This sequence was annotated by FlyBase (https://www.flybase.org). Please visit the

修正方案及优化后的脚本 #

问题根源 #

原脚本的until循环逻辑错误：先执行$line++再将行添加到数组，导致当遇到COMMENT行时，已经将该行加入了结果中。同时一次性读取整个文件到数组的方式对大文件不友好。

优化后的脚本 #

#!/usr/bin/perl
use strict;
use warnings;

# 读取输入的文件名
my $genbank = <STDIN>;
chomp $genbank;

# 打开文件（使用词法文件句柄，更安全）
open my $fh, '<', $genbank or die "无法打开文件: $!";

my @annotation;
my $start_capture = 0;

while (my $line = <$fh>) {
    # 匹配首行LOCUS，开始捕获注释
    if ($line =~ /^LOCUS/) {
        $start_capture = 1;
        push @annotation, $line;
        next;
    }
    
    # 遇到COMMENT行，立即停止捕获
    if ($line =~ /^COMMENT/) {
        last;
    }
    
    # 处于捕获状态时，将行加入结果
    push @annotation, $line if $start_capture;
}

close $fh;

# 输出结果
print @annotation;

优化点说明 #

1. 启用严格模式：添加use strict; use warnings;，避免未声明变量、语法错误等问题，提升代码健壮性。
2. 逐行处理文件：无需一次性加载整个文件到内存，处理大GenBank文件时更高效。
3. 修正捕获逻辑：遇到COMMENT行直接终止循环，完全避免包含COMMENT相关内容。
4. 利用LOCUS特性：精准定位注释起始点，跳过文件开头可能存在的无关内容（如果有的话）。
5. 安全文件操作：使用词法文件句柄$fh替代旧的裸句柄，避免文件句柄冲突。

内容的提问来源于stack exchange，提问作者user21991719