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個のサブディビジョンのカウントで同じ問題に直面しています。

この計算に正しいコードを提供できれば幸いです。

私が書いたコードは以下の通りです。これは、ここでの私の以前の質問の1つに対するボロディンの答えから大きく派生しています：Perl：配列要素全体でパターンを検索する

彼の答えはここにあります：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

これは前の問題とほとんど同じですが、間隔がシーケンスの長さに依存しないため、シーケンスごとに変更するのではなく、一度だけ定義できる点が異なります。

このプログラムは、私の以前のソリューションを変更したものです。@offsets私が説明したように、それはからの固定された値のセットで始まり、100ステップ1000で始まり100、最終的な範囲はまたは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 2

以下はうまくいくようです。遊んでいる間、私はあなたが投稿したデータを__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 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