regex - 如何制作 FASTA 文件中特定模式出现的直方图？

Question

我为以下生物信息学问题编写了一个 Perl 脚本，但不幸的是，输出存在问题。

问题

1）从40,000个唯一序列的文件中，唯一的意思是序列ID号，提取以下模式

 $gpat = [G]{3,5}; $npat = [A-Z]{1,25};<br>
 $pattern = $gpat.$npat.$gpat.$npat.$gpat.$npat.$gpat;  

2）对于每个序列，查找是否$pattern出现在的值之间

• 0-100
• 100-200
• 200-300
• ...
• 900-1000

• 1000

如果某个序列的长度小于 1000 个字符，则即使如此，也必须保持除法，即 0-100,100-200 等。

问题

我遇到的主要问题是计算每个序列细分出现 $pattern 的次数，然后为所有序列添加它的计数。

例如，对于序列 1，假设 $pattern 在长度 >1000 处出现 5 次。对于序列 2，假设 $pattern 在长度>1000 时出现 3 次。那么总数应该是 5+3 =8。

相反，我的结果是这样的： (5+4+3+2+1) + (3+2+1) = 21 即累计总数。

对于每个 100 个字符的前 10 个细分，我面临同样的问题。

如果可以为此计算提供正确的代码，我将不胜感激。

我写的代码如下。它很大程度上来自 Borodin 对我之前的一个问题的回答：Perl: Search a pattern across array elements

他的答案在这里：https ://stackoverflow.com/a/11206399/1468737

代码：

use strict;
use warnings;

my $gpat = '[G]{3,5}';
my $npat = '[A-Z]{1,25}';
my $pattern = $gpat.$npat.$gpat.$npat.$gpat.$npat.$gpat; 
my $regex = qr/$pattern/i;

open my $fh, '<', 'small.fa' or die $!;

my ($id, $seq); 
my @totals = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0); #intialize the @total arrays...  
#..it should  contain 10 parts for 10 divisions upto 1000bp
my @thousandcounts =(0); #counting total occurrences of $pattern at >1000 length

while (<$fh>) {
  chomp;

  if (/^>(\w+)/) {
    process_seq($seq) if $id;
    $id = $1;
    $seq = '';
    print "$id\n";
  }
  elsif ($id) {
    $seq .= $_;
    process_seq($seq) if eof;
  }
}

print "Totals : @totals\n";
print "Thousand Counts total : @thousandcounts\n";

##**SUBROUTINE**    

sub process_seq {

  my $sequence = shift @_;   
  my $subseq = substr $sequence,0,1000;
  my $length = length $subseq;
  print $length,"\n";

  if ($length eq 1000) {

  my @offsets = map {sprintf '%.0f', $length * $_/ 10} 1..10;
  print "Offsets of 10 divisions: @offsets\n";

  my @counts = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
  my @count = (0); 

     while ($sequence =~ /$regex/g) {
     my $place = $-[0];
     print $place,"\n\n"; 

        if ($place <=1000){
        for my $i (0..9) { 
        next if $place >= $offsets[$i];                   
        $counts[$i]++;                                    
        last;
        }       

     }
      print "Counts : @counts\n\n";

      $totals[$_] += $counts[$_] for 0..9; 



        if ($place >1000){

        for my $i(0){
        $count[$i]++;
        last;
        }




        } print "Count greater than 1000 : @count\n\n"; 

         $thousandcounts[$_] += $count[$_] for 0;


  } 

} 

   #This region of code is for those sequences whose total length is less than 1000
   #It is working great ! No issues here
   elsif ($length != 1000) {

    my $substr = join ' ', unpack '(A100)*', $sequence;

    my @offsets = map {sprintf '%.0f', $length * $_/ ($length/100)} 1..10;
    print "Offsets of 10 divisions: @offsets\n";

    my @counts = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0,);

       while ($sequence =~ /$regex/g) {
       my $place = $-[0];
       print "Place : $place","\n\n"; 

         for my $i (0..9) { 
         next if $place >= $offsets[$i];                   
         $counts[$i]++;                                    .
         last;
        }
      }
       print "Counts : @counts\n\n";

       $totals[$_] += $counts[$_] for 0..9;

  }


}#subroutine ends

