获取文件中最短和最长的序列答案

【问题标题】：Get shortest and longest sequence in file获取文件中最短和最长的序列
【发布时间】：2016-05-29 08:31:38
【问题描述】：

我正在尝试在包含多个类似 genbank 的条目的文件中获取最短和最长的序列。文件示例：

LOCUS       NM_182854               2912 bp    mRNA    linear   PRI 20-APR-2016
DEFINITION  Homo sapiens mRNA.
ACCESSION   NM_182854
SOURCE      Homo sapiens (human)
  ORGANISM  Homo sapiens
            Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;
            Mammalia; Eutheria; Euarchontoglires; Primates; Haplorrhini;
            Catarrhini; Hominidae; Homo.

ORIGIN      
        1 gggcgatcag aagcaggtca cacagcctgt ttcctgtttt caaacgggga acttagaaag
       61 tggcagcccc tcggcttgtc gccggagctg agaaccaaga gctcgaaggg gccatatgac
      //

LOCUS       NM_001323410            6992 bp    mRNA    linear   PRI 20-APR-2016
DEFINITION  Homo sapiens  mRNA.
ACCESSION   NM_001323410
SOURCE      Homo sapiens (human)
  ORGANISM  Homo sapiens
            Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;
            Mammalia; Eutheria; Euarchontoglires; Primates; Haplorrhini;
            Catarrhini; Hominidae; Homo.

ORIGIN      
        1 actacttccg gcttccccgc cccgccccgt ccccgggcgt ctccattttg gtctcaggtg
       61 tggactcggc aagaaccagc gcaagaggga agcagagtta tagctacccc ggc
      //

我想从最短序列和最长序列中打印入藏号、生物的类型

到目前为止我的代码：

#!/usr/bin/perl

use strict;
use warnings;

print "enter file path\n";

while (my $line = <>){
    chomp $line;
    my @record = ($line);

    foreach my $file(@record){
    open(IN, "$file") or die "\n error opening file \n;/\n";

    $/="//";

    while (my $line = <IN>){
        my @gb_seq = split ("ORIGIN", $line);
        my $definition = $gb_seq[0];
        my $sequence = $gb_seq[1];

        $definition =~ m/ORGANISM[\s\t]+(.+)[\n\s\t]+/;
        my $organism = $1;

        if ($definition =~ m/ACCESSION[\s\t]+(\D\D_\d\d\d\d\d\d(\d*))[\n\s\t]+/){
        my $accession = $1;

            $sequence =~ s/\d//g;
            $sequence =~ s/[\n\s\t]//g;
            my $size = length($sequence);
            my @sorted_keys = sort { $a <=> $b } keys my %size;
            my $shortest = $sorted_keys[0];
            my $longest = $sorted_keys[-1];

            print "this is the shortest: $accession $organism size: $shortest\n";
            print "this is the longest: $accession $organism size: $longest\n";
    }
    }}}
    exit;

我想过将长度放入哈希中以获得最短和最长的，但那里有问题。我收到这些错误：

Use of uninitialized value $organism in concatenation (.) or string at test.pl line 39, <IN> chunk 1
Use of uninitialized value $shortest in concatenation (.) or string at test.pl line 39, <IN> chunk 1.
Use of uninitialized value $longest in concatenation (.) or string at test.pl line 40, <IN> chunk 1.

我应该改变哪一部分？谢谢

【问题讨论】：

我在数据中没有看到ORGANISM。也许你的意思是ORIGIN？
您的主要问题是您声明了一个新的空哈希 %size 用于排序命令，这与上面的 $size 标量无关。您需要在 while ($line) 循环上方声明 $biggest_sequence 和 $smallest_sequence 之类的内容，并为每个序列计算它是否应该取代旧的 $biggest_sequence 或 $smallest_sequence。
是的，抱歉，我剪掉了标题，因为它太大了，错过了有机体部分。
好的，我会努力的，谢谢

标签： perl

【解决方案1】：

我们需要找到极端长度的条目，同时能够识别它们所属的记录。通过// 阅读记录又是一个好主意。但是，每条记录都是一个字符串，直接从中提取序列比先将其分成几行更难。因此，我们不妨逐行进行，因为所有需要的东西都有明确的标记。

数据结构的选择很重要，取决于目的。在这里，我将数据组织成一个带有元素的散列，以便于使用它

%block = ( 'accession' => { 'type' => type, 'sequence' => sequence }, ... )

