急急急！提取最長序列

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]);
   應(yīng)該沒錯吧？

MMMIX

我還是覺得這個應(yīng)該用 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;

如何去學習模塊的知識? 我再具體的問一下怎樣才能知道我在執(zhí)行哪些操作的時候去學習哪些相應(yīng)的模塊? 看到您每次使用一些新的技術(shù)的時候總是感覺很新奇, 但是只能停留在欣賞的狀態(tài), 想學習根本無從下手的感覺, 可能現(xiàn)在的我即使您向我說明也可能不太明白, 您能否指點一下...
   

   

MMMIX

回復 15# sunzhiguolu


    這還真沒有什么速成的方法，只能慢慢積累。你對哪方面感興趣，就多看哪方面的（好）代碼、資料，然后自然的了解的相關(guān)知識、模塊就越來越多。


當然，對于具體的模塊，就用 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 +）

復制代碼