我还附上了我正在使用的文件的一小部分。这个文件的标题是small.fa，我一直在尝试使用这个文件，然后才转到包含超过 40,000 个序列的更大文件。

>NR_037701 1
aggagctatgaatattaatgaaagtggtcctgatgcatgcatattaaaca
tgcatcttacatatgacacatgttcaccttggggtggagacttaatattt
aaatattgcaatcaggccctatacatcaaaaggtctattcaggacatgaa
ggcactcaagtatgcaatctctgtaaacccgctagaaccagtcatggtcg
gtgggctccttaccaggagaaaattaccgaaatcactcttgtccaatcaa
agctgtagttatggctggtggagttcagttagtcagcatctggtggagct
gcaagtgttttagtattgtttatttagaggccagtgcttatttagctgct
agagaaaaggaaaacttgtggcagttagaacatagtttattcttttaagt
gtagggctgcatgacttaacccttgtttggcatggccttaggtcctgttt
gtaatttggtatcttgttgccacaaagagtgtgtttggtcagtcttatga
cctctattttgacattaatgctggttggttgtgtctaaaccataaaaggg
aggggagtataatgaggtgtgtctgacctcttgtcctgtcatggctggga
actcagtttctaaggtttttctggggtcctctttgccaagagcgtttcta
ttcagttggtggaggggacttaggattttatttttagtttgcagccaggg
tcagtacatttcagtcacccccgcccagccctcctgatcctcctgtcatt
cctcacatcctgtcattgtcagagattttacagatatagagctgaatcat
ttcctgccatctcttttaacacacaggcctcccagatctttctaacccag
gacctacttggaaaggcatgctgggtctcttccacagactttaagctctc
cctacaccagaatttaggtgagtgctttgaggacatgaagctattcctcc
caccaccagtagccttgggctggcccacgccaactgtggagctggagcgg
gagggaggagtacagacatggaattttaattctgtaatccagggcttcag
ttatgtacaacatccatgccatttgatgattccaccactccttttccatc
tcccagaagcctgctttttaatgcccgcttaatattatcagagccgagcc
tggaatcaaactgcctctttcaaaacctgccactatatcctggctttgtg
acctcagccaagttgcttgactattctcagtctcagtttctgcacctgtc
aaatagggtttatgttaacctaactttcagggctgtcaggattaaatgag
catgaaccacataaaatgtttggtgtatagtaagtgtacagtaaatactt
ccattatcagtccctgcaattctatttttcttccttctctacacagcccc
tgtctggctttaaaatgtcctgccctgctttttatgagtggataccccca
gccctatgtggattagcaagttaagtaatgacactcagagacagttccat
ctttgtccataacttgctctgtgatccagtgtgcatcactcaaacagact
atctcttttctcctacaaaacagacagctgcctctcagataatgttgggg
gcataggaggaatgggaagcccgctaagagaacagaagtcaaaaacagtt
gggttctagatgggaggaggtgtgcgtgcacatgtatgtttgtgtttcag
gtcttggaatctcagcaggtcagtcacattgcagtgtgtcgcttcacctg
gctccctcttttaaagattttccttccctctttccaactccctgggtcct
ggatcctccaacagtgtcagggttagatgccttttatgggccacttgcat
tagtgtcctgatagaggcttaatcactgctcagaaactgccttctgccca
ctggcaaagggaggcaggggaaatacatgattctaattaatggtccaggc
agagaggacactcagaatttcaggactgaagagtatacatgtgtgtgatg
gtaaatgggcaaaaatcatcccttggcttctcatgcataatgcatgggca
cacagactcaaaccctctctcacacacatacacatatacattgttattcc
acacacaaggcataatcccagtgtccagtgcacatgcatacacgcacaca
ttcccttcctaggccactgtattgctttcctagggcatcttcttataaga
caccagtcgtataaggagcccaccccactcatctgagcttatcaaccaat
tacattaggaaagactgtatttcctagtaaggtcacattcagtagtactg
agggttgggacttcaacacagctttttgggggatcataattcaacccatg
acagccactgagattattatatctccagagaataaatgtgtggagttaaa
aggaagatacatgtggtacaaggggtggtaaggcaagggtaaaaggggag
ggaggggattgaactagacacagacacatgagcaggactttggggagtgt
gttttatatctgtcagatgcctagaacagcacctgaaatatgggactcaa
tcattttagtccccttctttctataagtgtgtgtgtgcggatatgtgtgc
tagatgttcttgctgtgttaggaggtgataaacatttgtccatgttatat
aggtggaaagggtcagactactaaattgtgaagacatcatctgtctgcat
ttattgagaatgtgaatatgaaacaagctgcaagtattctataaatgttc
actgttattagatattgtatgtctttgtgtccttttattcatgaattctt
gcacattatgaagaaagagtccatgtggtcagtgtcttacccggtgtagg
gtaaatgcacctgatagcaataacttaagcacacctttataatgacccta
tatggcagatgctcctgaatgtgtgtttcgagctagaaaatccgggagtg
gccaatcggagattcgtttcttatctataatagacatctgagcccctggc
ccatcccatgaaacccaggctgtagagaggattgaggccttaagttttgg
gttaaatgacagttgccaggtgtcgctcattagggaaaggggttaagtga
aaatgctgtataaactgcatgatgtttgcaggcagttgtggttttcctgc
ccagcctgccaccaccgggccatgcggatatgttgtccagcccaacacca
caggaccatttctgtatgtaagacaattctatccagcccgccacctctgg
actccctcccctgtatgtaagccctcaataaaaccccacgtctcttttgc
tggcaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
aaa
>NR_002714 1
gttatacatctctaccattacctagcctgaaaagccacctcagattcagc
caacaagtaagtgggcattacaggagaagggtacctttcacaagggctgt
aatctaaaatcttggggaagatacagcgtcatctgtccaagaggtgtcag
cagtaacgaagcctcagtagaagccaaagttattttggattactgagcct
gtatagtttccagattctcaagagaaatatatgggaatgtagatatctca
gaggaccttcctgctgtcaggaattcagaggaggaaataaggaaggtaat
aggtgctctgctctcattctctcaaaccctcttccctgtgttttcctata
gagattgctgatttgctccttaagcaagagattcactgctgctcagcatg
gctcagaccaactcatgcttcatgctgatctcctgcctgatgttcctgtc
tctgagccaaggtgagattgttttccccacacatacctcccacaacccca
gccctgaagccctcactctatcctcatgcatatgagttcacttgagaaaa
agcagagtcaagttcaggggttgttttgtgttgttcagtgatatttattg
ctgatctcatcccattcaaaaacatcctgacctccctaaggagttagaga
tggaacttagcataaccctttatcagtgaccactgcagttggcattggtt
tgtcatattaacactactcatgatgggggtgttgaggatgtctgtttgta
gacagtcattagtggaatggggaactgaggggagctttgtgtgtagagaa
actggacaggcttgagaaagaagcctcagtccttcaaggaagaaaaagcc
ataagtaaaagggacaatggggacacttttcatgagcctattcattgtgt
gctcttgtcttgagcaaagacatcttgagagcctataggtaagatgcaga
agggcagaagtgaccaatcgcttcgtgacctataggatccttctattcct
ataaagaatcctcagaagctcctacctcatattttagcctttaccttgcc
ctgagggtctttcttaattgtctctcttttcccaggacaggaggcccatg
ctgagttgcccaaggcccagatcagctgcccagaaggcaccagtgcctaa
ggctcccactgctactactttaatgaagagcatgagacctgggtttatgc
agatgtgagtgaggagagcagtgtgggaagggaggctcacgaagggaggg
gaagctgccactctccagtgtgttcagtggctgatatgagatgagactaa
tcccctccctatccaatcatcagcccaaaactttccaatctactttatcc
catcattcagcacagagatgctggtggtcagtgacagcatcatcagggac
atttctgtgctgtcctttttctgttacatcctctgggagggctcaatatg
tctcccacactttcctccttcactgagtgctccattttcttctccaacag
ctctactgccagaacatgaattcaggtaacctggtgtctgtgctcaccca
ggctgagggtgcctttgtggcttcgctgattaaagagagtggcaccaagg
atagcaatgtctggattggcctccatgacccccaccggatcagtctgctg
catcttctacctcctgattatcaggttccagagggtctgatgtctggcac
ctcaagcatcagtttttactatattatgataaaagcaacctctctataaa
tcatataatgtaaaggatatcaaggttctccataggttcttcgagataag
cttaaagctgaatttcctgtgtgtttcaggcattcacagataaactcatt
ctctgtacttctagggtagcatctttatgtatctattatgtacctcttat
ctattgtgttatcatctctgttatagaagagccttctgtagaccatatag
aaaaagattatagaggaggagaatctactgctggcaattgggaaccgcaa
ggtatactaaataatatatcaacaactaatggccatctaatgctatgctg
gatatgaacttttggggcctcaggaaagaaaaaccaggaactagtttcaa
taatgaggtgtcatggttccctgtggcaaatttagaacgcttatcgtttg
gcaggacacagagaggtaggtgaacattccaggaaagaagcagcttagag
aaaatgtggaggaaataatatgacacttagagaaaaaggaaggtttattc
ttgtcttatgtcttgacctgtttctgagtgcgaacacaaaccaggtgttt
ctgtctctttctgagtcacgtctgcccctgttctggcccttccccatcta
gaactgccattatcagtggagtagtgggtccctggtctcctacaaatcct
gggacattggatccccaagctgtgccaatactgcctactgtgctagcctg
acttcaagctcaggtgaggggcacagaatccacacacttattgccatcct
ctcctatttatctctgaggatcgaccggggactgggatagaggaagggtg
agctcctcattcaggaaatagaggagtgtttcctctttatttttgctgag
tcctgcagccaggagggtaatacactctgatcccctcagtctgaatcttc
tcattgtcttataggattcaagaaatggaaggatgattcttgtaaggaga
agttctcctttgtttgcaagttcaaatactggaggcaattgtaaaatgga
cgtctagaattggtctaccagttactatggagtaaaagaattaaactgga
ccatctctctccatatcaatctggaccatctctcctctgctaaatttgca
tgactgatctttagtatctttacctacctcaatttctggagccctaaaca
ataaaaataaacatgtttcccccat
>NR_003569 1
ctgggacccacgacgacagaaggcgccgatggccgcgcctgctgagccct
gcgcggggcagggggtctggaaccagacagagcctgaacctgccgccacc
agcctgctgagcctgtgcttcctgagaacagcaggggtctgggtaccccc
catgtacctctgggtccttggtcccatctacctcctcttcatccaccacc
atggccggggctacctccggatgttccccactcttcaaagccaagatggt
gcttggattcgccctcatagtcctgtgtacctccagcgtggctgtcgctc
tttggaaaatccaacagggaacgcctgaggccccagaattcctcattcat
cctactgtgtggctcaccacgatgagcttcgcagtgttcctgattcacac
caagaggaaaaagggagtccagtcatctggagtgctgtttggttactggc
ttctctgctttgtcttgccagctaccaacgctgcccagcaggcctccgga
gcgggcttccagagcgaccctgtccgccacctgtccacctacctatgcct
gtctctggtggtggcacagtttgtgctgtcctgcctggcggatcaacccc
ccttcttccctgaagacccccagcagtctaacccctgtccagagactggg
gcagccttcccctccaaagccacgttctggtgggtttctggcctggtctg
gaggggatacaggaggccactgagaccaaaagacctctggtcgcttggga
gagaaaactcctcagaagaacttgtttcccggcttgaaaaggagtggatg
aggaaccgcagtgcagcccgggggcacaacaaggcaatagcatttaaaag
gaaaggcggcagtggcatggaggctccagagactgagcccttcctacggc
aagaagggagccagtggcgcccactgctgaaggccatctggcaggtgttc
cattctaccttcctcctggggaccctcagcctcgtcatcagtgatgtctt
caggttcactgtccccaagctgctcagccttttcctggagtttattggtg
atcccaagcctccagcctggaagggctacctcctcgccgtgctgatgttc
ctctcggcctgcctgcaaacgctgtttgagcagcagaacatgtacaggct
caaggtgctgtagatgaggctgcggtcggccatcactggcctggtgtaca
gaaaggcatccacagcatatctgaagaaatattcagaagttaactaatct
cagatgatttcagcaggagtaaagaagagaaacagactcagaaatgccat
tacaacagttaattatgtcaaatttatcaccctgattgatcacgcagcat
taacctcaagaacgccaagccaagtttttttgacaaatgtgagccaaggt
ttccgaaaaactagcagatatgactgtgacttacaaaatggaaaaagtaa
acgagaaacacaatttgatatgatttaataaaagatttgtttccaccact
tctcctgggaacctcagcacattttctttccactgacagttattatctct
acctttattgaacaaagacacccggaacacagctgctgaggatcagtaaa
gaaaatcattcttttattaataagactgttattagcaggaaaaaaaaatc
catgtttgggagtttgcactgaagttacaggccattttgaagaaatatgg
ctgactagtgccaacattatttcaggcaatttcatgatcaaatgtcttat
taggttgtttaaaatttttatagagattgtaaatcagaactattttctat
ttgccctaaatatttagatgctacagggaaagcagatcaaattaaagggt
actgtgcacatttttttactgggaactcccagggatataaatcatttcgc
ctgcagcatggaattcttcagtacacatgcttgtggaaacattccacgct
ccgccagcacgctcattaaagtgatgatttgggttgcaacaacagtgcca
agtacttcctgtgttcaactggggaccatgtggcaagacccaaagcttcc
ccagagatcctatgggaataagttttttgagccaccatattccattattt
cagcctaaaataacaccatgggacaagaatcagaagacagaggagcagac
aaatgtgtgtagacatgctggaaggaatctttctttttagaaacagggtc
aatatctattaaactttaagatgtgtatctcttgacctggcagtttctgt
atttgagttttaacctactgatatacccatgcatgtgaataaagtatctt
cctgcatgtaacaggatatttaatgtaaccttgattatagttgcaaatgc
tgggaaacgatccaaatgtctttcaatatggcactgattaaataaattat
ggcacagtctcacaatgaaaaacaaatgtagccattaaacagaatgaaat
gggtctagctaaattgaaataggactacctctaagatatgttgttaaaaa
gaaaaaaaagaaagtgcagaggaacaagtatgataccattttgtattttt
taacatatgcaagcgtgattgtgcccacacagaatacctttgaaaataaa
ctcagtatttgcctcagtggataaaaacaagaaccagccttattttcact
gttatatcttttggtgccactttttgaactttttaccatatgtgcatatg
taactttctaaataaattttgtaaaaaaaaaaaaaaaaaa
>NR_002817 2
aactcggtctccactgcactgctggccagacgagggatgttattttgggc
agtgcatctggacttggttcaagtggcaccagccaaatccctgccttact
gacctctcccctggaggagcaggagcagtgctcaaggccgccctgggagg
gctgagaggcaggctctggactggggacacagggatagctgagccccagc
tgggggtggaagctgagccagggacagtcacagaggaacaagatcaagat
gcgctttaactgagaagcccccaaggcagaggctgagaatcagaagacat
ttcagcagacatctacaaatctgaaggacaaaacatggttcaagcatctg
ggcacaggcggtccacccgtggctccaaaatggtctcctggtccgtgata
gcaaagatccaggaaatatggtgcgaggaagatgagaggaagatggcgcg
agagttcctggccgagttcatgagcacatatgtcatgatggagtggctga
ccgggatgctccagctgtgtctcttcgccatcgtggaccaggagaacaac
ccagcactgccaggaacacacgcactggtgataggcatcctcgtggtcat
catcagggtgtaccatggcatgaacacaggatatgccatcaatccgtccc
gggacctgccccccccccccgcatcttcaccttcattgctggttggggca
aactggtcttcaggtactgcccctgcccaggcccattcctttgagatttt
ctgtggggcccctgtgtgttgaggtgtggggggtgatgtgaggggcagca
caggagggtcctgcagagcccccaggtggcctggggagcaggagtgagtc
ccaacatttccccaggccagtagagatacagatcctgcacctgcactgag
tgtcaaccctgtccctgagtcgggctgaggctgaccagggccccgggttg
ggggtgtttcctgggttagcctgaggatgactcctctgctcaaccagtct
tggcccgaggtggatgagggtgctgtcctgggcatcagccccctcagccg
gcctctgcctcttgcctgcagcgatggggagaacttgtggtgggtgccag
tggtggcaccacttctgggtgcctctctaggtggcatcatctacctggtc
ttcattggctccaccatcccacgggagcccctgaaattggaggactctgt
ggcatatgaagaccacgggataaccgtattgcccaagatgggatctcatg
aacccatgatctctccccttaccctcatctccgtgagccctgccaacaga
tcttcagtccaccctgccccacccttacatgaatccatggccctagagca
cttctaagcagagattatttgtgatcccatcccttccccaataaagagaa
gcttgtcccacagcagtacccccacttcctgggggcctcctgtggttggg
cttccctcctgggttcttccaggagctctagggctatgtcttagcccaag
gtgtagaggtgaggcacctcaagtctttcatgccctgggaactggggtgc
cccagggggagaatggggaagagctgacctgcgccctcagtaggaacaag
gtaagatgaaagaatgacagaaacagaatgagggattttcaggcaagggg
gaaggaagggcagttttggtgaaaggactgtagctgactggtggggggct
ggctttggaaatactttgaggggatcctgagactggactctagactctcc
cctggttgttcccttccccgagttctggccggttcttggaccagacaagg
catggcccaagaaggtagatcagaattttttagcctttttttcattagtg
ccttccctagtataattccagattttttttcttaatcacatgaaatttta
ataccacagatatactatacatctgtttatgttctgtatatgttctgtgc
tttatacgtaaaaaagagtaagattttttttcacctccccttttaagaat
cagttttaattcccttgagaatgcttgttatagattgaaggctggtaagg
ggttgggctcctctttcttcttcctggtgccagagtgctcccacatgaag
gaataggaaaggaagatgcaaagagggaaatccttcgaacacatgaagac
acaggaagaggcctcttagggctccaagggctccagggaagcagctgcag
aggttgggtggggtgaggggccaggatccactgaccctggggccaggcag
gaatcactctgttgcctggggctcagaaggcagtatcacccatggttcct
gtcattgctcatgtattttgcctttcaacaattattgtgcacctactgtg
tgcaggccctgcctggacactggggatgcgcagtggatgcactgggctct
gcctttgagggttgcagtttaatgggtgacaggtaattataaggaagaag
gtgagtgcagagtgggaggcttggaggctgtggggcttggggtgggggag
ctcacatccagcctctgggccaaggccaggaggcttcccagagcaggaga
cagagcagggtattgtggtggggggtgtcctttttggggctgggatctgc
actttacagtttgaggggatgggcagaggaggctgggcttcattctggag
gtggggacatggtgaggtgaggtttagaaagcacacctgagccgcagtgt
gtaggatgctggaaatggtggagatgggcctgcgaagagagtgctgggaa
gtgatgacccaggagcagcagccgggcacctaacaatgggtcagcaccgt
gggcgtggagacaaaggccgggattgatcaatacccgagaagtacaatgt
acaggacttgggctccatttggatggagtgggtgagggaggagtcagaaa
tggcttccgatttccagcttgggcctggggattggagatgtccccactga
gagtagggcacaagtgaggaaatggtttggagaggaagatgataagttac
atcatggatgtgctgagtctgagttgcctatgggacttggaatggggggt
ggcaaaaggtgtgtgatcttgagcaagatattcaactcttctgggccttg
gtcttctcatttgtaaaacggtgataagaatattacttcccatttgtgtt
gctgtgaatattaaatgcgctaccacatgt

