Уточнение Perl-кода для подхода скользящего окна - PullRequest
Новая
       10

Уточнение Perl-кода для подхода скользящего окна

0 голосов
/ 29 марта 2012

Я работаю с файлом (файл fasta), вот формат - 

>chr1
AACCCCCCCCTCCCCCCGCTTCTGGCCACAGCACTTAAACACATCTCTGC
CAAACCCCAAAAACAAAGAACCCTAACACCAGCCTAACCAGATTTCAAAT
TTTATCTTTAGGCGGTATGCACTTTTAACAAAAAANNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
GCCCATCCTACCCAGCACACACACACCGCTGCTAACCCCATACCCCGAAC
CAACCAAACCCCAAAGACACCCCCCACAGTTTATGTAGCTTACCTCNNNN
>chrM
GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCAT
TTGGTATTTTCGTCTGGGGGGTGTGCACGCGATAGCATTGCGAGACGCTG
GAGCCGGAGCACCCTATGTCGCAGTATCTGTCTTTGATTCCTGCCTCATT
CTATTATTTATCGCACCTACGTTCAATATTACAGGCGAACATACCTACTA
AAGTGTGTTAATTAATTAATGCTTGTAGGACATAATAATAACAATTGAAT
GTCTGCACAGCCGCTTTCCACACAGACATCATAACAAAANAATTTCCACC

Я хочу использовать метод скользящих окон (не перекрывающиеся окна, размер = 50). Я хочу найти координаты для каждого символа в окне 50bp, но не включая N. Выход для первого chr1 должен быть -> 

chr1   0       50
chr1   50      100
chr1   100     215 
chr1   215     265

Код - 

use warnings; 
*ARGV or die "No input file specified";
open *first, '<',$ARGV[0] or die "Unable to open input file: $!";
$start=1;
while(<first>) {
    chomp;
    if ( /(>)(\w)/ ) {   #taking lines which have name of chromosome
    @arr=split(">");  #splitting at ">" character and in $arr[1], there is chr name now


        if (defined @array){

            foreach (@array){
            $length++;      

                if($_ ne N){
                    $non++;
                    if ($non == 50){

                    print $chr,"\t",$start,"\t",$length,"\n";
                    $start=$length;
                    $non=0;

                    }
                }       
            }
        }

        undef @array;  
        $length=0;
        $non=0;
        $start=0;
    }

    else {

        @count=split(//, $_); #splitting each character in line

        push(@array,@count);  #storing each character in array till we find next chromosome

        $chr=$arr[1];
    } 

}




foreach (@array){
        $length++;

          if($_ ne N){
          $non++;
               if ($non == 50){

        print $chr,"\t",$start,"\t",$length,"\n";
        $start=$length;
        $non=0;

               }
          }

}

Дело в том, что мой файл fasta большой, и этот код занимает много памяти и времени. Можете ли вы дать предложения, как я могу сделать это быстрее, используя меньше памяти.

Спасибо

Ответы [ 3 ]

4 голосов
/ 29 марта 2012

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

Откуда у вас появилась идея использовать typeglobs таким образом?*ARGV is всегда true, поэтому бесполезно проверять, является ли @ARGV пустым, и использование *first в качестве дескриптора файла будет работать, но это очень необычно.Лучше всего использовать лексический дескриптор файла, например, 

open my $first, '<', $ARGV[0] or die $!;

Однако нет необходимости явно открывать файл, указанный в качестве параметра: Perl сделает это неявно для вас, если вы прочитаете из нулевого дескриптора файла <>.

Эта программа, кажется, делает то, что вам нужно. 

use strict;
use warnings; 

use constant WINDOW => 50;

@ARGV or die "No input file specified";

my ($key, $pos, $start, $size);

while (<>) {

  if ( /^>(.+?)\s/ ) {
    $key = $1;
    $pos = $size = 0;
    undef $start;
    next;
  }

  chomp;

  for (split //) {
    next unless /[ATGC]/;
    $start //= $pos;
    $size++;
    if ($key and $size == WINDOW) {
      printf "%-6s %4d %4d\n", $key, $start, $pos + 1;
      undef $start;
      $size = 0;
    }
  }
  continue {
    $pos++;
  }
}

вывод 

chr1      0   50
chr1     50  100
chr1    100  215
chr1    215  265
chrM      0   50
chrM     50  100
chrM    100  150
chrM    150  200
chrM    200  250
1 голос
/ 29 марта 2012

Поскольку вам нужен код для вывода данных дважды, я переместил его в подпрограмму. 

#!/usr/bin/perl

use strict ;
use warnings ;

if( ! @ARGV  ) {
  die "No input file specified";
}

open my $file , '<', $ARGV[0] or die "Unable to open input file: $!";
my ( $chromosome , $start ) = ( undef , 1 ) ;

my @array = () ;
while(<$file>) {
  chomp;

  if ( m/^>(\w+)/ ) { # New chromosome
    my $new_chromosome = $1 ; # Save the new chromosome name temporarily
    if( @array ) {
      split_sequence( $chromosome , \@array ) ;
    }
    @array = () ;
    $chromosome = $new_chromosome ;
  } else {
    push @array , split( // ) ;
  }
}

split_sequence( $chromosome , \@array ) if @array ;

sub split_sequence {
  my ( $chromosome , $arrayref ) = @_ ;

  printf "%-10.10s  %d (total length)\n" , $chromosome , $#{ $arrayref } ;
  my ( $start , $nonN ) = ( 0 , 0 ) ;
  for( my $i = 0 ; $i <= $#{ $arrayref } ; $i++ ) {
    if( $arrayref->[$i] ne 'N' ) {
      $nonN++ ;
      if( $nonN == 50 ) {
        printf "%-10.10s  %8d  %8d\n" , $chromosome , $start , $i ;
        $start = $i + 1 ;
        $nonN = 0 ;
      }
    }
  }
  if( $#{ $arrayref } > $start ) { # Incomplete window leftover ...
                                   # less than 50 bases long
    printf "%-10.10s  %8d  %8d **\n" , $chromosome , $start , $#{ $arrayref } ;
  }
}

Вывод: 

perl SO002.pl SO002.fasta
chr1             299 (total length)
chr1               0        49
chr1              50        99
chr1             100       214
chr1             215       264
chr1             265       299 **
chrM             300 (total length)
chrM               0        49
chrM              50        99
chrM             100       149
chrM             150       199
chrM             200       249
chrM             250       300
0 голосов
/ 31 марта 2012

Вот решение, которое использует модуль Bio :: SeqIO для анализа файла fasta. 

#!/usr/bin/perl
use strict;
use warnings;
use Bio::SeqIO;

use constant WINDOW => 50;

my $in  = Bio::SeqIO->new(-file   => "fasta.txt" ,
                          -format => 'Fasta');

while ( my $seq = $in->next_seq() ) {
    my $count = 0;
    my $beg_pos = 0;
    local $_ = $seq->seq;
    while (/(.)/g) {
        ++$count if $1 =~ /[TAGC]/;

        if ($count == WINDOW) {
            $count = 0;
            printf "%s %d %d\n", $seq->id, $beg_pos, pos() - 1;
            $beg_pos = pos();
        }
        elsif (pos == length) { # have read last char in string
            printf "%s %d %d\n", $seq->id, $beg_pos, pos() - 1;
        }
    }
}
...