Untitled

#!/usr/bin/perl -w
#
#	egq.pl - EfficientGenomeQuery
#	copyright Erich Blume, 2006 
#	Erich Blume - eblume@ucsc.edu
#
$|=1;

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.
foreach $genfilename (@genfiles) { #every file in the list
	$start = 0;
	$lastlinelength = 0;
	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();
		$_=uc($_);
		if (!(/^(>|;)/)) { #omit FASTA comments and headers
			if ($start == 0) { 
				$buffer=$_;
				$start++;
				$lastlinelength=length($_);
			} else {
				$buffer.=$_;
				for my $pos (0..length($buffer)-20) {
					my $query = substr($buffer,$pos,20);
					#DebugHit
					#print "${query}\n";
					if ($query=~/^[ATGC]{20}$/) {
						if(exists($smallquerylist{$query})) {
################Re-Indent################################
foreach $big (@{$smallquerylist{$query}}) {
	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############################
						}
					}
				}
				substr($buffer,0,$lastlinelength) = ''; 
				$lastlinelength = length($_);
			}
		}
	}
	close(GENFILE);
}

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

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