感谢您花时间解决我的问题。

任何帮助和输入将不胜感激。

感谢您花时间解决我的问题！

Answer 1

以下似乎有效。在玩的过程中，我将您发布在__DATA__脚本末尾的部分中的数据放在了这里。要将它与真实数据文件一起使用，您需要打开它，并将文件句柄传递给run.

#!/usr/bin/env perl

use strict; use warnings;
use Data::Dumper;
use List::MoreUtils qw( first_index );

if (@ARGV) {
    my ($input_file) = @ARGV;
    open my $input, '<', $input_file
        or die "Cannot open '$input_file': $!";
    run($input);
    close $input
        or die "Cannot close '$input_file': $!";
}
else {
    run(\*DATA);
}

sub run {
    my ($fh, $start_pat, $stop_pat) = @_;

    # These are your patterns. I changed $npat because I don't
    # think, e.g., q is a valid character in your input.
    my $gpat = '[g]{3,5}';
    my $npat = '[acgt]{1,25}';
    my $wanted = qr/$gpat$npat$gpat$npat$gpat$npat$gpat/;

    # These just tell us where a sequence begins and ends.
    my $start = qr/\A>([A-Za-z_0-9]+)/;
    my $stop = qr/[^acgt]/;

    # Set up the bins and labels for the histogram.
    my @bins = map 100 * $_, 1 .. 10;
    my @labels = map sprintf('%d - %d', $_ - 100, $_), @bins;

    # Initialize the histogram with all zero counts.
    my %hist = map { $_ => 0 } @labels;

    my $id;
    while (my $line = <$fh>) {
        # Whenever you see a new sequence, read it completely
        # and pass it to build_histogram.
        if (($id) = ($line =~ $start)) {
            print "Start sequence: '$id':\n";
            my $seq_ref;
            ($line, $seq_ref) = read_sequence($fh, $stop);

            my $hist = build_histogram(
                $seq_ref,
                $wanted,
                \@bins,
                \@labels,
            );

            # Add the counts from this sequence to the overall
            # histogram.

            for my $key ( keys %$hist ) {
                $hist{ $key } += $hist->{$key};
            }

            # exit loop if read_sequence stopped because of EOF.
            last unless defined $line;

            # else see if the line that stopped input is the start
            # of a new sequence.
            redo;
        }
    }

    print Dumper \%hist;
}

