急急急！提取最长序列

sunzhiguolu

回复 9# 赛萌
    抱歉, 在看问题的时候理解的有问题. 我将代码修改了一下, 再试试...

1. #!/usr/bin/perl
   
2. use 5.010;
   
3. use strict;
   
4. use warnings;
   
5. 
   
6. sub set_or_get_last_id{
   
7.         state $last_id;
   
8.         $last_id = shift if (@_);
   
9.         $last_id;
   
10. }
    
11. 
    
12. my %keys;
    
13. while (<>){
    
14.         chomp;
    
15.         my @a_line = split /\s+/;
    
16.         my (undef, $id) = (@a_line);
    
17.         if (m/^>/ and !(exists $keys{$id})){
    
18.                 set_or_get_last_id($id);
    
19.                 $keys{$id} = undef;
    
20.         }else{
    
21.                 my $current_str = $a_line[0];
    
22.                 my $last_id = set_or_get_last_id();
    
23.                 $keys{$last_id} = $current_str if (length($current_str) > length($keys{$last_id}));
    
24.         }
    
25. }
    
26. 
    
27. while (my ($id, $value) = each %keys){
    
28.         say "> $id\n$value";
    
29. }

复制代码

赛萌

回复 11# sunzhiguolu

我直接把这一行($name, $id) = ($a_line[0], $a_line[1]);改成了($name, $id) = ($a_line[1], $a_line[1]);
   应该没错吧？

MMMIX

我还是觉得这个应该用 Bio::Seq 来读取。

MMMIX

回复 1# 赛萌


    试试这个用 Bio::SeqIO 的版本：

#!/usr/bin/perl

use strict;
use warnings;

use v5.14;
use autodie;
use Data::Dumper;

use Bio::SeqIO;

my $in  = Bio::SeqIO->new( '-format' => 'fasta', '-file' => $ARGV[0] );
my $out = Bio::SeqIO->new( '-format' => 'fasta', '-fh'   => \*STDOUT );

my %longest;

while (my $seq = $in->next_seq()) {
  my $id = $seq->id();
  $id =~ s/\.[^.]+$//;
  if (not exists $longest{$id} or $seq->length() > $longest{$id}->length()) {
    $longest{$id} = $seq;
  }
}

for my $seq (values %longest) {
  $out->write_seq($seq);
}

sunzhiguolu

回复 14# MMMIX
    向您请教一个问题,

use Bio::SeqIO;

如何去学习模块的知识? 我再具体的问一下怎样才能知道我在执行哪些操作的时候去学习哪些相应的模块? 看到您每次使用一些新的技术的时候总是感觉很新奇, 但是只能停留在欣赏的状态, 想学习根本无从下手的感觉, 可能现在的我即使您向我说明也可能不太明白, 您能否指点一下...
   

   

MMMIX

回复 15# sunzhiguolu


    这还真没有什么速成的方法，只能慢慢积累。你对哪方面感兴趣，就多看哪方面的（好）代码、资料，然后自然的了解的相关知识、模块就越来越多。


当然，对于具体的模块，就用 perldoc 查看其文档，例如 perldoc Bio::SeqIO。

sunzhiguolu

回复 16# MMMIX
    我懂了, 感谢您的帮助及支持... 多谢多谢...

   

b114213903

回复 14# MMMIX


    if判断很好，哈哈，学习了！


有一个问题是这样输出的序列名还是原名，与楼主要求稍有差异！

可以加上两行代码：

1. $seq->id($id);                                #修改序列名（如: EM.1.1  => EM.1）
   
2. $seq->desc('');                                #删除序列注释信息（如：EM.1 LK029911.1 +）

复制代码