索引一段序列

清泉一边

我有一段序列：

>1
AAAATAAATCACCTTTCGGGCTAACTTTGCGGTCGAGAACTCATTACCCAAATCCTACAATCACAAATTTG
ACATACTTAGAATTAAAAAACCAAACCACACACAGAAACAAGGTTAAGATAAATGAACAAAGGAGAATGATTTAGTTAGTAACCTCAACATTAGAGAGCT
TTCCTTCGTTTAATACTTTGAAGATGGCAAAACCACCTGGCGTCTCAAACAGTATTAGCATTTTTACAAGCTCCTGAAAAGAAGAAAAAACAAGATTAAG
AGAAACCCTTGAATCAACTCAAAGTCACCAAACTTGCAAGTTAGTGTTTTACGAGCTAAAGCCATTAAGCTAGTACGAAACAATATTCAAGGTAAAACTT
TTTCCCTGTTTCGAAGTTTACAAATCAGAATAATAGTTAAGGCAGATACTCTGTTCAAATTCTTGAAAATCCGACAAAGAACAGAACTATACTTTGTCAA
AGTCTTAGAATTTTGGAAGCTTATTTGCCATAACAACGAAGAAAGAGAGGGCGAGAGAGATTCAGTACCTTAGAGGGTTAACGCACAACTACAGCTTTAG
AGAGAGCTAAAATCAGGAGCCGCGATTCTGCTAGGGTTTAAGATGGTTTTTA
>2
TAAAAAAAAAATAACATCATTATATAATATATAGAGTTTAAAACATCTCAAAAACAAATTCATCATATTTTGTGATTCGAAATTTTAAGAATGAACATAT
ATTAACTAATTGGCGAAAAATGCGTGGGTTCAACGTCCCGCAACGAATAAAATATTTTGACAATGATTCATAAACATATTATAAATAAGATCAACATTAA
TAAAATAAATAATTTTTTTTTGTGGATGGATTTGGTTTGGCAGGACGTTACTTAATAACAATTGTAAACTATAAAATAATTTACAAATTTT
ATATATATTAATTTAAAAAATGAATTGTCTACGCGGTGTACCGCATGTTAAAATTTAGTTTCTATATATTTTAGAAACAACTTTGAATTTATACTTTAAT
ATTGAATAAACAACACCAAACCCCCTATTATTCATGTTATCCATTTTTTGAAATAACAGAAAAATAGAAAATAATCATA
AGAAACCAAACAAAATATACACAACAAAAAATCAAATCATAAAGCTTTAAATACATATAAGTGAAAGATCAAATCATAAAACTATAAAGACATGAAGTAC
CTAAAACATAATATATGCTAAAAAGAAATTCAAAATACAAAATCTTCTACGTATTTGAATAATTCATCCAAACCTAAAACTGTATATCTGTTCACATATT
TGAATGAAAAATCAAATAAGCAAATCAAGCAAGAAATATTAATGATGTCTATGTGATTTTTTTTTTTGT
>3
TTAATCTGCTTTTTTTTTTTTTTTTTAATTTACTCATATTAGATTTAGCTTAATTTTGAGACTGTTAGCTTTCGGTGTGAACAAAAGAAATTTGTGAAAT
TTGATATTGTTGATACATTCTCTAGAAATTTTGGAAAGATTGTGTGTTTCTTTTCAAAATTCAAATATTAATAACGCACCAAAATATCTGAATAGAAAGA
ATAAATAATGCGCCAAAATATTGATATGATGAAAGGTTCCCGTCTCAATATGTTTTTAGACCCTAGGTAAAACTAAAT
TTACATATCCTTTTCACACGATTTTTTTTTTTTTTTTTTGACTCTTTTACTTAAAGGTTTTTTTAAAAAAATTTGCCATGCACCCTGGCAATGGCTTTTG
CCCCCACCTCCCCCACATTAAGCCAATCTTGTTGCATGGCCATCTCCCCTGACGAACACCATTAAAGATTCATCTATATGTGGTA
GCCACTGAGTAGACTTAATAGAGCATTAAATAAATGAAATTCGTGGATGCAAATTGTAGAAGAACTAGTATTTAACGGAGTGTTGCTTCATCACAAATTC
>4......等等好多
我想从每一条中挑选出GCGT**六个碱基（两个星号代表GCGT后面两个模糊匹配）的位置，如果可以最好把对应的六个碱基序列也拿出来。
结果应该像这样：>1 GCGTAT 65(起始位置) GCGTTT 99(起始位置)
                      >2 GCGTGA  87(起始位置) GCGTAA 654(起始位置) GCGTTT 756(起始位置)
                      >3 ........
                      我写的起始位置不一定对，是我自己随意写的。