sub build_histogram {
    my ($seq_ref, $wanted, $bins, $labels) = @_;

    my %hist;

    while ($$seq_ref =~ /$wanted/g) {
        # Whenever we find segment which matches what we want,
        # store the position,
        my $pos = $-[0];

        # and find the bin where it fits.
        my $idx = first_index { $_ > $pos } @$bins;

        # if you do not have List::MoreUtils, you should install it
        # however, the grep can be used instead of first_index
        # my ($idx) = grep { $bins->[$_] > $pos } 0 .. $#$bins;
        # $idx = -1 unless defined $idx;

        # if it did not fit in the bins, then the position must
        # be greater than the upper limit of the last bin, put
        # it in "> than upper limit of last bin".
        my $key = ($idx == -1 ? "> $bins->[-1]" : $labels->[$idx]);
        $hist{ $key } += 1;
    }

    # we're done matching, return the histogram for this sequence
    return \%hist;
}

sub read_sequence {
    my ($fh, $stop) = @_;

    my ($line, $seq);

    while ($line = <$fh>) {
        $line =~ s/\s+\z//;
        last if $line =~ $stop;
        $seq .= $line;
    }

    return ($line, \$seq);
}

__DATA__

-- Either paste your data here, or pass the name
-- of your input file on the command line

输出：

