extractProject.pl

#!/usr/bin/perl -w

use FindBin;                # Find the script location
use lib "$FindBin::Bin/lib";# Add the script libdir to libs
use Molmed::Sisyphus::Libpath;

use strict;
use Getopt::Long;
use Pod::Usage;
use Cwd qw(abs_path cwd);
use File::Basename;
use File::Find;
use File::Copy qw/cp move/;
use XML::Simple;
use XML::LibXSLT;
use XML::LibXML;

use Molmed::Sisyphus::Common qw(mkpath);
use Molmed::Sisyphus::QStat;
use Molmed::Sisyphus::Plot;

=pod

=head1 NAME

extractProject.pl - Extract all files for delivery to a specific project

=head1 SYNOPSIS

 extractProject.pl -help|-man
 extractProject.pl -runfolder <runfolder> -project <project> -outdir <outdir> [-debug]

=head1 OPTIONS

=over 4

=item -h|-help

prints out a brief help text.

=item -m|-man

Opens the manpage.

=item -runfolder

The runfolder to extract from.

=item -outdir

The directory to put the extracted project in. Defaults to <runfolder>/Projects/<project>

=item -project

The project to extract

=item -skip

Skip delivery of the given lane. Give argument multiple times to skip multiple lanes.
Two reports will be generated, one for deliver and one for usage at the platform.

=item -debug

Print debugging information


=back

=head1 DESCRIPTION

extractProject.pl is a script for extracting, verifying and, optionally,
compressing the files in a runfolder for delivery to a project. The files
to extract are defined by the project option. The files extracted are sequence
files and additional statistics files.

=cut

# Parse options
my($help,$man) = (0,0);
my($rfPath,$outDir,$project) = (undef,undef,undef);
our($debug) = 0;
my @skipLanes;

GetOptions('help|?'=>\$help,
	   'man'=>\$man,
	   'runfolder=s' => \$rfPath,
	   'outdir=s' => \$outDir,
	   'project=s' => \$project,
	   'debug' => \$debug,
	   'skip=i'=> \@skipLanes,
	  ) or pod2usage(-verbose => 0);
pod2usage(-verbose => 1)  if ($help);
pod2usage(-verbose => 2)  if ($man);

unless(defined $rfPath && -e $rfPath){
    print STDERR "Runfolder not specified or does not exist\n";
    pod2usage(-verbose => 1);
    exit;
}

unless(defined $project){
    print STDERR "Project must be specified\n";
    pod2usage(-verbose => 1);
    exit;
}

my %skipLanes;
if(@skipLanes){
  @skipLanes{@skipLanes} = @skipLanes;
}

# Create a new sisyphus object for common functions
my $sisyphus = Molmed::Sisyphus::Common->new(PATH=>$rfPath, DEBUG=>$debug);
$sisyphus->runParameters();
$rfPath = $sisyphus->PATH;
my $statDir = "$rfPath/Statistics/Project_${project}";
my $machineType = $sisyphus->machineType();
my $seqDir = "$rfPath/Unaligned/${project}";

unless(defined $outDir){
    $outDir = "$rfPath/Projects/$project";
}

my $sampleSheet = $sisyphus->readSampleSheet();
unless (exists $sampleSheet->{$project}){
    die "$project does not exists in SampleSheet";
}

my @lanesLeft;
foreach my $lane (keys %{$sampleSheet->{$project}}){
    push @lanesLeft, $lane unless(exists $skipLanes{$lane});
}

if(@lanesLeft > 0){
    die "$statDir does not exist. Did you enter the project correctly?\n" unless(-e $statDir);
    die "$seqDir does not exist. Did you enter the project correctly?\n" unless(-e $seqDir);
}else{
    # If no lanes can be delivered, just skip it.
    mkdir($statDir) unless(-e $statDir);
    print STDERR "No lanes left to process\n";
    system("touch $statDir/extractProject.complete") == 0 or die $!;
    # Calc & store the checksum for the file
    $sisyphus->saveMd5("$statDir/extractProject.complete", $sisyphus->getMd5("$statDir/extractProject.complete", -noCache=>1));
    exit;
}


