Как я могу сделать гистограмму появления определенных шаблонов в файле FASTA?

Я написал Perl-скрипт для следующего вопроса по биоинформатике, но, к сожалению, возникла проблема с выводом.

Вопрос

1) Из файла, состоящего из 40 000 уникальных последовательностей, уникальных, означающих номера идентификаторов последовательностей, извлеките следующий шаблон.

 $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 встречается 5 раз при длине >1000. Для последовательности 2, скажем, $pattern встречается 3 раза при длине> 1000. Тогда общее количество должно быть 5+3=8.

Вместо этого мой результат выглядит следующим образом: (5 + 4 + 3 + 2 + 1) + (3 + 2 + 1) = 21, то есть совокупный итог.

Я столкнулся с той же проблемой с подсчетом первых 10 подразделений по 100 символов в каждом.

Я был бы признателен, если бы правильный код мог быть предоставлен для этого расчета.

Код, который я написал, приведен ниже. Он в значительной степени основан на ответе Бородина на один из моих предыдущих вопросов здесь: 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

Спасибо, что нашли время, чтобы пройти через мою проблему.

Любая помощь и вклад будут глубоко оценены.

Спасибо, что нашли время разобраться с моей проблемой!


arrays regex bioinformatics perl fasta
person Neal    schedule 29.06.2012    source источник
comment
Похоже, вы упускаете смысл гиперссылок и форматирования кода, ну да ладно...   -  person Sinan Ünür    schedule 29.06.2012
comment
Здравствуй, Синан. Я только начал писать коды, я пытался сделать их читабельными. Прошу прощения, если код не выглядит хорошо отформатированным. Будучи биологом, мир программирования для меня совершенно нов, но я пытаюсь...   -  person Neal    schedule 30.06.2012
comment
Я имел в виду тот факт, что вы отменили мои правки.   -  person Sinan Ünür    schedule 30.06.2012
comment
О... должно быть, это произошло непреднамеренно, потому что, пока я редактировал, я получил уведомление о другом редактировании, но я не знал, что делать, будучи здесь новичком.. Я должен был проверить новые правки, прежде чем приступить к моему. собственный. Я буду помнить это впредь. Я также вижу, что название вопроса было полностью изменено.   -  person Neal    schedule 30.06.2012

Ответы (4)


arrow_upward
2
arrow_downward

Это почти то же самое, что и ваша предыдущая проблема, за исключением того, что интервалы не зависят от длины последовательности и поэтому могут быть определены просто один раз вместо того, чтобы менять их для каждой последовательности.

Эта программа является модификацией моего предыдущего решения. Как я описал, он начинается с фиксированного набора значений в @offsets от 100 до 1000 с шагом 100, а конечный диапазон > 1000 заканчивается на 2E9 или 2 миллиардах. Это близко к максимальному положительному 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
person Borodin    schedule 01.07.2012

arrow_upward
2
arrow_downward

Следующее, кажется, работает. Играя, я поместил данные, которые вы разместили, в разделе __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

Выход:

Start sequence: 'NR_037701':
Start sequence: 'NR_002714':
Start sequence: 'NR_003569':
Start sequence: '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
        };

Кроме того, вам следует воспользоваться советом Криса Чарли и использовать Bio::SeqIO для чтения последовательностей вместо моей самодельной функции read_sequence. Мне просто было лень устанавливать BioPerl только для того, чтобы ответить на этот вопрос.