开始序列：'NR_037701'：
开始序列：'NR_002714'：
开始序列：'NR_003569'：
开始序列：'NR_002817'：
$VAR1 = {
          '700 - 800' => 0,
          '> 1000' => 10,
          '200 - 300' => 1,
          '900 - 1000' => 1,
          '800 - 900' => 1,
          '500 - 600' => 0,
          '0 - 100' => 0,
          '100 - 200' => 1,
          '300 - 400' => 0,
          '400 - 500' => 0,
          '600 - 700' => 0
        };

此外，您应该接受Chris Charley 的建议并使用Bio::SeqIO来读取序列，而不是我的自制read_sequence函数。BioPerl我只是为了回答这个问题而懒得安装。

Answer 2

这与您之前的问题几乎相同，只是间隔与序列的长度无关，因此可以只定义一次，而不是为每个序列更改它们。

该程序是对我以前的解决方案的修改。正如我所描述的，它从一组固定的值开始，@offsets从100到1000，步长为100，最终范围以20 亿或 20 亿> 1000终止。2E9这接近最大的 32 位正整数，用于捕获 1000 以上的所有偏移量。我假设您不会处理比这更大的序列？

@totals和数组初始化为零，其@counts元素数量与@offsets数组相同。

否则功能和以前一样。

use strict;
use warnings;