if(-e "$statDir/extractProject.complete"||-e "$statDir/extractProject.complete.gz"){
    print STDERR "$statDir/extractProject.complete already exists, assuming processing is already complete\n";
    exit 0;
}

# Get the statistics generated by RTA/CASAVA
my($RtaLaneStats,$RtaSampleStats) = $sisyphus->resultStats();

my %checksums;

my $plotter = Molmed::Sisyphus::Plot->new();
my $plotDir = $sisyphus->tmpdir(5e7); # Require 50MB free space
$plotDir .= '/Plots';

#if($debug){
#    foreach my $stat (@sampleStats){
#	print STDERR $stat->sample(), "\n";
#    }
#}

open(my $reportFhDel, '>', "$statDir/reportDel.xml");
print $reportFhDel q(<?xml version="1.0"?>), "\n";
print $reportFhDel q(<?xml-stylesheet type="text/xsl" href="report.xsl"?>), "\n";
print $reportFhDel q(<SequencingReport xmlns="illuminareport.xml.molmed">), "\n";
open(my $reportFhAll, '>', "$statDir/report.xml");
print $reportFhAll q(<?xml version="1.0"?>), "\n";
print $reportFhAll q(<?xml-stylesheet type="text/xsl" href="report.xsl"?>), "\n";
print $reportFhAll q(<SequencingReport xmlns="illuminareport.xml.molmed">), "\n";
my $xs = XML::Simple->new(RootName=>undef);
cp("$FindBin::Bin/report.xsl", "$statDir/report.xsl") or die "Failed to copy '$FindBin::Bin/report.xsl' to '$statDir/report.xsl': $!\n";

# Compile per lane statistics
my %laneStat;
my $laneXmlDataDel = {};
my $laneXmlDataAll = {};
foreach my $lane (keys %{$sampleSheet->{$project}}){ # Lane is second key in sample sheet
    $laneStat{$lane} = undef;
#    $laneXmlDataAll->{Lane}->{$lane}={};
#    unless($skipLanes{$lane}){
#	$laneXmlDataDel->{Lane}->{$lane}={};
#    }

    find({wanted => sub{findLane(\%laneStat, $lane)}, no_chdir => 1}, "$rfPath/Statistics");
    foreach my $read (sort {$a<=>$b} keys %{$laneStat{$lane}}){
	my $stat = $laneStat{$lane}->{$read};
	my %metrics = $stat->metrics();
	$laneStat{$lane}->{$read}->METRICS(\%metrics);
	$metrics{Id} = $metrics{Read};
	delete($metrics{Read});
	delete($metrics{Lane});

	# Add the RTA/CASAVA metrics
	foreach my $key (keys %{$RtaLaneStats->{$lane}->{$read}}){
	    $metrics{$key} = $RtaLaneStats->{$lane}->{$read}->{$key};
	}

	my $plotTitle = $sisyphus->RUNFOLDER . ', Lane ' . $stat->LANE . ', Read ' . $stat->READ;

	my @qplot = $plotter->plotQval($stat,"$plotDir/LanePlots/L00$lane-R$read-Qscores", "Q-score distribution $plotTitle");
	($metrics{QscorePlot} = $qplot[0]) =~ s:^$plotDir/:Plots/:;
	($metrics{QscorePlotThumb} = $qplot[1]) =~ s:^$plotDir/:Plots/:;

	my @aplot = $plotter->plotAdapters($stat,"$plotDir/LanePlots/L00$lane-R$read-Adapters", "Adapter sequences $plotTitle");
	($metrics{AdapterPlot} = $aplot[0]) =~ s:^$plotDir/:Plots/:;
	($metrics{AdapterPlotThumb} = $aplot[1]) =~ s:^$plotDir/:Plots/:;

	my @bplot = $plotter->plotBaseComposition($stat,"$plotDir/LanePlots/L00$lane-R$read-BaseComp", "Base Composition $plotTitle");
	($metrics{BaseCompPlot} = $bplot[0]) =~ s:^$plotDir/:Plots/:;
	($metrics{BaseCompPlotThumb} = $bplot[1]) =~ s:^$plotDir/:Plots/:;

	my @gcplot = $plotter->plotGCdistribution($stat,"$plotDir/LanePlots/L00$lane-R$read-GCdist", "GC Distribution $plotTitle");
	($metrics{GCPlot} = $gcplot[0]) =~ s:^$plotDir/:Plots/:;
	($metrics{GCPlotThumb} = $gcplot[1]) =~ s:^$plotDir/:Plots/:;

	my @lplot = $plotter->plotQ30Length($stat,"$plotDir/LanePlots/L00$lane-R$read-Q30Length", "Q30Length $plotTitle");
	($metrics{Q30Plot} = $lplot[0]) =~ s:^$plotDir/:Plots/:;
	($metrics{Q30PlotThumb} = $lplot[1]) =~ s:^$plotDir/:Plots/:;

	my @dplot = $plotter->plotDuplications($stat,"$plotDir/LanePlots/L00$lane-R$read-Duplicate", "Duplications $plotTitle");
	($metrics{DupPlot} = $dplot[0]) =~ s:^$plotDir/:Plots/:;
	($metrics{DupPlotThumb} = $dplot[1]) =~ s:^$plotDir/:Plots/:;

	my @qpbplot = $plotter->plotQPerBase($stat,"$plotDir/LanePlots/L00$lane-R$read-QvaluePerBase", "Q value per base $plotTitle");
    	($metrics{QValuePerBase} = $qpbplot[0]) =~ s:^$plotDir/:Plots/:;
    	($metrics{QValuePerBaseThumb} = $qpbplot[1]) =~ s:^$plotDir/:Plots/:;

	push @{$laneXmlDataAll->{Lane}->{$lane}->{Read}}, \%metrics;
	unless($skipLanes{$lane}){
	    push @{$laneXmlDataDel->{Lane}->{$lane}->{Read}}, \%metrics;
	}
    }
}
undef(%laneStat);