坐等大神！！！感激！！！

rubyish

biru:

1. 1>      GCGTAT 0, GCGTCT 216
   
2. 2>      GCGTGG 121
   
3. 3>
   

复制代码

1. #!/usr/bin/perl
   
2. $/ = ">";
   
3. <DATA>;
   
4. 
   
5. while (<DATA>) {
   
6.     chop if />$/;
   
7.     my @str = split /\n/;
   
8.     my $num = shift @str;
   
9.     my $str = join '', @str;
   
10.     my @index;
    
11.    
    
12.     while ( $str =~ /GCGT../g ) {
    
13.         push @index, "$& $-[0]";
    
14.     }
    
15. 
    
16.     print "$num>\t", join ', ', @index;
    
17.     print "\n";
    
18. }
    
19. 
    
20. 
    
21. __DATA__
    
22. >1
    
23. GCGTAT
    
24. AAAATAAATCACCTTTCGGGCTAACTTTGCGGTCGAGAACTCATTACCCAAATCCTACAATCACAAATTTG
    
25. ACATACTTAGAATTAAAAAACCAAACCACACACAGAAACAAGGTTAAGATAAATGAACAAAGGAGAATGATTTAGTTAGTAACCTCAACATTAGAGAGCT
    
26. TTCCTTCGTTTAATACTTTGAAGATGGCAAAACCACCTGGCGTCTCAAACAGTATTAGCATTTTTACAAGCTCCTGAAAAGAAGAAAAAACAAGATTAAG
    
27. AGAAACCCTTGAATCAACTCAAAGTCACCAAACTTGCAAGTTAGTGTTTTACGAGCTAAAGCCATTAAGCTAGTACGAAACAATATTCAAGGTAAAACTT
    
28. TTTCCCTGTTTCGAAGTTTACAAATCAGAATAATAGTTAAGGCAGATACTCTGTTCAAATTCTTGAAAATCCGACAAAGAACAGAACTATACTTTGTCAA
    
29. AGTCTTAGAATTTTGGAAGCTTATTTGCCATAACAACGAAGAAAGAGAGGGCGAGAGAGATTCAGTACCTTAGAGGGTTAACGCACAACTACAGCTTTAG
    
30. AGAGAGCTAAAATCAGGAGCCGCGATTCTGCTAGGGTTTAAGATGGTTTTTA
    
31. >2
    
32. TAAAAAAAAAATAACATCATTATATAATATATAGAGTTTAAAACATCTCAAAAACAAATTCATCATATTTTGTGATTCGAAATTTTAAGAATGAACATAT
    
33. ATTAACTAATTGGCGAAAAATGCGTGGGTTCAACGTCCCGCAACGAATAAAATATTTTGACAATGATTCATAAACATATTATAAATAAGATCAACATTAA
    
34. TAAAATAAATAATTTTTTTTTGTGGATGGATTTGGTTTGGCAGGACGTTACTTAATAACAATTGTAAACTATAAAATAATTTACAAATTTT
    
35. ATATATATTAATTTAAAAAATGAATTGTCTACGCGGTGTACCGCATGTTAAAATTTAGTTTCTATATATTTTAGAAACAACTTTGAATTTATACTTTAAT
    
36. ATTGAATAAACAACACCAAACCCCCTATTATTCATGTTATCCATTTTTTGAAATAACAGAAAAATAGAAAATAATCATA
    
37. AGAAACCAAACAAAATATACACAACAAAAAATCAAATCATAAAGCTTTAAATACATATAAGTGAAAGATCAAATCATAAAACTATAAAGACATGAAGTAC
    
38. CTAAAACATAATATATGCTAAAAAGAAATTCAAAATACAAAATCTTCTACGTATTTGAATAATTCATCCAAACCTAAAACTGTATATCTGTTCACATATT
    
