Symian/EGQ

#!/usr/bin/perl -w
#
#	egq.pl - EfficientGenomeQuery
#	copyright Erich Blume, 2006
#

use FindBin qw($Bin);
use lib "$Bin/../lib";


 
# preProcessQuery ($line) - determines $line is a valid query and 
#                           returns it with any necessary conversions.
sub preProcessQuery {
	my $line=uc(shift); #converts lowercase to uppercase
	chomp($line);
	if ($line=~/^[ATGC]{25}$/) { #Valid queries are 25 A,T,G,C only!
		return $line;
	} else {
		print STDERR "Query '${line}' at line $. of $ARGV ".
		"is invalid! Skipping...\n";
	}
}

# createSmallQueries ($line) - chop valid 25-mer into 6 20-mers.
#
sub createSmallQueries {
	my $line=shift;
	return (substr($line,0,20),
		substr($line,1,20),
		substr($line,2,20),
		substr($line,3,20),
		substr($line,4,20),
		substr($line,5,20)
		);
}

##TODO: Uncomment for production release! PROCESSES ~900MB TEXT
#@genfiles=("chr1.fa","chr2.fa","chr3.fa","chr4.fa","chr5.fa","chr6.fa",
#	   "chr7.fa","chr8.fa","chr9.fa","chr10.fa","chr11.fa","chr12.fa",
#	   "chr13.fa","chr14.fa","chr15.fa","chr16.fa","chr17.fa","chr18.fa",
#	   "chr19.fa","chr20.fa","chr21.fa","chr22.fa","chrX.fa","chrY.fa");

#TODO: Delete debug mode, below.
@genfiles=("test.fa");

#####################MAIN#############################


###Begin parse of query file(s).

while (<>) {
	$_=preProcessQuery($_);
	if (exists($bigquerylist{$_})) {
		print STDERR "Query $_ at line $. of $ARGV already ".
		"exists! Skipping...\n";
	} else {
		$bigquerylist{$_}=0;
		@queries=createSmallQueries($_);
		foreach $querysmall (@queries) {
			$smallquerylist{$querysmall} = [] 
				unless exists $smallquerylist{$querysmall};
			push @{$smallquerylist{$querysmall}}, $_;
		}
	}
}


###DebugDumps:
#foreach $value (%bigquerylist) { print "${value}\n"}
#foreach $aref (%smallquerylist) {print "@{$aref}\n"}

###Begin peicewise parse of chromosomes.
$query = "";

foreach $genfilename (@genfiles) { #every file in the list
	##TODO: Make sure length always returns num of characters. 
	print "Now opening ".$Bin."/HumGenome/".$genfilename."\n";
	open(GENFILE, $Bin."/HumGenome/".$genfilename)
		or die("Bad file: ".$Bin."/HumGenome/".$genfilename."\n");
	while (<GENFILE>) { #every line in the file
		chomp();
		if (!(/^(>|;)/)) { #omit FASTA comments and headers
			foreach $mer (split //) { #every letter in the line
				$mer=uc $mer;
				if ((length($query) < 20) && $mer=~/[ATGC]/) {
					$query=$query . $mer;
				} elsif ($mer!~/[ATGC]/){
					#TODO: error facility for irregular FASTA values.
				} elsif (length($query) > 20) {
					die ("Fatal Error: query $query is too long!\n");
				} else {
					
					if (exists($smallquerylist{$query})) {
						$bigref = $smallquerylist{$query};
						foreach $big (@{$bigref}) {
############################Re-Indent####################

if (exists($bigquerylist{$big})) {
	#DebugHit
	print "I found a match for query $query in 25-mer $big !\n";
	$bigquerylist{$big}++;
} else { #shouldn't ever happen, hence the die
	die("Can't find 25-mer query $big that the 20-mer $query says should exist!\n");
}
###########################End Re-Indent#################
						}
					}

					if ($mer=~/[ATGC]/) {
						$query = substr($query,1,19) . $mer;
					} else {
						$query = substr($query,1,19);
						#TODO: error facility for bad FASTA stuff.
					}
				}
			}
		}
	}
	close(GENFILE);
	$query = ""; #Don't wrap searches around chromosomes! It just doesn't happen.
}

print "\nLookup completed.\nQuery: Hits\n";

foreach $key (keys %bigquerylist) {
	$value = $bigquerylist{$key} / 6;
	print "${key}: ${value}\n";

}