print $reportFhAll $xs->XMLout($laneXmlDataAll, RootName=>'LaneMetrics', KeyAttr => {Lane => 'Id'});
print $reportFhDel $xs->XMLout($laneXmlDataDel, RootName=>'LaneMetrics', KeyAttr => {Lane => 'Id'});

# Get all the available samples
my @samples;
opendir(my $sDirFh, $statDir) or die "Failed to open '$statDir': $!\n";
foreach my $sampleDir (readdir($sDirFh)){
    next unless (-d "$statDir/$sampleDir");
    opendir(my $saDirFh, "$statDir/$sampleDir") or die "Failed to open '$statDir/$sampleDir': $!\n";
    foreach my $statfile (grep /\.statdump\.zip$/, readdir($saDirFh)){
	# Dual index tags contains a hyphen
        if($statfile =~ m/(.*)_(S\d+|Undetermined|NoIndex|[ACGT]+-?[ACGT]*)_L00(\d)_R(\d)\.statdump.zip/){
            my $sample = $1;
            my $index = $2;
            my $lane = $3;
            my $read = $4;
           $index =~ s/^S//;
           push @samples, {SAMPLE=>$sample,
                           TAG=>$sisyphus->getIndexUsingSampleNumber($lane, $project, $sample, $index, $sampleSheet),
                           LANE=>$lane,
                           READ=>$read,
                           PATH=>"$statDir/$sampleDir/$statfile",
                           SAMPLE_NUMBER=> $index };
        }
    }
    closedir($saDirFh);
}
closedir($sDirFh);

# Compile per sample statistics
my $sampleXmlDataDel = {};
my $sampleXmlDataAll = {};

# assume the same offset for all files
my $offset = 0;

