-1

我试图循环一个数组，递增 200。我在 200 内构建一个子数组，对子数组运行计算，然后将阅读框移动到下一组 200 并重复。但是，每当我进行第二次计算时，新的子数组都会包含第一个子数组的值。就像我无法重置或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的 step 部分没有任何问题，for但内容@trxValues永远不会改变。因此，在第一个循环结束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++ ) {

Though the inner loop should be checking against both the smoothing window and length($gene) (unless you know it is always an even multiple of 200?)

It also seems odd that you are going up to length($gene) one position at a time but your substr is extracting two characters.

于 2013-10-27T18:53:50.313 回答