#!usr/bin/perl
#paths must be specified for input/output not in active directory

#---->Get a directory of sam files
print "expects single line\nenter the path:\n";
$path=<stdin>;chomp $path;$path=~s/^file:\/\///;
opendir(DIR, $path) or die "can't open $path!";
while (defined($file=readdir(DIR)))	{next if $file =~ /^\.\.?$/;push(@files,$file);} 
@sorted=sort @files;closedir DIR;print "these are the files:\n@sorted\n";
#----<


#----> Provide thresholds for mapping
$threshG = 3;#skips are counted as gap of size 1, so front and back skips combined are a base gap penalty of 2. 
$threshM = 3;#M here is mismatch, not match
$threshL = 3;#how much shorter than full-length a read can be (ie what total skip in bp is allowed). 
#----<


#----> build hash of reference lengths from an input tab-delimited file
open(DICT,"/some/path/lengths.txt") || die "can't open lengths.txt\n";@dict=<DICT>;close DICT; 
foreach $dline (@dict) {chomp $dline;@info=split("\t",$dline);$len{$info[0]}=$info[1];} 
#----<


#loop through batch files
foreach $file (@sorted) {chomp $file;@libname=split('\/',$file);$lib=pop(@libname);

#make output files
open(SAM,">/some/path/$file-mapped.sam") || die "can't make samfile\n";
open(COUNTS, ">/some/path/$file-counts.txt") || die "can't make counts file\n";
open(LOG, ">/some/path/$file-log") || die "can't make logfile\n";

open(FILE,"<$path/$file") || die "file $file not found\n";
#SPECIFY ANY file-specific VARIABLES TO RESET
%string=();%count=();$c=0;$biff=$none=0; 
$badL=$badG=$badM=0;

while (<FILE>)	{if ($_=~m/^@/) {print SAM $_;} else {@arr=split("\t",$_); $c++;
if($arr[2] =~m/^\*/) {$none++;} else {print LOG "\n\n###\nworking on read ",$arr[0],"\n";$voteno=0; 
	if($arr[2]=~/\./) {@r=split('\.',$arr[2]);$ref=$r[0];} else {$ref=$arr[2];} #removes version numbers from accessions

#----> TABULATE GAPS, MATCH LENGTH, FRONT SKIP
#count indels from cigar string, assumes standard bowtie2-output sam format; capture the number of gap positions, skips count as 1

#different steps if S = 0,1, or 2 in order to isolate the cigar string correctly
$cigar=$arr[5];$snum= ($cigar=~tr/S//); #counts number of S characters
$msum=0;$gsum=0;$gap=0;$match=0; #resets match length sum (msum) and gap sum

if($cigar=~/^[0-9]{1,}S/) {$frontskip=$&;$frontskip=~s/S//g;} else {$frontskip=0;} 
$gsum+= ($cigar=~tr/S//);$modcig=$cigar;$modcig=~s/^[0-9]{1,}S//;$modcig=~s/[0-9]{1,}S$//;
$gapcig=$modcig;$gapmax=($modcig=~tr/[ID]//);$gapcount=0;
while($gapcount<$gapmax) {if($gapcig=~/[0-9]{1,}[ID]/) {$gap=$&;$gap=~s/[ID]//g;$gsum+=$gap;$gapcig=~s/[0-9]{1,}[ID]?//;} $gapcount++;} 
$matchcig=$modcig;$mmax=($matchcig=~tr/M//);$mcount=0;print LOG "matchcig is $matchcig and mmax is $mmax\n";
while($mcount<$mmax) {if($matchcig=~/[0-9]{1,}M/) {$match=$&;$match=~s/^[0-9]{1,}[ID]//;$match=~s/M//g;$msum+=$match;$matchcig=~s/^[0-9]{1,}[ID]//;$matchcig=~s/[0-9]{1,}M?//;print LOG "match is $match and matchcig is $matchcig\n";} $mcount++;}

print LOG "cigar is $cigar modcig is $modcig msum is $msum and gsum is $gsum\n";
#----< gap and match tabulation


#----> length check
print LOG "ref length is ",$len{$ref},"\n";
if($len{$ref}-$msum>$threshL) {$voteno++;print LOG "msum is $msum. voting no on LENGTH\n";$badL++;}
#----<


#----> gap check
if($gsum>$threshG) {$voteno++;print LOG "cigar is $cigar gsum is $gsum voting no\n";$badG++;} 
#----<


#----> Global mismatch check  #count mismatches from mdz tag, no reset needed if @mdz always exists, true for standard bowtie2 sam format
$mdz="";
@mdz=split('\:',$arr[18]);$mdz=$mdz[2];$mdz=~s/\^[ACGT]{1,}/\?/g;#deletions in the mdz are simply replaced with a dummy character.
@mdz2=split(/[ACGT]/,$mdz);$mismatch=@mdz2-1;@alt = ($mdz=~/[ACGT]/g);
if($mismatch>$threshM) {$voteno++;print LOG "voting no on threshM\n";$badM++;} 
#----<


#keep or pass read based on thresholds
if($voteno == 0) {$string{$lib}.="@".$ref;print SAM $_;$voteyes++;print LOG "adding alignment to kept sam\n";} else {$biff++;} 
}}}

close FILE;
print COUNTS "ref\tcount\tlib\n";

while (($k,$v)= each %string)	{
print "library is $k\nout of $c reads, biffed on $biff nomatch is $none and kept is ",$c-$biff-$none,"\n";
print "the bad breakdown for G L M was ",$badG," ",$badL," ",$badM,"\n";
%count=();$v=~s/^@//;#leading @ removed;
@hits=split("@",$v);
foreach $h (@hits) {$count{$h}++;}
while (($l,$u)=each %count) {print COUNTS "$l\t$u\t$k\n";}
} 

close OUTFILE;close COUNTS;
}

print "$voteyes matches were kept in total\n";exit;



