亚洲av成人无遮挡网站在线观看,少妇性bbb搡bbb爽爽爽,亚洲av日韩精品久久久久久,兔费看少妇性l交大片免费,无码少妇一区二区三区

<label id="7ucdm"></label>

Chinaunix

標(biāo)題: 急急急！提取最長(zhǎng)序列 [打印本頁(yè)]

---

作者: 賽萌    時(shí)間: 2015-09-28 13:42
標(biāo)題: 急急急！提取最長(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    時(shí)間: 2015-09-28 15:27


回復(fù) 1# 賽萌
您的代碼好多地方, 我還沒有接觸到. 所以不能為您指出問(wèn)題, 抱歉.  
在沒有更好的解決方案前, 可以試下我這個(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    時(shí)間: 2015-09-28 15:38


    另外, 借此機(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)... 多謝大家...

---

作者: 賽萌    時(shí)間: 2015-09-28 15:52
回復(fù) 2# sunzhiguolu


   艾，你寫的我也看不懂

---

作者: sunzhiguolu    時(shí)間: 2015-09-28 15:56


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

---

作者: jason680    時(shí)間: 2015-09-28 16:01


回復(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

   

---

作者: 賽萌    時(shí)間: 2015-09-28 16:21
回復(fù) 5# sunzhiguolu


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

---

作者: 賽萌    時(shí)間: 2015-09-28 16:22
回復(fù) 6# jason680


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

---

作者: 賽萌    時(shí)間: 2015-09-28 17:24
回復(fù) 2# sunzhiguolu


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

---

作者: sunzhiguolu    時(shí)間: 2015-09-28 17:32
回復(fù) 9# 賽萌
再試下...

   

---

作者: sunzhiguolu    時(shí)間: 2015-09-28 18:08


回復(fù) 9# 賽萌
    抱歉, 在看問(wèn)題的時(shí)候理解的有問(wèn)題. 我將代碼修改了一下, 再試試...

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. }

復(fù)制代碼

---

作者: 賽萌    時(shí)間: 2015-09-28 18:14
回復(fù) 11# sunzhiguolu

我直接把這一行($name, $id) = ($a_line[0], $a_line[1]);改成了($name, $id) = ($a_line[1], $a_line[1]);
   應(yīng)該沒錯(cuò)吧？

---

作者: MMMIX    時(shí)間: 2015-09-28 21:09
我還是覺得這個(gè)應(yīng)該用 Bio::Seq 來(lái)讀取。

---

作者: MMMIX    時(shí)間: 2015-09-28 21:48


回復(fù) 1# 賽萌


    試試這個(gè)用 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    時(shí)間: 2015-09-28 22:49


回復(fù) 14# MMMIX
    向您請(qǐng)教一個(gè)問(wèn)題,

use Bio::SeqIO;

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

   

---

作者: MMMIX    時(shí)間: 2015-09-28 22:59


回復(fù) 15# sunzhiguolu


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


當(dāng)然，對(duì)于具體的模塊，就用 perldoc 查看其文檔，例如 perldoc Bio::SeqIO。

---

作者: sunzhiguolu    時(shí)間: 2015-09-28 23:03
回復(fù) 16# MMMIX
    我懂了, 感謝您的幫助及支持... 多謝多謝...

   

---

作者: b114213903    時(shí)間: 2015-09-29 09:55


回復(fù) 14# MMMIX


    if判斷很好，哈哈，學(xué)習(xí)了！


有一個(gè)問(wèn)題是這樣輸出的序列名還是原名，與樓主要求稍有差異！

可以加上兩行代碼：

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

復(fù)制代碼

---

歡迎光臨 Chinaunix (http://www.72891.cn/)

Powered by Discuz! X3.2

感谢您访问我们的网站，您可能还对以下资源感兴趣：

亚洲av成人无遮挡网站在线观看

<dfn id="g9djc"></dfn>