use List::MoreUtils 'firstval';

my $gpat = '[G]{3,5}';
my $npat = '[A-Z]{1,25}';
my $pattern = $gpat.$npat.$gpat.$npat.$gpat.$npat.$gpat;
my $regex = qr/$pattern/i;

open my $fh, '<', 'small.fa' or die $!;

my @offsets = map $_*100, 1 .. 10;
push @offsets, 2E9;
my @totals = (0) x @offsets;

my ($id, $seq);

while (<$fh>) {

  chomp;

  if (/^>(\w+)/) {
    process_seq($seq) if $id;
    $id = $1;
    $seq = '';
    print "$id\n";
  }
  elsif ($id) {
    $seq .= $_;
    process_seq($seq) if eof;
  }
}

print "Total: @totals\n";



sub process_seq {

  my $sequence = shift;

  my @counts = (0) x @offsets;

  while ($sequence =~ /$regex/g) {
    my $place = $-[0];
    my $i = firstval { $place < $offsets[$_] } keys @offsets;
    $counts[$i]++;
  }

  print "Counts: @counts\n\n";
  $totals[$_] += $counts[$_] for keys @totals;
}

输出

对你的新数据文件运行这个程序small.fa会产生

Total: 1 1 0 0 0 0 0 1 0 1 10

