perl 比對，出不來結(jié)果，不知如何修改

1002279625

<IN1>:>AT1G50920.1 | [1-104] [403-416] | chr1:18870139-18870554 FORWARD LENGTH=118
GGTTGGGTCAGAAAACAGTCGACCGTCATTTGGTTGGGTCGGCTCGCTTCTCTTATCAGCTAACCCTAAGTACGCCGACTTTCAAATTCGC
CACTCCCGTTCATCTTTTCTCTGCATCGATAAGCGAAG
>AT1G15970.1 | [1-374] | chr1:5488495-5488868 REVERSE LENGTH=374
TCGTTTCGTCGTCGATCAGAAAAAGCTTCTTTCTTTCTTTTTTTTTTCTGGGTTTTTCTGAAGTCAAAAGCTTTTTTTTT
CTTCATTGTTCAAAGATTTTATTTTGTTGCTCGCTTGGGTCAACCTTGCTTTTCTCGTGAACACAGCAAGCAACTTGGTT
TAGTTTTACTAGAAGATTCCTAAAAATGAGATATTTAGCGGTGACAAAAAAGTTTCATTTTTGGTGCTAATGTAGCTGTA
AAGTCTAGGCCTTCCGATATATAATCTCTTGAAATAAGTTTCGAGACTGTTGGAAATAGCAATTTTGTTGTTGTTGTGTG
ATCCAGTGATAGATCTGTGAAATCAAAAAGCAGTTTTTTGTGAGAATCGGGTCA
>AT1G73440.1 | [1-55] | chr1:27611363-27611417 FORWARD LENGTH=55
ATTGAAAAGAAAACACATCCCAACTCTGGAGAAACAAATCCCTAATTAGCTAACA
>AT1G75120.1 | [1-146] | chr1:28198657-28198802 REVERSE LENGTH=146
CCGTGTTAAAAAACACGTTGGGTCAAAATCAATGGTGGGTGGGTCTTAACTTGCACTACTAGGATTTGATCACAGAAGAAA
GAAAGGGGTTTATTAATTTGGTTACTGTATCTTACTCTTTTTACTTTCGGATAACTTCAATTGGCA
<IN>:GGTTGGGTC
[ATC]GTTGGGTC
G[ATC]TTGGGTC
GG[ACG]TGGGTC
GGT[ACG]GGGTC
GGTT[ATC]GGTC
GGTTG[ATC]GTC
GGTTGG[ATC]TC
GGTTGGG[ACG]C
GGTTGGGT[AGT]
目的：將<IN>中的每行序列各自比對到<IN1>中所有基因序列上，將比對到的輸出，輸出位置信息及相應(yīng)的比對序列
我的程序：open IN1,"<seq_list.txt" ;
while (<IN1>
{
chomp;
my @aa1=split /\s+/,$_;
my $a=$aa1[0];
nn($a);
}
close IN1;
sub nn{
        my @a=$_;
open IN,"<test.fa" ;
open OUT,">>test_out.txt";
my %hash;
my $seq;
$/=">";
while (<IN>
{
chomp;
my @aa=split /\n/;
my $gene=(split /\s+/,shift @aa)[0];
$hash{$gene}=$_;
$seq=join('',@aa[0..$#aa]);
while ($seq =~ m/$a[0]/g)
{
my $pos = pos $seq;
push (@pos,$pos-;
pos($seq)=$pos-8;
if (pos($seq)){
print OUT "$gene\t", pos $seq, "\n";}
}
}
close IN;
}
出來的結(jié)果和我想要的不太一樣：AT1G50920.1        1        GGTTGGGTC
                                          AT1G50920.1        32        GGTTGGGTC
實際上還有結(jié)果沒有輸出

1002279625

額，沒有人會嗎

huang6894

回復(fù) 2# 1002279625

1. 
   
2. my($fa, $seq) = @ARGV;
   
3. 
   
4. my @all_seq;
   
5. open IN,"< $seq" or die "$!";
   
6. while(<IN>){
   
7.     chomp;
   
8.     if(/^\[([ATCG]+)\]([ATCG]+)$/){
   
9.         map{push @all_seq,$_.$2}(split //,$1);
   
10.     }elsif(/^([ATCG]+)\[([ATCG]+)\]([ATCG]+)$/){
    
11.         map{push @all_seq,$1.$_.$3}(split //,$2);
    
12.     }elsif(/^([ATCG]+)\[([ATCG]+)\]$/){
    
13.         map{push @all_seq,$1.$_}(split //,$2);
    
14.     }else{
    
15.         push @all_seq,$_;
    
16.     }
    
17. }
    
18. close IN;
    
19. $/ = '>';
    
20. open IN1,"< $fa" or die "$!";
    
21. <IN1>;
    
22. while(<IN1>){
    
23.     chomp;
    
24.     my ($key, @seq) = split /\n/;
    
25.     my ($gene) = $key =~/^([^\|]+)(\s+)?|/;
    
26.     my $s =  join("",@seq);
    
27.     foreach my $key2(@all_seq){
    
28.         my ($pos, $now) = (0, -1);
    
29.         until ( $pos eq -1 ) {
    
30.             $pos = index( $s, $key2, $now + 1 );
    
31.             $now = $pos;
    
32.             print $gene."\t".$pos."\t".$key2."\n" unless $pos < 0;
    
33.         }
    
34.     }
    
35. }
    
36. close IN2;
    
37. $/ = "\n";
    

復(fù)制代碼

結(jié)果：

1. AT1G50920.1     0       GGTTGGGTC
   
2. AT1G50920.1     31      GGTTGGGTC
   
3. AT1G15970.1     112     GCTTGGGTC
   
4. AT1G75120.1     15      CGTTGGGTC
   
5. AT1G75120.1     36      GGGTGGGTC
   

復(fù)制代碼

如果數(shù)據(jù)量小也就罷了，太大，可以考慮用kmp算法，或者直接blat

1002279625

回復(fù) 3# huang6894


    我是seq_list這個比較長，有116行，考慮到9bp中最多有兩個錯配。用過blat，但因為9bp太短，出不來結(jié)果，所以才自己寫程序比對。

1002279625

回復(fù) 3# huang6894
我發(fā)現(xiàn)，這樣寫，有漏掉的


   

chenyuluoyan沉魚落

回復(fù) 3# huang6894
留著以后慢慢研究大神程序~