foreach my $sampleFile (sort sortSamples @samples){
    my $sample = Molmed::Sisyphus::QStat->new(DEBUG=>$debug);
    $sample->loadData($sampleFile->{PATH});
    print STDERR "$sample->{SAMPLE}\tLane $sample->{LANE}\tRead $sample->{READ}\n";
    my %metrics = $sample->metrics();
    die "Could not get any metrics\n" unless(%metrics);
    $metrics{Id} = $metrics{Read};
    delete($metrics{Read});
    delete($metrics{Lane});
    my $lane = $sample->{LANE};
    my $smpl = $sample->{SAMPLE};
    my $tag = $sampleFile->{TAG};
    my $read = $sample->{READ};

    # Add the RTA/CASAVA metrics
    foreach my $key (keys %{$RtaSampleStats->{$smpl}->{$lane}->{$read}->{$tag}}){
	$metrics{$key} = $RtaSampleStats->{$smpl}->{$lane}->{$read}->{$tag}->{$key};
    }

    $offset = $sample->OFFSET unless($offset);

    my $barcode = $tag;
    # Entries without barcode are named Undetermined in SampleSheet hash
    if($tag !~ m/\w/){
	$barcode = 'Undetermined';
    }
    $metrics{FragmentSize} = exists($sampleSheet->{$project}->{$lane}->{$barcode}->{FRAGMENT_SIZE}) ?
      $sampleSheet->{$project}->{$lane}->{$barcode}->{FRAGMENT_SIZE} : '';
    $metrics{FragmentLower} = exists($sampleSheet->{$project}->{$lane}->{$barcode}->{FRAGMENT_LOWER}) ?
      $sampleSheet->{$project}->{$lane}->{$barcode}->{FRAGMENT_LOWER} : '';
    $metrics{FragmentUpper} = exists($sampleSheet->{$project}->{$lane}->{$barcode}->{FRAGMENT_UPPER}) ?
      $sampleSheet->{$project}->{$lane}->{$barcode}->{FRAGMENT_UPPER} : '';
    $metrics{LibraryName} = exists($sampleSheet->{$project}->{$lane}->{$barcode}->{LIBRARY_NAME}) ?
      $sampleSheet->{$project}->{$lane}->{$barcode}->{LIBRARY_NAME} : '';

    my $plotTitle = $sample->SAMPLE . '\n' . $sisyphus->RUNFOLDER . ', Lane ' . $sample->LANE . ', Read ' . $sample->READ . ", Tag $tag";
    $plotTitle .= ", Sample Number S" . $sampleFile->{SAMPLE_NUMBER};
    my $tagOrSampleNumber = "S".$sampleFile->{SAMPLE_NUMBER};

    my @qplot = $plotter->plotQval($sample,"$plotDir/$smpl/$tagOrSampleNumber-L00$lane-R$read-Qscores", "Q-score distribution $plotTitle");
    ($metrics{QscorePlot} = $qplot[0]) =~ s:^$plotDir/:Plots/:;
    ($metrics{QscorePlotThumb} = $qplot[1]) =~ s:^$plotDir/:Plots/:;

    my @aplot = $plotter->plotAdapters($sample,"$plotDir/$smpl/$tagOrSampleNumber-L00$lane-R$read-Adapters", "Adapter sequences $plotTitle");
    ($metrics{AdapterPlot} = $aplot[0]) =~ s:^$plotDir/:Plots/:;
    ($metrics{AdapterPlotThumb} = $aplot[1]) =~ s:^$plotDir/:Plots/:;

    my @bplot = $plotter->plotBaseComposition($sample,"$plotDir/$smpl/$tagOrSampleNumber-L00$lane-R$read-BaseComp", "Base Composition $plotTitle");
    ($metrics{BaseCompPlot} = $bplot[0]) =~ s:^$plotDir/:Plots/:;
    ($metrics{BaseCompPlotThumb} = $bplot[1]) =~ s:^$plotDir/:Plots/:;

    my @gcplot = $plotter->plotGCdistribution($sample,"$plotDir/$smpl/$tagOrSampleNumber-L00$lane-R$read-GCdist", "GC Distribution $plotTitle");
    ($metrics{GCPlot} = $gcplot[0]) =~ s:^$plotDir/:Plots/:;
    ($metrics{GCPlotThumb} = $gcplot[1]) =~ s:^$plotDir/:Plots/:;

    my @lplot = $plotter->plotQ30Length($sample,"$plotDir/$smpl/$tagOrSampleNumber-L00$lane-R$read-Q30Length", "Q30Length $plotTitle");
    ($metrics{Q30Plot} = $lplot[0]) =~ s:^$plotDir/:Plots/:;
    ($metrics{Q30PlotThumb} = $lplot[1]) =~ s:^$plotDir/:Plots/:;

    my @dplot = $plotter->plotDuplications($sample,"$plotDir/$smpl/$tagOrSampleNumber-L00$lane-R$read-Duplicate", "Duplications $plotTitle");
    ($metrics{DupPlot} = $dplot[0]) =~ s:^$plotDir/:Plots/:;
    ($metrics{DupPlotThumb} = $dplot[1]) =~ s:^$plotDir/:Plots/:;

    my @qpbplot = $plotter->plotQPerBase($sample,"$plotDir/$smpl/$tagOrSampleNumber-L00$lane-R$read-QvaluePerBase", "Q value per base $plotTitle");
    ($metrics{QValuePerBase} = $qpbplot[0]) =~ s:^$plotDir/:Plots/:;
    ($metrics{QValuePerBaseThumb} = $qpbplot[1]) =~ s:^$plotDir/:Plots/:;

    $metrics{SampleNumber} = "S" . $sampleFile->{SAMPLE_NUMBER};

    push @{$sampleXmlDataAll->{Sample}->{$smpl}->{Tag}->{$tag}->{Lane}->{$lane}->{Read}}, \%metrics;

    unless($skipLanes{$lane}){
	push @{$sampleXmlDataDel->{Sample}->{$smpl}->{Tag}->{$tag}->{Lane}->{$lane}->{Read}}, \%metrics;
    }
}