但是使用上一个问题的数据，sample.fa更有趣

Total: 5 4 1 0 0 2 2 1 0 0 1

Answer 3

通常，在 Perl 中，您可以通过以下方式计算模式的出现次数：

 $_ = $input;
 my $c = 0;
 $c++ while s/pattern//s;

Answer 4

我终于能够弄清楚我的代码哪里出了问题。原来是循环问题。以下代码完美运行。我已经在评论中标记了我进行修改的地方。

#!/usr/bin/perl -w

use strict;
use warnings;

my $gpat    = '[G]{3,5}';
my $npat    = '[A-Z]{1,25}';
my $pattern = $gpat . $npat . $gpat . $npat . $gpat . $npat . $gpat;
my $regex   = qr/$pattern/i;

open OUT, ">Quadindividual.refMrna.fa" or die;
open my $fh, '<', 'refMrna.fa' or die $!;

my ( $id, $seq );    # can be written as my $id; my $seq;
my @totals = ( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, );    #intialize the @total arrays.
my @thousandcounts = (0);

while (<$fh>) {

  chomp;

  if (/^>(\w+)/) {
    process_seq($seq) if $id;
    $id  = $1;
    $seq = '';
    print "$id\n";
    print OUT "$id\n";
  }
  elsif ($id) {
    $seq .= $_;
    process_seq($seq) if eof;
  }
}