person Sinan Ünür    schedule 29.06.2012
comment
Здравствуй, Синан! Спасибо за ваш ответ. Я новичок в Perl с менее чем 2-месячным опытом программирования, и ну... код для меня очень продвинутый. На самом деле, я даже не могу запустить его :| а тем более понять... :( - person Neal; 30.06.2012
comment
У меня установлен модуль BioPerl. Я установил его через Perl Package Manager. Но почему-то он не установился должным образом и продолжает выдавать сообщения об ошибках, в которых говорится, что он не может найти Bio/SeqIO.pm. - person Neal; 30.06.2012
comment
Конечно, невозможно понять, почему вы не можете запустить его, если не дадите никакой информации. Во всяком случае, я добавил некоторые комментарии и перешел на использование Data::Dumper вместо YAML и print вместо say. - person Sinan Ünür; 30.06.2012
comment
Здравствуй еще раз, Синан! За последние пару часов произошло замечательное событие. Наконец-то я смог увидеть, где я ошибался в своем собственном коде. Оказалось, проблема в цикле. Я также разместил свою исправленную версию. - person Neal; 30.06.2012
comment
Вы также представляете очень интересный взгляд на решение этой проблемы с помощью гистограммы. Я не мог подумать об этом. Так что спасибо вам за то, что помогли мне лучше узнать о замечательном языке Perl. Мне нужно многому научиться! - person Neal; 30.06.2012
comment
Что касается проблемы с модулем Bioperl, я получаю следующее сообщение об ошибке: Can't locate Bio/SeqIO.pm in @INC <@INC contains C:/Strawberry/perl/site/lib C:/Strawberry/perl/vendor/lib C:/Strawberry/perl/lib .> at C:\perltests\Transcripts\26thjune\quadatwhatposn2.pl line 4. BEGIN failed--compilation aborted at C:\perltests\Transcripts\26thjune\quadatwhatposn2.pl line 4. - person Neal; 30.06.2012
comment
@Neal При всем уважении, вы не имеете смысла: подсчет вещей в корзинах - это определение гистограммы. Это то, над чем вы работали все это время. Кроме того, я не уверен, что вы делаете это нарочно, но полностью игнорировать мой ответ, не предоставлять по нему никаких отзывов и даже не объяснять, почему он не работает, очень и очень грубо. В конце концов, я приложил некоторые усилия, чтобы понять ваш вопрос, а не публиковать одноразовый ответ, как единственный другой, который вы получили в ответ на свой запрос. - person Sinan Ünür; 01.07.2012
comment
Привет! Синан, мне очень жаль, если это так выглядит. Как я уже упоминал, я новичок в Perl, у меня вообще нет опыта программирования. Я прохожу стажировку по биоинформатике в течение 2 месяцев, и мне посоветовали эту книгу «Начальный Perl для биоинформатики». В объем книги не входило использование встроенных модулей. Возможно, это не очень продвинутая вещь, чтобы использовать модули, но я не смог найти время, чтобы прочитать об этом, мне нужно пару дней. Такие команды, как qw, qq, qr, являются новыми, и я столкнулся с ними только в последние пару дней. Как и функции map, sprintf. - person Neal; 01.07.2012
comment
А функции определения и повтора совершенно новые. Я, конечно, очень высоко оценил ваш ответ, я проголосовал за него, потому что он помог мне узнать о многих новых функциях Perl. А про гистограмму я знал только в статистике, а не в компьютерных языках. И только теперь, когда вы упомянули об этом, я вижу это таким образом. Я также не знаю, что имеется в виду под мусорными баками. На самом деле, я вовсе не игнорирую или что-то в этом роде, но, пожалуйста, дайте мне немного времени, чтобы ответ усвоился. Это действительно ново для меня. У меня есть только эта книга и онлайн-документация Perl, а учителя нет. Я стараюсь.. - person Neal; 01.07.2012

arrow_upward
0
arrow_downward

Как правило, в Perl вы можете подсчитать появление шаблона следующим образом:

 $_ = $input;
 my $c = 0;
 $c++ while s/pattern//s;
person Ωmega    schedule 29.06.2012
comment
Привет Омега. Спасибо за ваш ответ. Я попытался вставить приведенные выше строки кода в свою программу, но это не сработает. Не могли бы вы предложить альтернативу? Я указал точную область моего кода выше, где проблема. - person Neal; 30.06.2012
comment
Хм, я только что заметил, что весь мой код дает ошибочные результаты. Окончательная оценка суммируется для всех частей, а не только для частей длиной более 1000 символов. - person Neal; 30.06.2012

arrow_upward
0
arrow_downward

Наконец-то я смог понять, где я ошибался в своем коде. Оказалось, проблема в цикле. Следующий код работает отлично. Я отметил в комментариях места, где я сделал модификацию.

#!/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
person Neal    schedule 30.06.2012