2

世界中の Perl マスターにこんにちは。

私はプログラミングに関して別の問題を抱えています。特定の入力番号を持つプロテオム fasta ファイルからランダムなシーケンスを選択するプログラムをコーディングしています。

一般的な fasta ファイルは次のようになります。

>seq_ID_1 説明など

>seq_ID_2 説明など ASDGDSAHSAHASDFRHGSDHSDGEWTSHSDHDSHFSDGSGASGADGHHAH ASDSADGDASHDASHSAREWAWGDASHASGASGASG

等々.......

文字はアミノ酸ペプチドを表します。

したがって、1000 シーケンスの fasta ファイルがあり、それらの 63.21% (632.1 シーケンス) を取得したいと考えています。ただし、シーケンスは浮動小数点数にすることはできないため、0.5 を超える場合は切り上げ、0.5 未満の場合は切り捨てます。

これは、ランダム シーケンス サブセットを生成するための私のコードですが、少しうまく動作しません。

#!/usr/bin/perl

#Selecting 63.21% of random sequnces from a proteom file.
use strict;
use warnings;
use List::Util qw(shuffle);

#Give the first argument as a proteom file.
if (@ARGV != 1)
{
    print "Invalid arguments\n";
    print "Usage: perl randseq.pl [proteom_file]";
    exit(0);
}

my $FILE = $ARGV[0];
my $i = 0;
my %protseq = {};
my $nIdx = 0;

#Extraction and counting of the all headers from a proteom file.
open(LIST,$FILE);
open(TEMP1, ">temp1");
while (my $line = <LIST>){
    chomp $line;
    if ($line =~ />(\S+) (.+)/){
        $i++;
        print TEMP1 $1,"\n";
    }
}
close(LIST);
close(TEMP1);

#Selection of random headers for generating a random subset of the proteom file.
my $GET_LINES = RoundToInt ($i*0.6321);

my @line_starts;
open(my $FH,'<','temp1');
open(TEMP2, ">temp2");

do {
     push @line_starts, tell $FH
} while ( <$FH> );

my $count = @line_starts;

my @shuffled_starts = (shuffle @line_starts)[1..$GET_LINES+1];

for my $start ( @shuffled_starts ) {

     seek $FH, $start, 0
         or die "Unable to seek to line - $!\n";

     print TEMP2 scalar <$FH>;
}
close(TEMP2);

#Assigning the sequence data to randomly generated header file.
open(DATA,'<','temp2');
while(my $line = <DATA>)
{
    chomp($line);
    $line =~ s/[\t\s]//g;
    if($line =~ /^([^\s]+)/)
    {
        $protseq{$1}++;
    }
}
close(DATA);

open(DATA, "$FILE");
open(OUT, ">random_seqs.fasta");
while(my $line = <DATA>)
{
    chomp($line);
    if($line =~ /^>([^\s]+)/)
    {
        if($protseq{$1} ne "")
        {

            $nIdx = 1;
            print OUT "$line\n";
        }
        else
        {
            $nIdx = 0;
        }
    }
    else
    {
        if($nIdx == 1)
        {
            print OUT "$line\n";
        }
    }
}
close(DATA);
close(OUT);

#subroutine for rounding
sub RoundToInt {
  int($_[0] + .5 * ($_[0] <=> 0));
}

system("erase temp1");
system("erase temp2");
exit;

ただし、適切な数のシーケンスが得られることもあれば、もう 1 つのシーケンスが得られることもあります。どうすればそれを取り除くことができますか...何かアイデアをお願いします?

またはおそらくより短いコードですか？

ここでは、75 酵母プロテオム ファイルを取得できます。[http://www.peroxisomedb.org/Download/Saccharomyces_cerevisiae.fas][1]

これをすぐに修正できることを願っています... :(

4

2 に答える 2

4

あなたのアプローチは問題ないように見えますが、不必要に複雑です。私は次のようにします：

use strict;
use warnings;

# usage: randseq.pl [fraction] < input.fasta > output.fasta
my $fraction = (@ARGV ? shift : 0.6321);

# Collect input lines into an array of sequences:
my @sequences;
while (<>) {
    # A leading > starts a new sequence. (The "\" is only there to
    # avoid confusing the Stack Overflow syntax highlighting.)
    push @sequences, [] if /^\>/;
    push @{ $sequences[-1] }, $_;
}

# Calculate how many sequences we want:
my $n = @sequences;
my $k = int( $n * $fraction + 0.5 );
warn "Selecting $k out of $n sequences (", 100 * $k / $n, "%).\n";

# Do a partial Fisher-Yates shuffle to select $k random sequences out of $n:
foreach my $i (0 .. $k-1) {
    my $j = $i + int rand($n-$i);
    @sequences[$i,$j] = @sequences[$j,$i];
}

# Print the output:
print @$_ for @sequences[0 .. $k-1];

このコードは、入力ファイルの内容全体をメモリに読み込むことに注意してください。入力ファイルが大きすぎて、その小さなサブセットのみが必要な場合は、リザーバー サンプリングを使用して、それ以上保持せずに任意の大きなコレクションからk 個のランダム シーケンスを選択することができます。

use strict;
use warnings;

my $k = (@ARGV ? shift : 632);  # sample size: need to know this in advance

# Use reservoir sampling to select $k random sequences:
my @samples;
my $n = 0;  # total number of sequences read
my $i;      # index of current sequence
while (<>) {
    if (/^\>/) {
        # Select a random sequence from 0 to $n-1 to replace:
        $i = int rand ++$n;
        # Save all samples until we've accumulated $k of them:
        $samples[$n-1] = $samples[$i] if $n <= $k;
        # Only actually store the new sequence if it's one of the $k first ones:
        $samples[$i] = [] if $i < $k;
    }
    push @{ $samples[$i] }, $_ if $i < $k;
}

warn "Only read $n < $k sequences, selected all.\n" if $n < $k;
warn "Selected $k out of $n sequences (", 100 * $k / $n, "%).\n" if $n >= $k;

# Print sampled sequences:
print @$_ for @samples;

ただし、入力シーケンスの特定の部分が本当に必要な場合は、最初にファイルの別のパスでそれらをカウントする必要があります。

上記の両方のプログラムは、副作用として、サンプリングされたシーケンスを均一にシャッフルします。(実際、リザーバー サンプリング アルゴリズムを意図的に微調整して、nとkのすべての値に対してシャッフルが均一になるようにしました。) それが望ましくない場合は、印刷する前に、好きな基準に従っていつでもシーケンスを並べ替えることができます。

于 2013-03-21T13:39:06.777 に答える
1

一時ファイルの代わりにラウンド数と配列に spritf 関数を使用しました

#!/usr/bin/perl

use strict;

if (@ARGV != 1)
{
    print "Invalid arguments\n";
    print "Usage: perl randseq.pl [proteom_file]";
    exit(0);
}

my $FILE = $ARGV[0];

open(LIST,"<$FILE");

my @peptides;
my $element;
while (my $line = <LIST>){
      if ($line =~ />.*/) {
      push (@peptides, $element);
      $element=$line;
      }
      else {
      $element.=$line;
      }
}
close(LIST);

my $GET_LINES = sprintf("%.0f",$#peptides*0.6321);

my @out;
for (0..$GET_LINES) {
    my $index=$#peptides;
    push (@out, $peptides[int(rand($index))]);
    splice(@peptides, $index, 1);
}

open (OUT, '>out.fasta');
foreach (@out) {
  print OUT $_."\n";
}

exit;
于 2013-03-21T13:56:57.097 に答える