#!/usr/bin/perl  

use FindBin qw($Bin);
use lib "$Bin";
use Getopt::Long;
use Cwd;
use Cwd 'abs_path';
use Storable;
use File::Basename;
use Spiders::Dta;
use Spiders::Params;
use Spiders::ProcessingRAW;
use Spiders::ProcessingMzXML;
use Spiders::Job;
use Spiders::Path;
use Parallel::ForkManager;
use Spiders::PeakDetection;
use Spiders::Decharge;
use Spiders::MS2_scoring;
use Spiders::MassCorrection;
use Spiders::MathUtils;
use Spiders::Summary;


my $library = $Bin;
my $R_library = $Bin . "/R-3.1.0/";
if(!-e($R_library))
{
	print "\n\n  Please build the R (version: 3.1.0) under $library \n\n";
	exit;
}
my $VERSION = 1.02;

my $progname = $0;
$progname =~ s@(.*)/@@i;

my ($help,$parameter,$raw_file);
GetOptions('-help|h'=>\$help,
			'-p=s'=>\$parameter,
		);

if(!-e ($parameter))
{
	print "please input the paramter file\n\n";
	exit;
}	
usage() if ($help || !defined($parameter));

my $log_temp = time();
open($LOG,">./$log_temp");

print "\n\n  Initializing JUMPm program (version $VERSION)\n\n";
print $LOG "\n\n  Initializing JUMPm program (version $VERSION)\n\n";
my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst)=localtime(time); 

my $p = Spiders::Params->new('-path'=>$parameter);
my $params=$p->parse_param();
$p->parameter_check($params);

## Create the path for multiple raw files
my %rawfile_hash;
print "  Using the following rawfiles:\n";
print $LOG "  Using the following rawfiles:\n";
foreach $arg (sort @ARGV)
{
    my @suffixlist=();
    push @suffixlist,".raw";
    push @suffixlist,".RAW";
    push @suffixlist,".mzXML";

    if($arg=~/.[raw|RAW|mzXML]/)
    {
		print "  $arg","\n";
		print $LOG "  $arg","\n";		
	}
}

foreach $arg (sort @ARGV)
{
    my @suffixlist=();
    push @suffixlist,".raw";
    push @suffixlist,".RAW";
    push @suffixlist,".mzXML";

    if($arg=~/.[raw|RAW|mzXML]/)
    {	

		my ($filename, $directory, $suffix) = fileparse($arg,@suffixlist);		
        system(qq(mkdir $directory/$filename >/dev/null 2>&1));
        system(qq(mv $arg $directory/$filename >/dev/null 2>&1));
        my $datafile = "$directory/$filename";
		if(!-e($datafile))
		{
			print "  There is no such raw file ($arg)!\n";
			exit;
		}		
        my $path = new Spiders::Path($datafile);

        my $list = $path->make_directory_list();
        if(@$list)
        {
            $newdir = $path->choose_dir_entry($list,"  Choose a .out file directory",$newdir);
        }
        else
        {
			$newdir	= $filename . ".1";
          #  $newdir = $path->ask("  Choose a .out file directory",$newdir);
        }
        print "  Using: $newdir\n";
        print $LOG "  Using: $newdir\n";		
        $path->add_subdir($newdir);
        my $dir =  $path->basedir() . "/$newdir";
		my $rawfile = "$datafile/$arg";
		$rawfile_hash{$rawfile} = $dir; 
    }
}

foreach my $raw_file(sort keys %rawfile_hash)
{		
	print "\n  Searching data: $raw_file\n";
	print $LOG "\n  Searching data: $raw_file\n";	
	my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst)=localtime(time); 
	print "  Start: ";
	print $LOG "  Start: ";	
	printf "%4d-%02d-%02d %02d:%02d:%02d\n",$year+1900,$mon+1,$mday,$hour,$min,$sec;	
	printf $LOG "%4d-%02d-%02d %02d:%02d:%02d\n",$year+1900,$mon+1,$mday,$hour,$min,$sec;	
	
	my $curr_dir = getcwd;
	my $dta_path = $rawfile_hash{$raw_file};
	  
	if($raw_file =~/\.mzXML/)
	{
		$raw_file =~s/\.mzXML//g;
	}

	my $proc_raw = new Spiders::ProcessingRAW();
	$proc_raw->set_log_file($LOG);		
	$proc_raw->set_raw_file($raw_file);

