2

私は Perl とスクリプト作成に非常に慣れていませんが、研究用のコードが必要です。multiFASTA ファイルに保存されている DNA 配列の 11-mer の頻度を計算しようとしています。私が見つけたいくつかのスクリプトを一緒にマージすることによって、私はこれを書きました:

#!/usr/bin/perl

$k = 11;  @bases = ('A','C','G','T');
@words = @bases; open FILE1, ">kmers.txt" or die $!;
for $i (1..$k-1)  {
   undef @newwords;
   foreach $w (@words)
   {
       foreach $b (@bases)
       {
          push (@newwords,$w.$b);
       }
   }
   undef @words;
   @words = @newwords;  
}
foreach $w (@words) {  
   print FILE1 "$w \n"; 
} 
close FILE1;   
my $input=$ARGV[0]; 
my $output=$ARGV[1];
open(IN,"<$input") || die ("Error opening $input $!"); 
open OUT, ">$output" or die $|; my $line = <IN>;  
print OUT $line; 
while ($line = <IN>) { 
   chomp $line; 
   if ($line=~m/^>/) { 
      print OUT  "\n",$line,"\n"; 
   } else { 
      print OUT $line; 
   } 
} 
print OUT "\n";

chomp $seq; chomp $k;
#obtain all distinct kmers open FILE2, ">out.txt" or die $!;

for $line (@lines) { 
   if ($line=~m/^>/) { next; } 
}
foreach($i=1; length($line) >= $k; $i++)    {   
   $line =~ m/(^.{$k})/;  
   $w{$1}{cnt}++;
   push @{$w{$1}{pos}}, $i;  
   $line= substr($seq, 1, length($line)-1);
   foreach $line (keys %kmers)    {
      print FILE2 "$kmers\n";
   }
   close FILE2; 
   close OUT;    
}

基本的に、それはファイルを読み取り、すべてのシーケンス行を別のファイルに 1 行に配置し、すべての 11mer を書き留めて「out.txt」ファイルを作成します。このファイルには、11mer 頻度のシーケンス ヘッダーを保存してもらいます。ここが (私にとって) 難しい部分です: 各シーケンスの 11mer 頻度と共にシーケンス ヘッダーを書き込むようにスクリプトに指示するにはどうすればよいでしょうか?

4

3 に答える 3

0

user2029917、宣言されていない変数に問題がありました。これにより、スクリプトをuse strict;オンにすると実行できなくなりました。いくつかの変更を加えて、少しきれいにしました。 

#!/usr/bin/perl
use strict;
use warnings;

my $in_file = $ARGV[0];
my $out_tvir = $ARGV[1];
my $k = $ARGV[2];

my %seq_hash; # key = seq_name, value = seq;
{
    # redefine the record separator
    local $/ = ">";
    open IN, "<", $in_file or die "Can't open ${in_file}: $!";
    my $in_line = <IN>; # toss the first record
    while ( $in_line = <IN> ) {
        chomp $in_line; # remove the ">" character in the end 
        my ( $seq_name, $seq ) = split( /\n/, $in_line, 2 );
        $seq =~ tr/\t\n\r//d;    # Remove whitespace
        $seq_hash{$seq_name} = uc $seq;
    }
    close IN;
}

open OUT, ">", $out_tvir or die "Can't open ${out_tvir}: $!";
foreach my $seq_name ( sort keys %seq_hash ) {
    chomp $k;
    my %kmers;
    while (length($seq_hash{$seq_name}) >= $k) {
        $seq_hash{$seq_name}=~ m/(^.{$k})/;
        $kmers{$1}++;
        $seq_hash{$seq_name}= substr($seq_hash{$seq_name}, 1, length($seq_hash{$seq_name})-1);
    }
    my $num_kmers = keys %kmers;
    my $px;
    my $logpx;
    my $H;
    foreach my $str (keys %kmers) {
        my $px=$kmers{$str}/$num_kmers;
        $logpx=log($px);
        $H -= $px * log($px);
        if ($H <= 18) {print OUT ">$seq_name\t$H\n";}
    }
}

close OUT;

exit;

これで実行されるはずですが、このスクリプトが目的の出力を生成するかどうかはわかりません。たとえば、特定の k-mer について、それが出現するすべての FASTA エントリの H' 値を出力します (これは、FASTA エントリに関係なく、常に同じ値になります。これは、出現総数と総数で計算されるためです)。 k-merの）。現在、どの k-mer が参照されているかは表示されません。これは、最後のビットを に変更することで修正できるものですがprint OUT ">$seq_name\t$str\t$H\n";、それがあなたが求めている動作であるかどうかはわかりません. 必要な出力についてさらに詳細をお知らせいただければ、さらにお役に立てる可能性があります。

于 2013-02-18T15:32:29.767 に答える
0

コードをいじった後、私はこれを作成しました：

use strict;
use warnings;
my $in_file = $ARGV[0];
my $out_tvir = $ARGV[1];
my $k = $ARGV[2];
my %seq_hash; # key = seq_name, value = seq;
{
# redefine the record separator
local $/ = ">";
open IN, "<$in_file";
my $in_line = <IN>; # toss the first record
while ( $in_line = <IN> ) {
    chomp $in_line; # remove the ">" character in the end 
    my ( $seq_name, $seq ) = split( /\n/, $in_line, 2 );
    $seq =~ tr/ \t\n\r//d;    # Remove whitespace
    $seq_hash{$seq_name} = uc $seq;
}
close IN;
}

open OUT, ">$out_file";
open OUT2, ">$out_tvir";
foreach my $seq_name ( sort keys %seq_hash ) {
chomp $k;
%kmers = ();
while (length($seq_hash{$seq_name}) >= $k)
    {
    $seq_hash{$seq_name}=~ m/(^.{$k})/;
    $kmers{$1}++;
    $seq_hash{$seq_name}= substr($seq_hash{$seq_name}, 1,         length($seq_hash{$seq_name})-1);
    }
    $num_kmers = keys %kmers;
$px=();
$logpx=();
my $H=();
foreach $str (keys %kmers)
{
    my $px=$kmers{$str}/$num_kmers;
    $logpx=log($px);
    $H -= $px * log($px);
    if ($H <= 18) {print OUT2 ">$seq_name\t$H\n";}
}
}
close OUT;

...最後の「if ($H...」部分を省略し、各シーケンスに関連付けられたすべての H 値をリストすることでスクリプトにジョブを実行させると、どのような作業が行われますか?私には理由がわかりません。けれど。

于 2013-02-01T16:47:19.287 に答える