print "Totals : @totals\n";
print OUT "Totals : @totals \n";

print "Thousand Counts total : @thousandcounts\n";
print OUT "Thousand Counts total : @thousandcounts\n";

sub process_seq {

  my $sequence = shift @_;

  my $subseq = substr $sequence, 0, 1000;
  my $length = length $subseq;
  print $length, "\n";

  my @offsets = map { sprintf '%.0f', $length * $_ / 10 } 1 .. 10;
  print "Offsets of 10 divisions: @offsets\n";

  my @counts = ( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, );
  my @count = (0);

  # *MODIFICATION*
  # This if loop was intialized from my @offsets above earlier
  if ( $length eq 1000 ) {
    while ( $sequence =~ /$regex/g ) {
      my $place = $-[0];
      print $place, "\n\n";

      if ( $place <= 1000 ) {
        for my $i ( 0 .. 9 ) {
          next if $place >= $offsets[$i];
          $counts[$i]++;
          last;
        }

      }

      if ( $place > 1000 ) {

        for my $i (0) {
          $count[$i]++;
          last;
        }
      }

    }    #*MODIFICATION*
         #The following commands were also subsequently shifted to ..
         #...properly compute the total

    print "Counts : @counts\n\n";

    $totals[$_] += $counts[$_] for 0 .. 9;

    print "Count : @count\n\n";

    $thousandcounts[$_] += $count[$_] for 0;
  }

  elsif ( $length != 1000 ) {

    my $substr = join ' ', unpack '(A100)*', $sequence;

    my @offsets =
        map { sprintf '%.0f', $length * $_ / ( $length / 100 ) } 1 .. 10;
    print "Offsets of 10 divisions: @offsets\n";

    my @counts = ( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, );

    while ( $sequence =~ /$regex/g ) {
      my $place = $-[0];
      print "Place : $place", "\n\n";
      for my $i ( 0 .. 9 ) {
        next if $place >= $offsets[$i];
        $counts[$i]++;
        last;
      }
    }
    print "Counts : @counts\n\n";

    $totals[$_] += $counts[$_] for 0 .. 9;

  }

}    #subroutine ends