print $reportFhAll $xs->XMLout($sampleXmlDataAll, RootName=>'SampleMetrics', KeyAttr => {Lane => 'Id', Sample=>'Id', Tag=>'Id'});
print $reportFhDel $xs->XMLout($sampleXmlDataDel, RootName=>'SampleMetrics', KeyAttr => {Lane => 'Id', Sample=>'Id', Tag=>'Id'});

# Collect run information
my %metaData;
$metaData{Project} = $project;
$metaData{RunFolder} = $sisyphus->runfolder();
$metaData{SisyphusVersion} = $sisyphus->version();
$metaData{CsVersion} = $sisyphus->getCSversion();
$metaData{InstrumentModel} = $sisyphus->machineType();
$metaData{RtaVersion} = $sisyphus->getRTAversion();
$metaData{FlowCellVer} = $sisyphus->getFlowCellVersion();
$metaData{FlowCellId} = $sisyphus->fcId();
$metaData{SBSversion} = $sisyphus->getSBSversion();
$metaData{bcl2fastqVersion} = $sisyphus->getBcl2FastqVersion();
$metaData{ClusterKitVersion} = $sisyphus->getClusterKitVersion();

$metaData{Qoffset} = $offset;
my $runInfo = $sisyphus->getRunInfo();
for(my $i=0; $i<@{$runInfo->{reads}}; $i++){
    my $read = $runInfo->{reads}->[$i];
    my $cycles = $read->{last} - $read->{first}; # The last cycle is discarded, so do not add one here
    if($runInfo->{indexed}){
	if($read->{id} != 2){
	    my $id = $read->{id};
	    if($runInfo->{indexed} && $id>1){ # Do not count the index read
		$id -= 1;
	    }
	    push @{$metaData{Read}}, {Id=>$id, Cycles=>$cycles};
	}
    }else{
	push @{$metaData{Read}}, {Id=>$read->{id}, Cycles=>$cycles};
    }
}
print $reportFhAll $xs->XMLout(\%metaData, RootName=>'MetaData');
print $reportFhDel $xs->XMLout(\%metaData, RootName=>'MetaData');

print $reportFhAll q(</SequencingReport>), "\n";
print $reportFhDel q(</SequencingReport>), "\n";

close($reportFhAll);
close($reportFhDel);


# Transform xml + xsl to html
my $xslt = XML::LibXSLT->new();
my $stylesheet = $xslt->parse_stylesheet(XML::LibXML->load_xml(location=>"$statDir/report.xsl", no_cdata=>1));

my $xmlData = stripXmlNameSpace("$statDir/report.xml");

