view release on metacpan or  search on metacpan

lib/Bio/SDRS.pm  view on Meta::CPAN

	$low = shift @$range;
	$high = $low;
    }
    return "$high-$low";
}

sub _scan_dose_point {
    my $self = shift;
    my $dose = shift;
    my $ecstring = '';
    my $m = $self->{FDISTR_M};
    my $n = $self->{FDISTR_N};
    foreach my $assay (keys %{$self->{VALUES}}) {
	my ($al, $ah, $astep) = split (/\:/, $self->{ARANGE}->{$assay});
	my ($bl, $bh, $bstep) = split (/\:/, $self->{BRANGE}->{$assay});
	my ($a, $b, $d, $best_sse) = _compute_best_sse($al, $ah, $astep,
						       $bl, $bh, $bstep,
						       $dose,
						       $self->{VALUES}->{$assay},
						       $self->{DOSES});
	if ($best_sse == 0 or $n == 0) {
	    carp sprintf("best_sse($best_sse) = 0 or n($n) = 0 for $assay. arange = %s  brange = %s\n",
			 $self->{ARANGE}->{$assay},
			 $self->{BRANGE}->{$assay});
	    $ecstring .= "$assay\t" .
		sprintf("%.5f", $dose) .
		    "\t-1\t-1\t-1\t-1\n";
	}
	else {
	    my $f_score = ($self->{SST}->{$assay} - $best_sse) * $m / ($best_sse * $n);
	    $f_score = sprintf("%.3f", $f_score);
	    $a = sprintf("%.3f", $a);
	    $b = sprintf("%.3f", $b);
	    $d = sprintf("%.3f", $d);
	    $ecstring .= "$assay\t" .
		sprintf("%.5f", $dose) .
		    "\t$f_score\t$a\t$b\t$d\n";
	}
    }
    return $ecstring;
}

sub _define_ab_boundary {
    my $self = shift;
    my ($assay, $data) = @_;
    my @sorted = sort {$a<=>$b} @$data;
    my @sorted_copy = @sorted;
    my @brange = splice (@sorted, $#sorted-5, 6);
    $self->_find_step($assay, 'b', \@brange);
    my @arange = splice (@sorted_copy, 0, 6);
    $self->_find_step($assay, 'a', \@arange);
}

sub _find_step {

    my $self = shift;
    my ($assay, $pam, $data) = @_;
    my $stdev = Math::NumberCruncher::StandardDeviation($data);
    my $mean = Math::NumberCruncher::Mean($data);
    my $cutoff = $mean / 10;
    $stdev = $cutoff if ($stdev eq 'NaN' || $stdev < $cutoff);
    my ($l, $h, $step);
    $step = $stdev / 2.5;
    $step = sprintf("%.3f", $step);
    if ($step == 0.0) {
	carp "Data range too small for $assay -- step size raised to 0.001.\n";
	$step = "0.001";
    }
    if ($pam eq 'a') {
	$h = $mean + 2.3*$stdev;
	$l = $mean - 2*$stdev;
	$l = $self->{MIN}->{$assay} if ($l <= 0);
	$h = sprintf("%.3f", $h);
	$l = sprintf("%.3f", $l);
	$self->{ARANGE}->{$assay} = "$l:$h:$step";
    } elsif ($pam eq 'b') {
	$h = $mean + 2*$stdev;
	$l = $mean - 2.3*$stdev;
	$l = $self->{MIN}->{$assay} if ($l <= 0);
	$h = sprintf("%.3f", $h);
	$l = sprintf("%.3f", $l);
	$self->{BRANGE}->{$assay} = "$l:$h:$step";
    }
}

sub _compute_sst {
    my $data = shift;
    my $sst = 0;
    my $mean = Math::NumberCruncher::Mean($data);
    
    foreach my $data (@$data) {
	$sst += ($data - $mean)**2;
    }
    return $sst;
}

sub _system_with_check {

    my $command = shift;
    my $echo = shift;

    if (defined $echo and $echo) {
	&_dated_mesg($command);
    }
    my $status = system("$command");
    if ($status != 0) {
	croak "Shell command: $command returned $status error code.\n";
    }
}

sub _dated_mesg {

    # Print a dated message to STDERR
    my $mesg = shift;

    my $date = &_date;
    if (substr($mesg, -1, 1) ne "\n") {
	$mesg .= "\n";
    }
    print STDERR "At $date: $mesg";
}