急急急！提取最長(zhǎng)序列

賽萌

>EM.1.1 EM.1 LK029911.1 +
CTGTCAATCCCGCACATTAAAGCATCTCTCTCCCCTTGCCTATTTTCTCGCCCCTTCCccatcctccttctctctcgctctctcctaGCATAGCTCCAAGCGCGAGGAACA
>EM.1.2 EM.1 LK029911.1 +
TCCTTTTAAGTACTGCAAACATACCTGGCGGCTTAAGCTTGCAGTATAATGTCAAGTGCTTCCCTACCACTTATTGCTTATAATTGTGATCACCCGTCAAATTGGTCGTTCTTGGAATATAGGCTCTTACTTTTAGAGCATGGTTATTTAAACAGGGGGACGAAAATCCCTCTTTGCTCAATGTAAGTCAGCTAGTTGTACTTACGTTTTGAACTCCACGTCCACACTGGAACAGCATTCGCTATGGACTGTGTCACACGACAAGCAGGTATGTTCATTT
>EM.1.3 EM.1 LK029911.1 +
CTGTCAATCCCGCACATTAAAGCATCTCTCTCCCCTTGCCTATTTTCTCGCCCCTTCCccatcctccttctctctcgctc
很多上面的這些個(gè)序列，要求是 EM.1相同的情況下，EM.1.1、EM.1.2、EM.1.3下面對(duì)的序列哪個(gè)長(zhǎng)就把那個(gè)提取出來(lái)，像這個(gè)結(jié)果就要是
> EM.1
TCCTTTTAAGTACTGCAAACATACCTGGCGGCTTAAGCTTGCAGTATAATGTCAAGTGCTTCCCTACCACTTATTGCTTATAATTGTGATCACCCGTCAAATTGGTCGTTCTTGGAATATAGGCTCTTACTTTTAGAGCATGGTTATTTAAACAGGGGGACGAAAATCCCTCTTTGCTCAATGTAAGTCAGCTAGTTGTACTTACGTTTTGAACTCCACGTCCACACTGGAACAGCATTCGCTATGGACTGTGTCACACGACAAGCAGGTATGTTCATTT
下面是我寫的程序，寫的有問(wèn)題，我自己看不出，因?yàn)榈贸龅慕Y(jié)果有25個(gè)是重復(fù)的，就是提出來(lái)有25個(gè)是提取了兩次
#!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

回復(fù) 1# 賽萌
您的代碼好多地方, 我還沒(méi)有接觸到. 所以不能為您指出問(wèn)題, 抱歉.  
在沒(méi)有更好的解決方案前, 可以試下我這個(gè), 等有了更好的方法再更換好了! (笨方法)

1. #!/usr/bin/perl
   
2. use 5.010;
   
3. use strict;
   
4. use warnings;
   
5. 
   
6. #函數(shù) set_or_get_last_id 用來(lái)存放或獲取上一次的 $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 并行存儲(chǔ) ($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. }
    

復(fù)制代碼

sunzhiguolu

    另外, 借此機(jī)會(huì)我也向大家請(qǐng)教一個(gè)問(wèn)題. 我希望將 while 循環(huán)封裝成一個(gè)函數(shù), 怎么將兩個(gè)哈希變量返回給調(diào)用方. 像下面的這種方式嗎?

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

復(fù)制代碼

這種方法返回的列表, 隨著方法的結(jié)束數(shù)據(jù)是否會(huì)受到影響? 還是有更好的方法?
還請(qǐng)大家給予幫助及指點(diǎn)... 多謝大家...

賽萌

回復(fù) 2# sunzhiguolu


   艾，你寫的我也看不懂

sunzhiguolu

回復(fù) 4# 賽萌
比方說(shuō)? 具體一點(diǎn)... 另外, 代碼試過(guò)了嗎? (我們是希望探討代碼的內(nèi)容嗎? 還是先測(cè)試下代碼能否得到您的期望的結(jié)果.)
   

jason680

回復(fù) 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

   

賽萌

回復(fù) 5# sunzhiguolu


    恩，試過(guò)了，是正確的，萬(wàn)分感謝，辛苦您了。我今天狀態(tài)不對(duì)，沒(méi)有認(rèn)真學(xué)習(xí)的心思，老板今天又死揪著我不放，非常抱歉。

賽萌

回復(fù) 6# jason680


    非常不好意思，您指的問(wèn)題我會(huì)一一改正，謝謝。

賽萌

回復(fù) 2# sunzhiguolu


    有一個(gè)小問(wèn)題，就是最后的得到的每個(gè)序列的表頭，后面都多了.1

sunzhiguolu

回復(fù) 9# 賽萌
再試下...