open(my $htmlFhAll, '>', "$statDir/report.html") or die "Failed to open '$statDir/report.html' for writing: $!\n";
print $htmlFhAll
  $stylesheet->output_as_bytes(
			       $stylesheet->transform(
						      XML::LibXML->load_xml(
									    string => $xmlData
									   )
						     )
			      );

close($htmlFhAll);

my $xmlDataDel = stripXmlNameSpace("$statDir/reportDel.xml");

open(my $htmlFhDel, '>', "$statDir/reportDel.html") or die "Failed to open '$statDir/reportDel.html' for writing: $!\n";
print $htmlFhDel
  $stylesheet->output_as_bytes(
			       $stylesheet->transform(
						      XML::LibXML->load_xml(
									    string => $xmlDataDel
									   )
						     )
			      );

close($htmlFhDel);

# Copy the plots from temporary dir into Stats dir
if(-e "$statDir/Plots"){
  $plotter->sysWrap(0, 'rm', '-rf', "$statDir/Plots")==0 or die "Failed to remove old $statDir/Plots";
}
$plotter->sysWrap(0, 'cp', '-r', $plotDir, "$statDir/Plots")==0 or die "Failed to move $plotDir to $statDir/Plots";

# Get the checksums from the files on local disk to circumvent gulo problems
# saveMd5 will force absolute path, so the files must exist also
# in the runfolder, but this should be OK since the cp above
my %checksumCache;
opendir(my $plotDirFh1, $plotDir) or die "Failed to open dir '$plotDir': $!\n";
foreach my $sample (grep /^[^\.]/, readdir($plotDirFh1)){
    next unless(-d "$plotDir/$sample");
    opendir(my $samplePlotDirFh1, "$plotDir/$sample") or die "Failed to open dir '$plotDir/$sample': $!\n";
    foreach my $plot (grep /^[^\.]/, readdir($samplePlotDirFh1)){
	my $checksum = $sisyphus->getMd5("$plotDir/$sample/$plot", -skipMissing=>1, -noCache=>1);
	if(defined $checksum){
	    $sisyphus->saveMd5("$statDir/Plots/$sample/$plot", $checksum);
	    $checksumCache{"Plots/$sample/$plot"} = $checksum;
	}
    }
}

# Copy the README
cp("$FindBin::Bin/README.1", "$statDir/README") or die "Failed to copy '$FindBin::Bin/README.1' to '$statDir/README': $!\n";

# Calculate all checksums in the statdir, just to make sure there are no old cached checksums
$sisyphus->md5Dir($statDir, -noCache=>1, -save=>1);

# All complete
`touch $statDir/extractProject.complete`;
print "Extraction of $project completed\n";

sub sortSamples{
    return($a->{SAMPLE} . $a->{LANE} . $a->{TAG} . $a->{READ} cmp $b->{SAMPLE} . $b->{LANE} . $a->{TAG} . $b->{READ});
}

sub findLane{
    my $laneStat = shift;
    my $lane = shift;
    my $file = $_;
    if($file =~ m/_L00${lane}_R\d\.statdump.zip/){
#	print STDERR "Adding $file\n";# if($debug);
	my $stat = Molmed::Sisyphus::QStat->new(DEBUG=>$debug);
	$stat->loadData($file);

	my $read = $stat->READ();
#	print $stat->{SAMPLE}, " ", $stat->{Q30LENGTH}->{0}, "\n" if($read==1);

	if(defined $laneStat->{$lane}->{$read}){
	    my $tmp = $laneStat->{$lane}->{$read}->add($stat);
	    $laneStat->{$lane}->{$read} = $tmp;
	}else{
	    $laneStat->{$lane}->{$read} = $stat;
	}
    }
}

sub stripXmlNameSpace{
    my $file = shift;
    # Strip the namespace from the xml data, adding it to the xsl is a mess
    $/='';
    open(my $xmlFh, $file) or die;
    my $xmlData = <$xmlFh>;
    close($xmlFh);
    $/="\n";
    $xmlData=~ s/xmlns="illuminareport.xml.molmed"//;
    return $xmlData;
}