perl - スライディングウィンドウアプローチのための perl コードの改良

Question

私はファイル（fastaファイル）を扱っています。フォーマットは次のとおりです-

>chr1
AACCCCCCCCTCCCCCCGCTTCTGGCCACAGCACTTAAACACATCTCTGC
CAAACCCCAAAAACAAAGAACCCTAACACCAGCCTAACCAGATTTAAAAT
TTTATCTTTAGGCGGTATGCACTTTTAACAAAAAANNNNNNNNNNNNNNNNNNNNNNNN
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
GCCCATCCTACCCAGCACACACACACACCGCTGCTAACCCCATACCCCGAAC
CAACCAAACCCCAAAGACACCCCCCCAAGTTTATGTAGCTTACCTCNNNN
>chrM
GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCAT
TTGGTATTTTCGTCTGGGGGGTGTGCACGCGATAGCATTGCGAGACGCTG
GAGCCGGAGCACCCTATGTCGCAGTATCTGTCTTTGATTCCTGCCTCATT
CTATTATTTATCGCACCTACGTTCAATATTTACAGGCGAACATACCTACTA
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 ファイルが大きく、このコードが大量のメモリと時間を消費していることです。より少ないメモリを使用して高速化する方法を教えてください。

ありがとう

score 4 · Accepted Answer

いつでも use strict、use warningsプログラムの開始時に、特に助けを求めているときに。多くの単純な間違いを見つけることで、多くの時間を節約できます。

このように型グロブを使用するというアイデアはどこから得たのですか? *ARGVは常に@ARGVtrue であるため、が空かどうかをテストするのには役に立たず*first、ファイルハンドルとして使用しても機能しますが、これは非常に珍しいことです。最適なのは、このようなレキシカルファイルハンドルです。

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

ただし、パラメータとして指定されたファイルを明示的に開く必要はありません: null ファイルハンドルから読み取る場合、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

score 1 · Accepted Answer

データを 2 回出力するコードが必要なので、サブルーチンに移動しました。

#!/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

score 0 · Accepted Answer

これは、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;
        }
    }
}

perl - スライディング ウィンドウ アプローチのための perl コードの改良

3 に答える 3

Related

Reference

perl - スライディングウィンドウアプローチのための perl コードの改良