-1

配列を 200 ずつ増やしてループしようとしています。その 200 内にサブ配列を作成し、サブ配列で計算を実行してから、読み取りフレームを次の 200 のセットに移動して繰り返します。ただし、2 番目の計算を行うたびに、新しいサブ配列には最初のサブ配列の値が含まれます。ループundefの次の繰り返しで再利用するためにサブアレイをリセットまたは再利用できないようです。for

完全なコードは github で入手できます: https://github.com/bsima/yeast-TRX

コードを実行するには、スクリプトと同じディレクトリに次のゲノム ファイルが必要です: https://raw.github.com/bsima/yeast-TRX/master/data/bayanus/genome.csv

関連するコードは次のとおりです。

#!/usr/bin/perl

use strict;
use warnings;
use diagnostics;
use Statistics::Descriptive;

my $species = "bayanus";

open(SPECIES, "<./genome.csv") || die "Cannot open file: $!\n";
my @text = <SPECIES>;

my $geneNameRe  = qr/(eY\w{5}[CW])/;
my $geneRe      = qr/,([atgcATGC-]+)/;

foreach my $line (@text) {
    chomp($line);
    if ( defined $line && $line =~ m/e.{4,},[acgtACGT-]+/ ) {

        my $gene     = match($geneRe,$line);
        my $geneName = match($geneNameRe,$line);

        # Smoothing function
        #
        # This moves through the gene data and counts the position until it arrives
        # at the end of the smoothing window (e.g. 200). Then it calculates the average
        # of the selected data set and outputs it into the respecive `smooth.csv` file 
        # in the following format, to be read later by R's graphing functions:
        #       gene,position,trx.mean,energy.mean
        open(my $smooth, ">>./smooth.csv") || die "Cannot open file $!";
        my $smoothingWindow = 200;
        # Loop through every 200 characters
        for ( my $smoothing = 0; $smoothing < length($gene); $smoothing=$smoothing+$smoothingWindow ) {

            my @trxValues;
            my @energyScores;

            # Loop through every character
            for ( my $position = 0; $position <= $smoothingWindow; $position++ ) {
                my $dinucleotide = substr($gene,$position,2); 

                if ( $dinucleotide =~ m/[actgACTG]{2}/ ) {     
                    my $trxValue = trxScore($dinucleotide);
                    my $energyScore = energyScore($dinucleotide);

                    #print "trxValue = $trxValue\n"; 
                    push @trxValues, $trxValue;
                    push @energyScores, $energyScore;
                }                    
            }

            # Now run the calculations and print to SMOOTH
            my $trxStat = Statistics::Descriptive::Full->new();
            $trxStat->add_data(@trxValues);
            my $trxMean = $trxStat->mean();
            my $energyStat = Statistics::Descriptive::Full->new();
            $energyStat->add_data(@energyScores);
            my $energyMean = $energyStat->mean();

            print "$geneName $smoothing: TRX Values @trxValues\n"; 
            print "$geneName $smoothing: TRX Mean is $trxMean.\n";

            print $smooth $geneName . "," . $smoothing . "," . $trxMean . "," . $energyMean . "\n";

        }
        close $smooth;
    }
}

# This just makes it easy to do regex matches
sub match {
        my ( $re, $text ) = @_;
        if ( $text =~ $re ) {
                return $1;
        }
}

# @name trxScore
# @description Calculates and returns the TRX value of a given phosphate linkage
# @param $dinucleotide {string} The nucleotide to be checked
# @return {integer} The TRX value
sub trxScore {

    my ( $dinucleotide ) = @_;

        my %trxScores = (
            qr/(CG)/ => 43,
            qr/(CA)/ => 42,
            qr/(TG)/ => 42,
            qr/(GG)/ => 42,
            qr/(CC)/ => 42,
            qr/(GC)/ => 25,
            qr/(GA)/ => 22,
            qr/(TC)/ => 22,
            qr/(TA)/ => 14,
            qr/(AG)/ =>  9,
            qr/(CT)/ =>  9,
            qr/(AA)/ =>  5,
            qr/(TT)/ =>  5,
            qr/(AC)/ =>  4,
            qr/(GT)/ =>  4,
            qr/(AT)/ =>  0,
            qr/(.-)/ =>  0,
            qr/(-.)/ =>  0   # These last two are necessary for dealing with missing values in the genomes
        );

    foreach my $re (keys %trxScores) {
        if ( match($re,$dinucleotide) ) {
            return $trxScores{$re};
        } 
    }
    return 0;
}

# @name energyScore
# @description Calculates and returns the delta-E value of a given phosphate linkage
# @param $dinucleotide {string} The nucleotide to be checked
# @return {integer} The delta-E value
sub energyScore {

    my ( $dinucleotide ) = @_;

        my %energyScores = (
            qr/(AA)/ => -18.5,
            qr/(AC)/ => -19.0,
            qr/(AG)/ => -23.6,
            qr/(AU)/ => -15.7,
            qr/(CA)/ => -20.0,
            qr/(CC)/ => -21.4,
            qr/(CG)/ => -26.9,
            qr/(CU)/ => -17.2,
            qr/(GA)/ => -23.7,
            qr/(GC)/ => -22.9,
            qr/(GG)/ => -24.3,
            qr/(GU)/ => -18.9,
            qr/(UA)/ => -19.6,
            qr/(UC)/ => -28.2,
            qr/(UG)/ => -23.3,
            qr/(UU)/ => -15.8,
            qr/(.-)/ =>     0,
            qr/(-.)/ =>     0   # These last two are necessary for dealing with missing values in the genomes
        );

    foreach my $re (keys %energyScores) {
        if ( match($re,$dinucleotide) ) {
            return $energyScores{$re};
        } 
    }
    return 0;
}

そのコード ブロックが含まれるより大きなスクリプトを実行すると、その内容が@trxValues端末に出力され、変数が適切にインクリメントされていることがわかるので、ループ$smoothingのステップ部分に問題はありません。forしかし、内容は@trxValues決して変わりません。そのため、最初のループの終わりと 2 番目のループの始まりの間のどこかでfor、配列@trxValuesが適切にクリアされません。に対しても同じことが起こり@energyScoresます。

これに関するアイデアはありますか？undefループの最後で使用して、forループの最後でゼロに設定しようとしました@trxValuesが、どちらも機能しませんでした。私は完全に困惑しています。

編集：コードブロックを単独で実行できるように変更しました。

4

1 に答える 1

0

内側の for ループで常に同じ最初の平滑化ウィンドウを使用しているようです。

変更してみてください:

for ( my $position = 0; $position <= $smoothingWindow; $position++ ) {

に

for ( my $position = $smoothing; $position < $smoothing + $smoothingWindow; $position++ ) {

内側のループはスムージング ウィンドウとlength($gene)(常に 200 の偶数倍であることがわかっている場合を除きますか?) の両方に対してチェックする必要がありますが、

length($gene)また、一度に1 つの位置に移動しているのに、substr が 2 つの文字を抽出していることも奇妙に思えます。

于 2013-10-27T18:53:50.313 に答える