39. TGAATGAAAAATCAAATAAGCAAATCAAGCAAGAAATATTAATGATGTCTATGTGATTTTTTTTTTTGT
    
40. >3
    
41. TTAATCTGCTTTTTTTTTTTTTTTTTAATTTACTCATATTAGATTTAGCTTAATTTTGAGACTGTTAGCTTTCGGTGTGAACAAAAGAAATTTGTGAAAT
    
42. TTGATATTGTTGATACATTCTCTAGAAATTTTGGAAAGATTGTGTGTTTCTTTTCAAAATTCAAATATTAATAACGCACCAAAATATCTGAATAGAAAGA
    
43. ATAAATAATGCGCCAAAATATTGATATGATGAAAGGTTCCCGTCTCAATATGTTTTTAGACCCTAGGTAAAACTAAAT
    
44. TTACATATCCTTTTCACACGATTTTTTTTTTTTTTTTTTGACTCTTTTACTTAAAGGTTTTTTTAAAAAAATTTGCCATGCACCCTGGCAATGGCTTTTG
    
45. CCCCCACCTCCCCCACATTAAGCCAATCTTGTTGCATGGCCATCTCCCCTGACGAACACCATTAAAGATTCATCTATATGTGGTA
    
46. GCCACTGAGTAGACTTAATAGAGCATTAAATAAATGAAATTCGTGGATGCAAATTGTAGAAGAACTAGTATTTAACGGAGTGTTGCTTCATCACAAATTC
    

复制代码

rubyish

or overlap - 交叠
>1
GCGTGCGTAA

1. 1>      GCGTGC 0, GCGTAA 4

复制代码

result:

1. 1>      GCGTGC 0, GCGTAA 4, GCGTCT 220
   
2. 2>      GCGTGG 121
   
3. 3>
   

复制代码

1. #!/usr/bin/perl
   
2. 
   
3. $/ = ">";
   
4. <DATA>;
   
5. 
   
6. while (<DATA>) {
   
7.     chop if />$/;
   
8.     my @str = split /\n/;
   
9.     my $num = shift @str;
   
10.     my $str = join '', @str;
    
11.     my @index;
    
12.     $str =~ /GCGT..(?{ push @index, "$& $-[0]" })(*F)/;
    
13.     print "$num>\t", join ', ', @index;
    
14.     print "\n";
    
15. }
    
16. 
    
17. __DATA__
    
18. >1
    
19. GCGTGCGTAA
    
20. AAAATAAATCACCTTTCGGGCTAACTTTGCGGTCGAGAACTCATTACCCAAATCCTACAATCACAAATTTG
    
21. ACATACTTAGAATTAAAAAACCAAACCACACACAGAAACAAGGTTAAGATAAATGAACAAAGGAGAATGATTTAGTTAGTAACCTCAACATTAGAGAGCT
    
22. TTCCTTCGTTTAATACTTTGAAGATGGCAAAACCACCTGGCGTCTCAAACAGTATTAGCATTTTTACAAGCTCCTGAAAAGAAGAAAAAACAAGATTAAG
    
23. AGAAACCCTTGAATCAACTCAAAGTCACCAAACTTGCAAGTTAGTGTTTTACGAGCTAAAGCCATTAAGCTAGTACGAAACAATATTCAAGGTAAAACTT
    
24. TTTCCCTGTTTCGAAGTTTACAAATCAGAATAATAGTTAAGGCAGATACTCTGTTCAAATTCTTGAAAATCCGACAAAGAACAGAACTATACTTTGTCAA
    
25. AGTCTTAGAATTTTGGAAGCTTATTTGCCATAACAACGAAGAAAGAGAGGGCGAGAGAGATTCAGTACCTTAGAGGGTTAACGCACAACTACAGCTTTAG
    
26. AGAGAGCTAAAATCAGGAGCCGCGATTCTGCTAGGGTTTAAGATGGTTTTTA
    
27. >2
    
28. TAAAAAAAAAATAACATCATTATATAATATATAGAGTTTAAAACATCTCAAAAACAAATTCATCATATTTTGTGATTCGAAATTTTAAGAATGAACATAT
    
29. ATTAACTAATTGGCGAAAAATGCGTGGGTTCAACGTCCCGCAACGAATAAAATATTTTGACAATGATTCATAAACATATTATAAATAAGATCAACATTAA
    
