Summary

G::Seq::Markov - Perl extension for blah blah blah

Package variables

Included modules

Inherit

Exporter

Synopsis

  use G::Seq::Markov;
  blah blah blah

Description

Stub documentation for G::Seq::Markov was created by h2xs. It looks like the
author of the extension was negligent enough to leave the stub
unedited.
Blah blah blah.

Methods

Methods description

None available.

Methods code

sub DESTROY {

    my $self = shift;

}

sub markov {

    &opt_default(length=>6, mincount=>10, filename=>"markov.csv",output=>"stdout");
    my @args = opt_get(@_);
    my $gb = opt_as_gb(shift @args);
    my $filename = opt_val("filename");
    my @aSortedTable = ( );
    my $iTotalNucs = 0;
    my @ahNucsTable = ( );
    my %oe;

    for (my $iCounter = 0; $iCounter <= 32; $iCounter++) {
	$ahNucsTable[$iCounter] = { };
    }

    my $rhTmp;
    foreach $rhTmp (@ahNucsTable) { undef %$rhTmp; }
    $iTotalNucs = 0;

    my($nucs, $char);
    $nucs = '';
    foreach $char (split(//, $gb->{SEQ})) {
	$iTotalNucs++;
	$nucs .= $char;
	if (opt_val("length") < $iTotalNucs) {
	    substr($nucs, 0, 1) = '';
	}
	;# Now $nucs contains tail of sequence.
	my $iLoopEnd = opt_val("length");
	if ($iTotalNucs < $iLoopEnd) {
	    $iLoopEnd = $iTotalNucs;
	}
	my $iLen;
	for ($iLen = 1; $iLen <= $iLoopEnd; $iLen++) {
	    $ahNucsTable[$iLen - 1]->{substr($nucs, -$iLen, $iLen)}++;
	}
    }

    {
	my @aTmpTable1 = ( );
	my @aTmpTable2 = ( );
	my @aTmpTable3 = ( );
	my $sKey;
	foreach $sKey (keys(%{$ahNucsTable[opt_val("length") - 1]})) {
	    my $iTmp = $ahNucsTable[opt_val("length") - 1]->{$sKey};
	    if (opt_val("mincount") <= $iTmp) {
		my $sTmp = sprintf("%08d %s", $iTmp, $sKey);
		if ($iTmp == 1) {
		    push(@aTmpTable1, $sTmp);
		} elsif ($iTmp == 2) {
		    push(@aTmpTable2, $sTmp);
		} else {
		    push(@aTmpTable3, $sTmp);
		}
	    }
	}
	@aSortedTable = sort {$b cmp $a;} @aTmpTable3;
	push(@aSortedTable, @aTmpTable2);
	push(@aSortedTable, @aTmpTable1);
    }

    if (opt_val("output") eq "f"){
	mkdir ('data', 0777);
	open(TABLEFILE, '>data/' . $filename) || die;
	print TABLEFILE "oligomer,O-value,E-value,";
	my $i;
	for ($i = 1; $i <= opt_val("length") - 2; $i ++){
	    printf TABLEFILE "%d degree Markov,", $i;
	}
	print TABLEFILE "O/E value\n";
    }

    foreach my $sRecord (@aSortedTable) {
	my($iOVal, $sKey) = split(' ', $sRecord);
	my $klen = length($sKey);
	$iOVal =~ s/^0+//;
	my ($order, $iEVal);

	if (opt_val("output") eq "f"){
	    printf TABLEFILE "%s,%d,", $sKey, $iOVal;
	}elsif(opt_val("output") eq "stdout"){
	    &msg_send(sprintf("%s %5d", $sKey, $iOVal));
	}

	if (opt_val("length") == 1){
	    if (opt_val("output") eq "f"){
		printf TABLEFILE "\n";
	    }elsif(opt_val("output") eq "stdout"){
		&msg_send("\n");
	    }
	}else{
	    for ($order = 0; $order <= opt_val("length") - 2; $order++) {
		my $numerator = $iTotalNucs + 1 - $klen;
		my $denominator = 1.0;
		my $offset;
		for ($offset = 0; $offset <= $klen - $order - 1; $offset++) {
		    my $key = substr($sKey, $offset, $order + 1);
		    my $len = length($key);
		    $numerator *=  $ahNucsTable[$len - 1]->{$key} / 
			($iTotalNucs + 1 - $len);
		}
		if (1 <= $order) {
		    for ($offset = 1; $offset <= $klen - $order - 1; $offset++) {
			my $key = substr($sKey, $offset, $order);
			my $len = length($key);
			$denominator *=  $ahNucsTable[$len - 1]->{$key} / 
			    ($iTotalNucs + 1 - $len);
		    }
		} else {
		    $denominator = 1.0;
		}
		if ($denominator <= 0.0) {
		    $iEVal = 0.0;
		} else {
		    $iEVal = $numerator / $denominator;
		}
		
		if (opt_val("output") eq "f"){
		    printf TABLEFILE "%d,", $iEVal if (opt_val("output") eq "f");
		}elsif(opt_val("output") eq "stdout"){
		    &msg_send(sprintf(" %8d", $iEVal));
		}
	    }
	    if (opt_val("output") eq "f"){
		printf TABLEFILE "%.4f\n", $iOVal/$iEVal;
	    }elsif(opt_val("output") eq "stdout"){
		&msg_send(sprintf("   %3.4f\n", $iOVal/$iEVal));
	    }
	    $oe{$sKey} = $iOVal/$iEVal;
	}
    }
    close(TABLEFILE) if (opt_val("output") eq "f");

    return\% oe;

}

General documentation

A. U. Thor, a.u.thor@a.galaxy.far.far.away

perl(1).