G::Seq

Markov

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

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

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 description

Methods code

    my $self = shift;

    &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;

    my $pkg = shift;
    my $filename = shift;
    my $option = shift;
    my $this;

    return $this;

General documentation

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

perl(1).