通过“序列”（而不是“加入”）来组织这个数据将极大地帮助在读入数据后执行的搜索，但这会使其非常难以使用。我认为这最终可能会被更多地使用，并且速度的小幅损失并不重要。如果这里的唯一目标是以最佳性能回答特定问题，则其他方法会更合适。注释跟在代码后面。

use warnings;
use strict;
use feature qw(say);

my $file = 'data_seqs.txt';
open my $fh, '<', $file or die "Can't open $file -- $!";

# Hash, helper variables, flag (inside a sequence?), sequence-end marker
my (%block, $accession, $sequence);
my $is_seq = 0;
my $end_marker = qr(\s*//); # marks end of sequence: //

while (my $line = <$fh>) 
{
    chomp($line);
    next if $line =~ /^\s*$/;       # skip empty lines

    if ($line =~ /$end_marker/) {  # done with the sequence
        $is_seq = 0;
        $sequence = ''; 
        next;
    }   

    if ($line =~ /^\s*ACCESSION\s*(\w+)/) { 
        $accession = $1; 
    }   
    elsif ($line =~ /^\s*ORGANISM\s*(.+)/) {
        $block{$accession}{'type'} = $1; 
    }   
    elsif ($line =~ /^\s*ORIGIN/) {  # start sequence on next line
        $is_seq = 1;
    }   
    elsif ($is_seq) {                # read (and add to) sequence
        if ($line =~ /^\s*\d+\s*(.*)/) {
            $block{$accession}{'sequence'} .= $1; 
        }
        else { warn "Not sequence? Line: $line " }
    }   
}

# Identify keys for max and min lenght. Initialize with any keys
my ($max, $min) = keys %block;

foreach my $acc (keys %block) 
{
    my $current_len = length($block{$acc}{'sequence'});
    if ( $current_len > length($block{$max}{'sequence'}) ) { 
        $max = $acc;
    }
    if ( $current_len < length($block{$min}{'sequence'}) ) {
        $min = $acc;
   }
}

say "Maximum length sequence:  ACCESSION: $max, ORGANISM: " . $block{$max}{'type'};
say "Minimum length sequence:  ACCESSION: $min, ORGANISM: " . $block{$min}{'type'};

use Data::Dumper;
print Dumper(\%block);

这个打印（Dumper 的打印输出省略）

最大长度序列：ACCESSION: NM_182854, ORGANISM Homo sapiens 最小长度序列：ACCESSION: NM_001323410, ORGANISM Homo sapiens

关于搜索效率的评论

一种常见的方法是首先构建一个反向查找哈希，然后使用一个库（例如来自 List::Utils 的库）来查找最大值和最小值，然后查找它们所属的位置。为此，我们确实需要构建查找哈希，并且我们会使用该库两次，同时如上所述手动搜索它会遍历结构并且也更简单。另一种选择是将散列顶级键作为序列，然后直接找到最大值和最小值。但是，这样的散列会更难处理。

另一种方法是将数据组织成一个结构，可以更有效地检索此特定信息，可能基于数组。

但是，效率提升似乎并不能证明便利性的巨大损失是合理的。如果速度被证明是一个问题，那么应该考虑这一点。

如果您需要处理多个文件，只需将循环更改为 while (<>) 并在命令行上提交它们。然后将逐行读取所有这些行中的所有行，并且代码保持不变。

可能是我误解了一些术语。我不会从“序列”中删除空格，并且仅在第一行使用单词作为“类型”，只是为了命名几个候选人。这些很容易调整，请告诉我。

【讨论】：

@jnmf 很高兴它很有用。感谢您的反馈:):

【解决方案2】：

您声明您需要两条数据——种质和有机体——用于最长和最短的序列。这意味着您的哈希值需要存储两个元素。此外，当您使用“//”作为记录分隔符时，“//”仍会出现在每条记录的末尾。因此，当您从序列中过滤掉空格和数字时，最后仍然会留下“//”。当我通过调试器运行您的代码时，我发现长度都因此减少了 2。

其他几件事：

