急急急！提取最长序列

赛萌

>EM.1.1 EM.1 LK029911.1 +
CTGTCAATCCCGCACATTAAAGCATCTCTCTCCCCTTGCCTATTTTCTCGCCCCTTCCccatcctccttctctctcgctctctcctaGCATAGCTCCAAGCGCGAGGAACA
>EM.1.2 EM.1 LK029911.1 +
TCCTTTTAAGTACTGCAAACATACCTGGCGGCTTAAGCTTGCAGTATAATGTCAAGTGCTTCCCTACCACTTATTGCTTATAATTGTGATCACCCGTCAAATTGGTCGTTCTTGGAATATAGGCTCTTACTTTTAGAGCATGGTTATTTAAACAGGGGGACGAAAATCCCTCTTTGCTCAATGTAAGTCAGCTAGTTGTACTTACGTTTTGAACTCCACGTCCACACTGGAACAGCATTCGCTATGGACTGTGTCACACGACAAGCAGGTATGTTCATTT
>EM.1.3 EM.1 LK029911.1 +
CTGTCAATCCCGCACATTAAAGCATCTCTCTCCCCTTGCCTATTTTCTCGCCCCTTCCccatcctccttctctctcgctc
很多上面的这些个序列，要求是 EM.1相同的情况下，EM.1.1、EM.1.2、EM.1.3下面对的序列哪个长就把那个提取出来，像这个结果就要是
> EM.1
TCCTTTTAAGTACTGCAAACATACCTGGCGGCTTAAGCTTGCAGTATAATGTCAAGTGCTTCCCTACCACTTATTGCTTATAATTGTGATCACCCGTCAAATTGGTCGTTCTTGGAATATAGGCTCTTACTTTTAGAGCATGGTTATTTAAACAGGGGGACGAAAATCCCTCTTTGCTCAATGTAAGTCAGCTAGTTGTACTTACGTTTTGAACTCCACGTCCACACTGGAACAGCATTCGCTATGGACTGTGTCACACGACAAGCAGGTATGTTCATTT
下面是我写的程序，写的有问题，我自己看不出，因为得出的结果有25个是重复的，就是提出来有25个是提取了两次
#!perl
use strict;
use warnings;
use Data::Dumper;

die "usage $0 input output" unless @ARGV==2;
open IN,$ARGV[0];
my $zn;
my $n;
my %seq;
while(<IN>){
        if(s/^>//){               
                my @its=split;
                $n=$its[0];
                $zn=$its[1];
        }else{
                chomp;
                $seq{$zn}{$n}.=$_;               
        }       
}
close IN;
#print Dumper(%seq),"\n";
open OUT,">","$ARGV[1]";
my %seq1;
my %seq2;
my $seq;
foreach my  $zw (keys %seq){
#        my $seq=$seq{$zw};
        my @length;
        foreach my $w (keys %{$seq{$zw}}){
                $seq=$seq{$zw}{$w};
#                $seq="\U$seq";
                my $length=length($seq);
                push @length,$length;
                $seq1{$zw}{$w}{$length}{$seq}=1;
        }
#        print Dumper(%seq1);
        my $longest=$length[-1];
        foreach(@length){
                $longest=$_ if $_ > $longest;
        }
       
    foreach my $w (keys %{$seq1{$zw}}){
            foreach (keys %{$seq1{$zw}{$w}}){
                    if(exists $seq1{$zw}{$w}{$longest}){
                            my $seqs=Dumper(keys %{$seq1{$zw}{$w}{$longest}});
                            my @seqs=keys %{$seq1{$zw}{$w}{$longest}};
#                            print "@seqs\n";
                            print OUT ">$zw\n@seqs\n";
                    }                                      
            }           
    }      
}   

sunzhiguolu

回复 1# 赛萌
您的代码好多地方, 我还没有接触到. 所以不能为您指出问题, 抱歉.  
在没有更好的解决方案前, 可以试下我这个, 等有了更好的方法再更换好了! (笨方法)

1. #!/usr/bin/perl
   
2. use 5.010;
   
3. use strict;
   
4. use warnings;
   
5. 
   
6. #函数 set_or_get_last_id 用来存放或获取上一次的 $id 的值
   
7. sub set_or_get_last_id{
   
8.         state $last_id;
   
9.         $last_id = shift if (@_);
   
10.         $last_id;
    
11. }
    
12. 
    
13. #哈希变量 %keys_name, %keys_value 并行存储 ($id, $name) 及 ($id, $max_long_str)
    
14. my (%keys_name, %keys_value, $max_long_str);
    
