###############################################################################
# Name: min li
# Group:Clough Lab
# Date: 5/24/2005
# usage: this program is to caculate the average intensity of each sample based on the design file which is for R normalization, and the imputed file for SAS.
#         output file is a .txt file, the program will take design file and imputedfile for sas. it will generate a file which contains Genbank location,
#         ID each sample average intensity and the number of sample intensity used for average.
#
# Executable format: Perl caculatingAverageIntensity.pl
# Note: design file has 3 columns: array dye and sample in order.
#          file for sas analysis contains 4 columns of location information and GB ID.
#     11/29/05: stdev and sd/average across all samples for a gene is added.
###############################################################################
if($ARGV[0])
{
  die("\n!!!!!WARNING: you shouldn't have any parameter after the program name!!!!!\n");
}
############
while(1)
{
    print "\n\nPlease Enter the Design File Name which you used for Rmaanova Normalization\n";
    print "R_MAANOVA Design File(Extension is .txt)\n\n";
    $designfile=<STDIN>;
#check the output file formatf
    if($designfile=~m/.+\.txt$/)
    {
             last;
    }
}
##########################################
############
while(1)
{
    print "\n\nPlease Enter the File Name which you used for SAS Analysis\n";
    print "SAS Analysis File(Extension is .txt)\n\n";
    $fileforsas=<STDIN>;
#check the output file formatf
    if($fileforsas=~m/.+\.txt$/)
    {
             last;
    }
}
##########################################
############
while(1)
{
    print "\n\nPlease Enter the Average Intensity file\n";  #this is output file
    print "Average Intensity File(Extension is .txt)\n\n";
    $averageintensiyfile=<STDIN>;
#check the output file formatf
    if($averageintensiyfile=~m/.+\.txt$/)
    {
             last;
    }
}
##########################################
open DESIGNFILE, $designfile or die "could not open $designfile for reading\n";
$firstline=<DESIGNFILE>; ####this is from design file reading.
open OUTFILE, ">$averageintensiyfile" or die "could not open $averageintensiyfile for writing\n";

@arrayofsamples; #hold all the samples in an array in a linear style.
@arrayofarrays; #hold all the array number from design file in a linear style.
$linenumber=0;
@linenumberarray;
$howmanysample=0;
$emptyarrayofsample=0;

while($line=<DESIGNFILE>) #this is from second line because first line already was read.
{
	chomp $line;
	$linenumber++;
	push @linenumberarray, $linenumber;

	($array, $dye, $sample)=split/\t/, $line;

	$foundsamesample=0;

	if($emptyarrayofsample==1)
	{
		for($i=0; $i<@arrayofsamples; $i++)
		{
			if($sample eq $arrayofsamples[$i])
			{
				$foundsamesample=1;
				last;
			}
		}
	}
	else # this is meaning the arrayofsample is empty. actually only work for the first line or first sample id.
	{
		$emptyarrayofsample=1;
	}
	push @arrayofarrays, $array;
	push @arrayofsamples, $sample;
	if($foundsamesample==0)
	{
		$howmanysample++;
	}
}
#print "\nhow many samples: $howmanysample\n";
#####################################################################
for($i=0; $i<@linenumberarray; $i++)
{
	if($arrayofsamples[$i] ne "gone")
	{
		if($i==0)
		{
			$samplearraysindex=$i;
		}
		else
		{
			$samplearraysindex++;
		}
		push @{$samplearrays[$samplearraysindex]}, $linenumberarray[$i];
	}
	for($j=$i+1; $j<@arrayofsamples; $j++)
	{
		if(($arrayofsamples[$i] eq $arrayofsamples[$j]) and ($arrayofsamples[$j] ne "gone"))
		{
			push @{$samplearrays[$samplearraysindex]}, $linenumberarray[$j];
			$arrayofsamples[$j]="gone";
		}
	}
}
###due to the changing to "gone" in arrayofsamples, no redundancy samples left. here we build a title for later output file.
$sampletitle="";
foreach $item (@arrayofsamples)
{
	if($item ne "gone")
	{
		if($sampletitle eq "")
		{
			$sampletitle="Average_Sample".$item."\t"."NumberofSampleUsed";
		}
		else
		{
			$sampletitle=$sampletitle."\t"."Average_Sample".$item."\t"."NumberofSampleUsed";
		}
	}
}
#print $sampletitle;
#print "\n";
####################################################################
for($i=0; $i<$howmanysample; $i++)
{
      #	print join("\t", @{$samplearrays[$i]});
       #	print "\n";
}

####################################################################
open FILEFORSAS, $fileforsas or die "could not open $fileforsas for reading\n";
$firstline=<FILEFORSAS>; ####this is from sas input file reading.
#print $firstline,"\n\n";
#$firstline=~m/(\D+\t\D+\t\D+\t\D+\t\S+)(\t)(.+)/;
#print "first is $1\n";
$titleforAverage_SD="\taverage_across_allSlides\tStandardDeviation\tSD/Ave";
$finaltitle="Meta Row\tMeta Column\tGrid Row\tGrid Column\tGB ID\t".$sampletitle.$titleforAverage_SD;
@array_average_foronegene;
print OUTFILE $finaltitle, "\n";

while($line=<FILEFORSAS>)
{
	chomp $line;
	$line=~m/(\S+\t\S+\t\S+\t\S+\t\S+)(\t)(.+)/;
	$temp=$3;
        chomp $temp;
        $temp=~s/^\t//;
	@onelinerealdata=split/\t/, $temp;
      #  print "@onelinerealdata\n";
       # print "$temp\n";
	$inputforaveragefile=$1;
	for($i=0; $i<$howmanysample; $i++)
	{
		$numberofsampleused=0;
		$sum=0;
	      
        foreach $item (@{$samplearrays[$i]})
		{
			if($onelinerealdata[$item-1] ne ".")
			{
                            #	print "item is $onelinerealdata[$item-1]\n";
				$sum+=2**$onelinerealdata[$item-1];  #change back to the intensity
				$numberofsampleused++;
			}
		}
		if($numberofsampleused==0)
		{
			$average=0;
		}
		else
		{
			$average=$sum/$numberofsampleused;
		}
		$average=int $average;      ##################   this convert $average to an integer.
		#$Sumallsamplefor_onegene=$Sumallsamplefor_onegene + $average;
		push @array_average_foronegene, $average;
		$inputforaveragefile=$inputforaveragefile."\t".$average."\t".$numberofsampleused;
	}
	
	print OUTFILE $inputforaveragefile, cal_average_sd(@array_average_foronegene), "\n";
	
	@array_average_foronegene=();  #empty the array for processing the next line;
	
}
##################################################################################
sub cal_average_sd   #####find the stdev for samples for a gene.
{
	my @averages=@_;
	my $sum_averages=0;
	my $average_allsample_onegene;
	my $forsqure;
	my $sd;
	my $sd_average;
	foreach $item (@averages)
	{
		$sum_averages=$sum_averages + $item;
	}
	$average_allsample_onegene=int $sum_averages/($#averages+1);
	foreach $item (@averages)
	{
		$forsqure=$forsqure + ($item-$average_allsample_onegene)**2;
	}
	$sd=($forsqure/($#averages+1))**0.5;
	$sd_average=$sd/$average_allsample_onegene;
	return "\t$average_allsample_onegene\t$sd\t$sd_average";
}



