view release on metacpan or  search on metacpan

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

    return keys %{$self->{RESPONSES}};
}

=item C<calculate()>

Perform the SDRS calculation.

=cut

sub calculate {
    my $self = shift;
    return if $self->{STATE} eq 'calculated';
    $self->_make_tmp;
    $self->{STATE} = "calculating";
    #
    # Calculate F Test degrees of freedom.
    #
    my $ndoses = scalar(@{$self->{DOSES}});
    my $p = 4;	#number of parameters
    $self->{FDISTR_N} = $p - 1;
    $self->{FDISTR_M} = $ndoses - $p;
    $self->{CUTOFF} = Statistics::Distributions::fdistr($self->{FDISTR_N},
							$self->{FDISTR_M},
							$self->{SIGNIFICANCE});

    my $count = scalar(keys %{$self->{RESPONSES}});
    foreach my $assay (keys %{$self->{RESPONSES}}) {
	my @data = @{$self->{RESPONSES}->{$assay}};
	my ($max, $min) = Math::NumberCruncher::Range(\@data);
	$self->{MAX}->{$assay} = $max;
	$self->{MIN}->{$assay} = $min;
	my $value = _compute_sst(\@data);	#total sum of squares
	$self->_define_ab_boundary($assay, \@data);
	$self->{SST}->{$assay} = $value;
	$self->{VALUES}->{$assay} = \@data;
    }
    my $assayNum = int($count * (1 - $self->{TRIM}));
    my %expressedAssays;
    $count = 0;
    foreach my $g (sort {$self->{MAX}->{$b} <=> $self->{MAX}->{$a}} (keys %{$self->{MAX}})) {
	last if ($count >= $assayNum);
	$expressedAssays{$g} = 1;
	$count++;
    }

    my (%pvalues, %sortedprobes, $pid, @pids, $file, $iproc, $err);
    my $tmpdir = $self->{TMPDIR};
    for ($iproc = 0; $iproc < $self->{MAXPROC}; $iproc++) {
	if ($pid = fork) {
	    print STDERR "Process $pid forked.\n" if $self->{DEBUG};
	    push @pids, $pid;
	}
	elsif (defined $pid) {
	    sleep $iproc;
	    $file = "$tmpdir/sdrs.out.$iproc";
	    open (ECOUT, ">$file") || die "can not open EC50 output file: $!\n";
	    #
	    # The following statement causes the system to write each output line
	    # to disk immediately, rather than waiting for a buffer to fill.
	    #
	    select((select(ECOUT), $| = 1)[$[]); # Change buffering to line by line.
	    my $last_dose = $self->{LDOSE};
	    my $icount = 0;
	    for (my $dose = $self->{LDOSE};
		 $dose < $self->{HDOSE} * $self->{MULTIPLE};
		 $dose *= $self->{MULTIPLE}) {
		if ($dose - $last_dose > $self->{STEP}) {
		    $dose = $last_dose + $self->{STEP};
		}
		if (($icount++ % $self->{MAXPROC}) == $iproc) {
		    print ECOUT $self->_scan_dose_point($dose);
		}
		$last_dose = $dose;
	    }
	    close ECOUT;
	    exit 0;
	}
	else {
	    die "Can't fork: $!\n";
	}
    }
    $err = 0;
    foreach $pid (@pids) {
	waitpid($pid, 0);
	if ($? != 0) {
	    print STDERR "Process $pid return status was $?\n";
	    $err = 1;
	}
    }
    
    # Now collect all the outputs together.
    %{$self->{GSCORE}} = ();
    %{$self->{GPARAM}} = ();
    %{$self->{FSCORES}} = ();
    my $multiple = $self->{MULTIPLE};
    my $step = $self->{STEP};
    @{$self->{SCANDATA}} = ();
    
    for ($iproc = 0; $iproc < $self->{MAXPROC}; $iproc++) {
	my $infile = "$tmpdir/sdrs.out.$iproc";
	open (IN, "<$infile") || die "can not open input file $infile: $!\n";
	while (<IN>) {
	    push (@{$self->{SCANDATA}}, $_);
	    chomp;
	    my ($assay, $dose, $fscore, $a, $b, $d) = split (/\t/, $_);
	    $self->{GSCORE}->{$assay}->{$dose} = $fscore;
	    $self->{GPARAM}->{$assay}->{$dose} = "$a:$b:$d";
	    if (defined $expressedAssays{$assay}) {
		$self->{FSCORES}->{$dose}->{$assay} = $fscore;
	    }
	}
	close IN;
    }
    
    #sort probesets based on F scores for every selected doses, output P values for FDR calculation
    %{$self->{SORTED_DATA}} = ();
    %{$self->{PVAL_DATA}} = ();
    foreach my $dose (keys %{$self->{FSCORES}}) {
	my %fs = %{$self->{FSCORES}->{$dose}};
	foreach my $assay (sort {$fs{$b} <=> $fs{$a} || $a cmp $b } (keys %fs)) {
	    push (@{$self->{SORTED_DATA}->{$dose}}, $assay);