15. while(<>){
    
16.         my ($name, $id);
    
17.         chomp;
    
18.         my @a_line = split /\s+/;
    
19.         if (m/^>/){
    
20.                 ($name, $id) = ($a_line[0], $a_line[1]);
    
21.                 if (!(exists $keys_name{$id})){
    
22.                         $name =~ s/\.[0-9]+$//;    #加入这句
    
23.                         $keys_name{$id} = $name;
    
24.                         $keys_value{$id} = "";
    
25.                         set_or_get_last_id($id);
    
26.                 }
    
27.         }else{
    
28.                 my $last_id = set_or_get_last_id();
    
29.                 $max_long_str = $keys_value{$last_id};
    
30.                 my $current_str = $a_line[0];
    
31.                 if (length($current_str) > length($max_long_str)){
    
32.                         $keys_value{$last_id} = $current_str;
    
33.                 }
    
34.         }
    
35. }
    
36. 
    
37. foreach (keys %keys_name){
    
38.         print "$keys_name{$_}\n$keys_value{$_}\n";
    
39. }
    

复制代码

sunzhiguolu

    另外, 借此机会我也向大家请教一个问题. 我希望将 while 循环封装成一个函数, 怎么将两个哈希变量返回给调用方. 像下面的这种方式吗?

1. sub abc{
   
2.    ...
   
3.    (\%keys_name, \%keys_value);
   
4. }

复制代码

这种方法返回的列表, 随着方法的结束数据是否会受到影响? 还是有更好的方法?
还请大家给予帮助及指点... 多谢大家...

赛萌

回复 2# sunzhiguolu


   艾，你写的我也看不懂

sunzhiguolu

回复 4# 赛萌
比方说? 具体一点... 另外, 代码试过了吗? (我们是希望探讨代码的内容吗? 还是先测试下代码能否得到您的期望的结果.)
   

jason680

回复 1# 赛萌

$ perl max_length.pl
> EM.1
TCCTTTTAAGTACTGCAAACATACCTGGCGGCTTAAGCTTGCAGTATAATGTCAAGTGCTTCCCTACCACTTATTGCTTATAATTGTGATCACCCGTCAAATTGGTCGTTCTTGGAATATAGGCTCTTACTTTTAGAGCATGGTTATTTAAACAGGGGGACGAAAATCCCTCTTTGCTCAATGTAAGTCAGCTAGTTGTACTTACGTTTTGAACTCCACGTCCACACTGGAACAGCATTCGCTATGGACTGTGTCACACGACAAGCAGGTATGTTCATTT

$ cat max_length.pl
use strict;
use warnings;

my $sEM  = "";
my $sSeq = "";
my %hSeq;

sub check_length{
  my ($sEM, $sSeq) = @_;

  return if($sEM eq "");

  $hSeq{$sEM} = $sSeq if(! exists $hSeq{$sEM});
  if(length($hSeq{$sEM}) < length($sSeq)){
    $hSeq{$sEM} = $sSeq;
  }
}


while(<DATA>){
  chomp;
  if(m/^[^>]/){
    $sSeq .= $_;
    next;
  }
  check_length($sEM, $sSeq);

  $sEM = (split)[1];
  $sSeq = "";
}

check_length($sEM, $sSeq);

foreach(sort keys %hSeq){
  print "> $_\n$hSeq{$_}\n";
}

__DATA__
>EM.1.1 EM.1 LK029911.1 +
CTGTCAATCCCGCACATTAAAGCATCTCTCTCCCCTTGCCTATTTTCTCGCCCCTTCCccatcctccttctctctcgctctctcctaGCATAGCTCCAAGCGCGAGGAACA
>EM.1.2 EM.1 LK029911.1 +
TCCTTTTAAGTACTGCAAACATACCTGGCGGCTTAAGCTTGCAGTATAATGTCAAGTGCTTCCCTACCACTTATTGCTTATAATTGTGATCACCCGTCAAATTGGTCGTTCTTGGAATATAGGCTCTTACTTTTAGAGCATGGTTATTTAAACAGGGGGACGAAAATCCCTCTTTGCTCAATGTAAGTCAGCTAGTTGTACTTACGTTTTGAACTCCACGTCCACACTGGAACAGCATTCGCTATGGACTGTGTCACACGACAAGCAGGTATGTTCATTT
>EM.1.3 EM.1 LK029911.1 +
CTGTCAATCCCGCACATTAAAGCATCTCTCTCCCCTTGCCTATTTTCTCGCCCCTTCCccatcctccttctctctcgctc

   

赛萌

回复 5# sunzhiguolu


    恩，试过了，是正确的，万分感谢，辛苦您了。我今天状态不对，没有认真学习的心思，老板今天又死揪着我不放，非常抱歉。

赛萌

回复 6# jason680


    非常不好意思，您指的问题我会一一改正，谢谢。

赛萌

回复 2# sunzhiguolu


    有一个小问题，就是最后的得到的每个序列的表头，后面都多了.1

sunzhiguolu

回复 9# 赛萌
再试下...