使用正则表达式时，请使用“扩展模式”/x，这样您就可以包含空格以提高可读性
您在挖掘 $definition 时假定匹配成功 - 最好测试您的正则表达式并在匹配时分配，在不匹配时死亡
与其将长度存储在哈希中（并丢失序列本身），不如存储序列并稍后计算长度；
我将变量 $line 重命名为 $chunk，因为它包含几行
与计算最短和最长以及打印结果有关的所有事情都需要移出循环。取而代之的是，您只需要在哈希中输入一个条目。如上所述，哈希值需要是一个包含两个值的数组 - 加入和有机体。
您在一个命令中从序列中删除数字，然后在另一个命令中从序列中删除空格 - 不妨将它们一起执行。在我们这样做的时候，不妨删除记录末尾的“/”。

鉴于上面的模组，我明白了；

use v5.14;
use warnings;

print "Enter file path: ";
chomp(my $filename = <>);
open(IN, $filename) or die "\n error opening file \n;/\n";

$/ = "//" ;

my %organisms ;
while (my $chunk = <IN>)  {
    next if $chunk =~ /^\s*\n\s*$/ ;
    my ($definition , $sequence) = split "ORIGIN", $chunk ;

    my $organism ;
    $definition =~ m/ ORGANISM [\s\t]+ (.+) [\n\s\t]+ /x
        ? $organism = $1
        : die "Couldnt find ORGANISM line" ;

    my $accession ;
    $definition =~ m/ ACCESSION [\s\t]+ (\D\D _ \d{6} (\d*))  [\n\s\t]+ /x
        ? $accession = $1
        : die "Cant find ACCESSION line" ;

    $sequence =~ s/[\d\n\s\t\/]//g;
    $organisms{ $sequence } = [ $accession , $organism ] ;
}


my @sorted_keys = sort { length $a  <=>  length $b } keys %organisms ;
my $shortest = $sorted_keys[0];
my $longest  = $sorted_keys[-1];

say "this is the shortest: ",  $organisms{$shortest}->[0],
                        ", ",  $organisms{$shortest}->[1],
                   " size: ",  length $shortest, "\n",
               " sequence: ",  $shortest ;

say  "this is the longest: ",  $organisms{$longest}->[0],
                        ", ",  $organisms{$longest}->[1],
                   " size: ",  length $longest, "\n",
               " sequence: ",  $longest ;

exit;

在您的数据上运行时，它会产生；

$ ./sequence.pl
Enter file path: data.txt
this is the shortest: NM_001323410, Homo sapiens size: 113
 sequence: actacttccggcttccccgccccgccccgtccccgggcgtctccattttggtctcaggtgtggactcggcaagaaccagcgcaagagggaagcagagttatagctaccccggc
this is the longest: NM_182854, Homo sapiens size: 120
 sequence: gggcgatcagaagcaggtcacacagcctgtttcctgttttcaaacggggaacttagaaagtggcagcccctcggcttgtcgccggagctgagaaccaagagctcgaaggggccatatgac

更新上面代码的问题是，如果相同的序列出现在两个块中，那么数据将在哈希中被覆盖并丢失。下面是一个更新版本，它将数据存储在一个数组数组中，这将避免这个问题。它产生完全相同的输出：

use v5.14;
use warnings;

print "Enter file path: ";
chomp(my $filename = <>);
open(IN, $filename) or die "\n error opening file \n;/\n";

$/ = "//" ;

my @organisms ;
while (my $chunk = <IN>)  {
    next if $chunk =~ /^\s*\n\s*$/ ;
    my ($definition , $sequence) = split "ORIGIN", $chunk ;

    my $organism ;
    $definition =~ m/ ORGANISM [\s\t]+ (.+) [\n\s\t]+ /x
        ? $organism = $1
        : die "Couldnt find ORGANISM line" ;

    my $accession ;
    $definition =~ m/ ACCESSION [\s\t]+ (\D\D _ \d{6} (\d*))  [\n\s\t]+ /x
        ? $accession = $1
        : die "Cant find ACCESSION line" ;

    $sequence =~ s/[\d\n\s\t\/]//g;
    push @organisms, [$organism , $accession , $sequence] ;
}


my @sorted_organisms = sort { length $a->[2]  <=>  length $b->[2] }  @organisms ;

my ($organism , $accession , $sequence) = @{ $sorted_organisms[0] };
say "this is the shortest: $accession, $organism, size: ",
    length $sequence, "\n", " sequence: ",  $sequence ;

($organism , $accession , $sequence) = @{ $sorted_organisms[-1] };
say "this is the longest: $accession, $organism, size: ",
    length $sequence, "\n", " sequence: ",  $sequence ;

exit;

【讨论】：