###################### window part ##########################
	print "  Converting .raw into .mzXML file\n";
	print $LOG "  Converting .raw into .mzXML file\n";
	my $proc_xml = new Spiders::ProcessingMzXML();	
	$proc_xml ->set_dta_path($dta_path);
	$proc_xml ->set_log_file($LOG);	
	my (%ms_hash,%msms_hash,@mz_array);		
	my $MS1_peaks_sequenced;
	my ($msms_hash_corrected,$mz_array_corrected,$meanMassShift, $stdMassShift,)=(0,0,0,0);	
	if($params->{data_acquisition_mode} == 1)
	{	
		my $mzXML = $proc_raw->raw2mzXML("centroid");
	
		print "\n\n  Using centroid mode\n";
		print $LOG "\n\n  Using centroid mode\n";		
		print "  Extracting MS1 and MS2 peaks from .mzXML\n";
		print $LOG "  Extracting MS1 and MS2 peaks from .mzXML\n";	
		$proc_xml ->set_mzXML_file($mzXML);

################### preprocessing #########################
		$proc_xml->set_parameter($params);
		$proc_xml->generate_hash_dta(\%ms_hash, \%msms_hash, \@mz_array, $params);
		my $ms1N = scalar(keys %{$ms_hash{'surveyhash'}});
		my $ms2N = scalar(keys %msms_hash)-scalar(keys %{$ms_hash{'surveyhash'}});
		
		printf("\n  There are %d MS and %d MS/MS in the entire run\n", $ms1N , $ms2N);
		printf $LOG ("\n  There are %d MS and %d MS/MS in the entire run\n", $ms1N , $ms2N);

		if($params->{mass_correction})
		{
			print "\n\n  Mass correction\n";
			print $LOG "\n\n  Mass correction\n";	
			my $masscorr = new Spiders::MassCorrection();
			$masscorr->set_log_file($LOG);
			### No correction for MS1, but it will perform correction at the database search step
			($msms_hash_corrected,$mz_array_corrected,$meanMassShift, $stdMassShift) = $masscorr->massCorrection(\%ms_hash, \%msms_hash, \@mz_array, $params);
			printf $LOG ("  Calculated mass-shift: mean = %.5f ppm and SD = %.5f ppm\n", $meanMassShift, $stdMassShift); 
			%msms_hash = %$msms_hash_corrected;
			@mz_array = @$mz_array_corrected;	
		}
		
		my $decharge = new Spiders::Decharge();
		$decharge->set_parameter($params);
		$decharge->decharge(\%ms_hash, \%msms_hash, \@mz_array);

		print "  Generating MS1 dta files\n";	
		print $LOG "  Generating MS1 dta files\n";			
		$proc_xml->generate_MS1_file(\@mz_array);
		$MS1_peaks_sequenced = $proc_xml->generate_MS2_file(\%msms_hash,\@mz_array);		
	}
	elsif($params->{data_acquisition_mode} == 2)
	{
		my $mzXML = $proc_raw->raw2mzXML("profile");
		print "\n\n  Using profile mode\n";	
		print $LOG "\n\n  Using profile mode\n";		
		print "  Feature detection\n";
		print $LOG "  Feature detection\n";	


		$proc_xml ->set_mzXML_file($mzXML);

		$proc_xml->set_parameter($params);
		print $LOG "  Gathering scan information\n";			
		$proc_xml->generate_hash_dta(\%ms_hash, \%msms_hash, \@mz_array, $params);
		my $ms_hash_peak_detection = $proc_xml->{'_mshash4peakdetection'};
		my $ms1N = scalar(keys %{$ms_hash{'surveyhash'}});
		my $ms2N = scalar(keys %msms_hash)-scalar(keys %{$ms_hash{'surveyhash'}});
		
		printf("\n  There are %d MS and %d MS/MS in the entire run\n", $ms1N , $ms2N);
		printf $LOG ("\n  There are %d MS and %d MS/MS in the entire run\n", $ms1N , $ms2N);

		my $peakdetection = new Spiders::PeakDetection();
		$peakdetection->set_log_file($LOG);
		my ($ms,$peak) =$peakdetection->Gen_3Dpeaks($params,$mzXML,$dta_path);
		
		my $correction_raw = $peakdetection->generate_correction_file($peak,$dta_path);
		
		if($params->{mass_correction})
		{
			print "\n  Mass correction\n";
			print $LOG "\n  Mass correction\n";	
			my $masscorr = new Spiders::MassCorrection();
			my $mass_shift_file = $dta_path . ".mass_shift";
			$masscorr->set_log_file($LOG);
			($ms_hash_corrected,$msms_hash_corrected,$mz_array_corrected,$mass_shift_mean_all) = $masscorr->massCorrection(\%ms_hash, \%msms_hash, \@mz_array,$correction_raw, $params,$mass_shift_file);
			@mz_array = @$mz_array_corrected;
			#$peakdetection->generate_MS1_file($ms,$peak,$dta_path);
			$peakdetection->generate_MS1_file_based_feature($peak,$dta_path);
			$MS1_peaks_sequenced = $proc_xml->generate_MS2_file($msms_hash_corrected,\@mz_array);
			store $mass_shift_mean_all, "$dta_path/.mass_shift";
		}
		else
		{
			#$peakdetection->generate_MS1_file($ms,$peak,$dta_path);
			$peakdetection->generate_MS1_file_based_feature($peak,$dta_path);
			$MS1_peaks_sequenced = $proc_xml->generate_MS2_file(\%msms_hash,\@mz_array);		
		}
		print "\n";
		print $LOG "\n";

		
	}
	else
	{
		system(qq(cp $log_temp $dta_path/JUMPm.log));		
		print " please set the right data_acquisition_mode\n";
		exit;
	}
		
	my @dtafiles = glob("$dta_path/*.MS1");
	my $dta=new Spiders::Dta; 
	
	my $job = new Spiders::Job;
	$job->set_log_file($LOG);
	$job->set_library_path($library);
	$job->set_R_library_path($R_library);		
	$job->set_dta_path($dta_path);
	$job->set_parameter($params);
	
	my $summary = new Spiders::Summary();	
	$summary->set_library_path($R_library);
	$summary->set_log_file($LOG);
	$summary->generate_feature_table($dta_path);
	if($params->{'labeled_data'})
	{
		print "  You are searching a labeled data set\n";
		print $LOG "  You are searching a labeled data set\n";		
		$job->create_pscript;
		my $random = int(rand(100));
		my $ratio;
		my $defect;
		$ratio->{'NC'} = 0;
		$ratio->{'CC'} = 0;
		$parameter = abs_path($parameter);			
		if($params->{loading_normalization})
		{		
			print "  Estimating the labeling mix ratio (it takes time and please be patient!)\n";
			print $LOG "  Estimating the labeling mix ratio (it takes time and please be patient!)\n";
	
			$job->runjobs(\@dtafiles,$dta_path,"pair_${random}",0,0,$parameter);		
			($ratio,$defect) = $summary->calculate_intensity_ratio($dta_path);
			my $NC_ratio = 2**$ratio->{'NC'};
			my $CC_ratio = 2**$ratio->{'CC'};
			
			print "  Intensity ratio between N15 and C12 is $NC_ratio \n";
			print "  Intensity ratio between C13 and C12 is $CC_ratio \n";
			print $LOG "  Intensity ratio between N15 and C12 is $NC_ratio \n";
			print $LOG "  Intensity ratio between C13 and C12 is $CC_ratio  \n";
			print "  Mass defect between N15 and C12 is $defect->{'NC'}->[0] \n";
			print "  Mass defect between C13 and C12 is $defect->{'CC'}->[0] \n";
			print $LOG "  Mass defect between N15 and C12 is $defect->{'NC'}->[0] \n";
			print $LOG "  Mass defect between C13 and C12 is $defect->{'CC'}->[0] \n";
		}
		else
		{
			print "  No loading normalization\n";
			print $LOG "  Not loading normalization\n";
			$ratio->{'CC_std'} = 1;			
			$ratio->{'NC_std'} = 1;			
			$defect->{'NC'}->[0]=1.00335;
			$defect->{'CC'}->[0]=0.99703;
			$defect->{'NC'}->[1]=1;
			$defect->{'CC'}->[1]=1;			
		}
		print "\n  Pair scoring, searching formulas and structures (it takes time and please be patient!)\n";
		print $LOG "\n  Pair scoring, searching formulas and structures (it takes time and please be patient!)\n";				
		$job->create_mscript;
		$random = int(rand(100));
			
		$job->runjobs(\@dtafiles,$dta_path,"sch_${random}",$ratio,$defect,$parameter);
		print "\n\n  Summarizing pair results\n";
		print $LOG "\n\n  Summarizing pair results\n";		
		$summary->pair_summary($dta_path);
		print "\n  Summarizing missile results";
		print $LOG "\n  Summarizing missile results";	
		my $ms_hash_mol = $summary->missile_summary($dta_path);
		my $missile_scan = scalar keys (%$ms_hash_mol);
		
		if($missile_scan > 0)
		{	
			print "\n  Performing MS2 matching and scoring\n\n";
			print $LOG "\n  Performing MS2 matching and scoring\n\n";	

			my ($file_array,$MS1_MS2_matched,$ms2_scan_num)=$job->luanch_ms2_jobs(\%msms_hash,$ms_hash_mol,$dta_path,$parameter);
			

			print "\n  $ms2_scan_num MS2 scans match to MISSILES (one missile may have multiple MS2 matches)\n";
			print $LOG "\n  $ms2_scan_num MS2 scans match to MISSILES (one missile may have multiple MS2 matches)\n";
			if($ms2_scan_num>0)
			{
				$random = int(rand(100));
				$job->runjobs($file_array,$dta_path,"frag_${random}",0,0,$parameter);

######## Finally summary #############
				my $missile_structure = $summary->final_summary($ms_hash_mol,$MS1_MS2_matched,$dta_path);
			}
		}		
	}
	else	
	{
		$proc_xml->generate_MS1_file_unlabel($MS1_peaks_sequenced);		
		print "\n  Searching formulas and structures (it takes time and please be patient!)\n";		
		print $LOG "\n  Searching formulas and structures (it takes time and please be patient!)\n";			
		$job->create_mscript_unlabel();
		$random = int(rand(100));
		$parameter = abs_path($parameter);			
		$job->runjobs(\@dtafiles,$dta_path,"sch_${random}",0,0,$parameter);	
		my $ms_hash_mol = $summary->missile_summary_unlabel($dta_path);
		my ($file_array,$MS1_MS2_matched,$ms2_scan_num)=$job->luanch_ms2_jobs(\%msms_hash,$ms_hash_mol,$dta_path,$parameter);
		print "\n  $ms2_scan_num MS2 scans match to MISSILES (one missile may have multiple MS2 matches)\n";
		print $LOG "\n  $ms2_scan_num MS2 scans match to MISSILES (one missile may have multiple MS2 matches)\n";
		
		$random = int(rand(100));
		$job->runjobs($file_array,$dta_path,"frag_${random}",0,0,$parameter);

######## Finally summary #############
		my $missile_structure = $summary->final_summary($ms_hash_mol,$MS1_MS2_matched,$dta_path);
	}
	($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst)=localtime(time); 
	print "  JUMPm finished: ";
	print $LOG "  JUMPm finished: ";	
	printf "%4d-%02d-%02d %02d:%02d:%02d\n\n",$year+1900,$mon+1,$mday,$hour,$min,$sec;	
	printf $LOG "%4d-%02d-%02d %02d:%02d:%02d\n\n",$year+1900,$mon+1,$mday,$hour,$min,$sec;	
	system(qq(cp $log_temp $dta_path/JUMPm.log));
	system(qq(mkdir $dta_path.intermediate));

	system(qq(mv $dta_path.pair.all.comb $dta_path.intermediate/));	
	system(qq(mv $dta_path.pair.all.comb.uniq $dta_path.intermediate/));
	system(qq(mv $dta_path.tmp.feature $dta_path.intermediate/));
	system(qq(mv $dta_path.missile $dta_path.intermediate/));
	system(qq(mv $dta_path.missile.unique.pass $dta_path.intermediate/));
	system(qq(mv $dta_path.MS2.structure $dta_path.intermediate/));

}

sub usage {

print <<"EOF";
	
################################################################
#                                                              #
#       **************************************************     # 
#       ****                                          ****     # 
#       ****                 JUMPm                    ****     # 
#       ****        version $VERSION                  ****     # 
#       ****        Xusheng Wang / Junmin Peng        ****     # 
#       ****         Copyright (C) 2012 - 2013        ****     # 
#       ****            All rights reserved           ****     # 
#       ****                                          ****     # 
#       **************************************************     # 
#                                                              #
################################################################

Usage: $progname -p parameterfile rawfile.raw 
	or
       $progname -p parameterfile rawfile.mzXML
	
EOF
exit 1;
}