30. TAAAATAAATAATTTTTTTTTGTGGATGGATTTGGTTTGGCAGGACGTTACTTAATAACAATTGTAAACTATAAAATAATTTACAAATTTT
    
31. ATATATATTAATTTAAAAAATGAATTGTCTACGCGGTGTACCGCATGTTAAAATTTAGTTTCTATATATTTTAGAAACAACTTTGAATTTATACTTTAAT
    
32. ATTGAATAAACAACACCAAACCCCCTATTATTCATGTTATCCATTTTTTGAAATAACAGAAAAATAGAAAATAATCATA
    
33. AGAAACCAAACAAAATATACACAACAAAAAATCAAATCATAAAGCTTTAAATACATATAAGTGAAAGATCAAATCATAAAACTATAAAGACATGAAGTAC
    
34. CTAAAACATAATATATGCTAAAAAGAAATTCAAAATACAAAATCTTCTACGTATTTGAATAATTCATCCAAACCTAAAACTGTATATCTGTTCACATATT
    
35. TGAATGAAAAATCAAATAAGCAAATCAAGCAAGAAATATTAATGATGTCTATGTGATTTTTTTTTTTGT
    
36. >3
    
37. TTAATCTGCTTTTTTTTTTTTTTTTTAATTTACTCATATTAGATTTAGCTTAATTTTGAGACTGTTAGCTTTCGGTGTGAACAAAAGAAATTTGTGAAAT
    
38. TTGATATTGTTGATACATTCTCTAGAAATTTTGGAAAGATTGTGTGTTTCTTTTCAAAATTCAAATATTAATAACGCACCAAAATATCTGAATAGAAAGA
    
39. ATAAATAATGCGCCAAAATATTGATATGATGAAAGGTTCCCGTCTCAATATGTTTTTAGACCCTAGGTAAAACTAAAT
    
40. TTACATATCCTTTTCACACGATTTTTTTTTTTTTTTTTTGACTCTTTTACTTAAAGGTTTTTTTAAAAAAATTTGCCATGCACCCTGGCAATGGCTTTTG
    
41. CCCCCACCTCCCCCACATTAAGCCAATCTTGTTGCATGGCCATCTCCCCTGACGAACACCATTAAAGATTCATCTATATGTGGTA
    
42. GCCACTGAGTAGACTTAATAGAGCATTAAATAAATGAAATTCGTGGATGCAAATTGTAGAAGAACTAGTATTTAACGGAGTGTTGCTTCATCACAAATTC
    

复制代码

yestreenstars

use 5.010;
my ( %h, $k, $l );
while (<DATA>) {
        chomp;
        if (/>/) {
                $l = 0;
                $k = $_;
                $h{$_} = $_;
                next;
        }
        while (/GCGT../) {
                my $i = index( $_, $& );
                $h{$k} .= " $& " . ( $l + $i );
                s/GCGT/    /;
        }
        $l += length;
}
say $h{$_} for sort keys %h;

rubyish

回复 4# yestreenstars

maybe small bug~

1. __DATA__
   
2. >1
   
3. GCGTGCGTAAACGCGT
   
4. AAAATAAATCACGCGTAAAAGCGTA
   
5. >2
   
6. AAAATAAATCACGCGTAA
   
7. AAATCAC
   
8. >3
   
9. TAAATCACG

复制代码

1. 1>      GCGTGC 0, GCGTAA 4, GCGTAA 12, GCGTAA 28
   
2. 2>      GCGTAA 12
   
3. 3>

复制代码

1. >1 GCGTGC 0 GCGTAA 4 GCGTAA 28
   
2. >2 GCGTAA 12
   
3. >3

复制代码

yestreenstars

回复 5# rubyish

So lucky to reply you, I got a book badges.

Thanks for debug. I have modify my script.

The script after modify:

use 5.010; my ( %h, $k, $s ); while (<DATA>) {         chomp;         if (/>/) {                 undef $s;                 $k = $_;                 $h{$_} = $_;                 next;         }         $s .= $_;         while ( $s =~ /GCGT../ ) {                 my $i = index( $s, $& );                 $h{$k} .= " $& $i";                 $s =~ s/GCGT/    /;         } } say $h{$_} for sort keys %h;

清泉一边

学习了！感谢！回复 2# rubyish