quickReport.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 File::Basename;

use Molmed::Sisyphus::Common qw(mkpath);
use Molmed::Sisyphus::QStat;
=pod

=head1 NAME

quickReport.pl - Create a brief report and mail it to the specified address

=head1 SYNOPSIS

 quickReport.pl -help|-man
 quickReport.pl -runfolder <runfolder> -mail <email address> [-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 generate report on.

=item -mail

Send the report to this email address

=item -debug

Print debugging information

=back

=head1 DESCRIPTION

Compiles some quick statistics about the run and sends a mail to the specified address.

=cut

# Parse options
my($help,$man) = (0,0);
my($rfPath,$mail) = (undef,undef);
my $sender = 'seq@medsci.uu.se';
our($debug) = 0;

GetOptions('help|?'=>\$help,
	   'man'=>\$man,
	   'runfolder=s' => \$rfPath,
	   'mail=s' => \$mail,
       'sender=s' => \$sender,
	   'debug' => \$debug,
	  ) 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;
}

my $OFFSET = 33;

# Create a new sisyphus object for common functions
my $sisyphus = Molmed::Sisyphus::Common->new(PATH=>$rfPath, DEBUG=>$debug);
$rfPath = $sisyphus->PATH;
my $ignoreFastqs = 0; 

my $config = $sisyphus->readConfig();

if(defined($config->{IGNORE_MISSING_FASTQS})){
    $ignoreFastqs = $config->{IGNORE_MISSING_FASTQS};
}

# Initialize the random generator and use a random seed based upon the flowcell id
my $fcid = $sisyphus->fcId() || "ABC123CXX";
my $rseed = 0;
foreach my $c (split //, $fcid) {
    $rseed += ord($c);
}
srand($rseed);
print STDERR "Using random seed $rseed based on fcId $fcid\n";

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

my $sampleSheet = $sisyphus->readSampleSheet();
my $machineType = $sisyphus->machineType();

my %files;
my $baseQC;
sub findFastq{
    my $files = shift;
    my $file = $_;
    if($file !~ /Undetermined_indices/ && $file !~ /\/Data\// && $file !~ /\/IndexCheck\//){
        if($file =~ m/\.fastq(\.gz)?$/){
            my @path = split '/', $file;
            my $project = $path[-3];
            $project =~ s/^Project_//;
            my $sample = $path[-2];
            $sample =~ s/^Sample_//;
            if($file =~ m/.*\/(.+)_([ACTG]+-?[ACGT]*|NoIndex)_L(\d{3})_R(\d)_(\d{3})\.fastq\.gz$/){
                my ($sample,$index,$laneId,$read,$segment) = ($1,$2,$3,$4,$5);
                $index =~ s/NoIndex/Undetermined/ if($index =~ /NoIndex/);
                $files->{$laneId}{$sample}{$index}{$read}{$segment} = $file;
            }
            elsif($file =~ m/.*\/(.+)_S(\d*)_L(\d{3})_R(\d)_(\d{3})\.fastq\.gz$/){
                my ($sample,$index,$laneId,$read,$segment) = ($1,$2,$3,$4,$5);
                $sample = "Sample_" . $sample;
                $files->{$laneId}{$sample}{$index}{$read}{$segment} = $file;
            }
        }
    }
}
use File::Find;
find({wanted => sub{findFastq(\%files)}, no_chdir => 1, follow => 1}, "$rfPath/Unaligned");

my $laneQC;
my $samples;

foreach my $proj (keys %{$sampleSheet}){
    foreach my $lid (keys %{$sampleSheet->{$proj}}){
        foreach my $tag (keys %{$sampleSheet->{$proj}->{$lid}}){
            my $info = $sampleSheet->{$proj}->{$lid}->{$tag};
            my $laneId = ("0" x (3 - length($lid))) . $lid;
            my $fastQFilesFound = 0;
            my $indexOrSampleCounter = $info->{SampleNumber};
            foreach my $read (keys %{$files{$laneId}{$info->{SampleID}}{$indexOrSampleCounter}}){
                foreach my $pctLane (keys %{$files{$laneId}{$info->{SampleID}}{$indexOrSampleCounter}{$read}}){
                    my $stat = Molmed::Sisyphus::QStat->new(OFFSET=>$OFFSET, DEBUG=>$debug);
                    my $filehandle;
                    $samples->{$info->{SampleName}}->{$lid}->{$tag} = 1;
                    if($files{$laneId}{$info->{SampleID}}{$indexOrSampleCounter}{$read}{$pctLane} =~ /fastq.gz$/) {
                        open($filehandle, "zcat $files{$laneId}{$info->{SampleID}}{$indexOrSampleCounter}{$read}{$pctLane} |") or die "Failed to open $files{$laneId}{$info->{SampleID}}{$indexOrSampleCounter}{$read}{$pctLane}: $!";
                    } elsif($files{$laneId}{$info->{SampleID}}{$indexOrSampleCounter}{$read}{$pctLane} =~ /fastq$/) {
                        open($filehandle,'-|', "grep fastq.gz $files{$laneId}{$info->{SampleID}}{$indexOrSampleCounter}{$read}{$pctLane}") or die "Failed to open $files{$laneId}{$info->{SampleID}}{$indexOrSampleCounter}{$read}{$pctLane}: $!";
                    } 

                    $fastQFilesFound = 1;
                    my $seq = "";
                    my $qual = "";
                    my $counter = 0;

                    FASTQ: while(<$filehandle>) {
                        $seq = <$filehandle>;
                        chomp($seq);
                        <$filehandle>;
                        $qual = <$filehandle>;
                        chomp($qual);
                        if(rand() < $stat->{SAMPLING_DENSITY}) {
                            $counter = $counter + 1;
                            $stat->addQValuePerBaseAndPosition($seq,$qual);
                            last FASTQ if($counter > $stat->{SAMPLING_COUNTER});
                        }
                    }
                    close($filehandle);
                    $baseQC->{$info->{SampleID}}->{$lid}->{$read}->{$indexOrSampleCounter} = $stat->calculateQValuePerBase();
                    if(!defined($laneQC->{$lid}->{$read})) {
                        $laneQC->{$lid}->{$read} = Molmed::Sisyphus::QStat->new(OFFSET=>$OFFSET, DEBUG=>$debug);	
                    }
                    $laneQC->{$lid}->{$read} = $laneQC->{$lid}->{$read}->add($stat);
                }
            }
            if($fastQFilesFound == 0) {
                print "Couldn't find fastq files(s) for lane $laneId, $info->{SampleID}\n";
                unless($ignoreFastqs){ 
                    exit 1;
                }
            }
        }
    }
}

my %laneFrac;
my %laneUnknown;
my $reads;

foreach my $sample (keys %{$RtaSampleStats}){
    foreach my $lane (keys %{$RtaSampleStats->{$sample}}){
        foreach my $barcode (keys %{$RtaSampleStats->{$sample}->{$lane}->{1}}){
            if($barcode eq 'Undetermined' || ($barcode eq 'unknown' && $sample eq 'Undetermined')){
                $laneUnknown{$lane} = sprintf('%.1f', $RtaSampleStats->{$sample}->{$lane}->{1}->{$barcode}->{PctLane});
            }
            else{
                if(defined($samples->{$sample}->{$lane}->{$barcode})) {
                    push @{$laneFrac{$lane}}, sprintf(' %.1f:%s', $RtaSampleStats->{$sample}->{$lane}->{1}->{$barcode}->{PctLane}, $sample);
            }
            }
        }
    }
}

open(my $repFh, '>', "$rfPath/quickReport.txt");
print $repFh join("\t", "Lane", "Read", "ReadsPF (M)", "Yield Q30 (G)", "ErrRate", "Excluded", "Q per base (A/C/G/T)", "Sample Fractions", "Unidentified"), "\n";
foreach my $lane (sort {$a<=>$b} keys %{$RtaLaneStats}){
    foreach my $read (sort {$a<=>$b} keys %{$RtaLaneStats->{$lane}}){
        print $repFh join("\t", $lane, $read,
            sprintf('%.0f', (defined($RtaLaneStats->{$lane}->{$read}->{PF}) ? $RtaLaneStats->{$lane}->{$read}->{PF} : 0)/1e6),
            sprintf('%.1f', (defined($RtaLaneStats->{$lane}->{$read}->{YieldQ30}) ? $RtaLaneStats->{$lane}->{$read}->{YieldQ30} : 0)/1e9),
            defined($RtaLaneStats->{$lane}->{$read}->{ErrRate}) ? $RtaLaneStats->{$lane}->{$read}->{ErrRate} : '-',
            $RtaLaneStats->{$lane}->{$read}->{ExcludedTiles});

        my @fractionSorted = defined($laneFrac{$lane}) ? sort({sortLaneFrac($a,$b)} @{$laneFrac{$lane}}) : ('-');

        if(defined($laneQC->{$lane}) && defined($laneQC->{$lane}->{$read})) {
            my $result = $laneQC->{$lane}->{$read}->calculateQValuePerBase();

            print $repFh "\t";
            printf $repFh sprintf('%.1f±%.1f,', $result->{"a"}->{MEAN},$result->{'a'}->{STDV});
            printf $repFh sprintf('%.1f±%.1f,', $result->{"c"}->{MEAN},$result->{'c'}->{STDV});
            printf $repFh sprintf('%.1f±%.1f,', $result->{"g"}->{MEAN},$result->{"g"}->{STDV});
            printf $repFh sprintf('%.1f±%.1f', $result->{"t"}->{MEAN},$result->{'t'}->{STDV});
        }else {
            print $repFh "\t-,-,-,-";
        }
        print $repFh "\t", join(',', @fractionSorted);

        print $repFh "\t", (defined($laneUnknown{$lane}) ? $laneUnknown{$lane} : 0);
 
        print $repFh "\n";
    }
}
close($repFh);

if(defined $mail && $mail =~ m/\w\@\w/){
    open(my $repFh, '<', "$rfPath/quickReport.txt");
    my $msg = '<html><body><table>' . "\n";
    my $i=0;
    while(<$repFh>){
        $i++;
        s:\t:</td><td>:g;
        s:,:<br />:g;
        if($i==1){
            s/td/th/g;
            $msg .= '<tr><th>' . $_ . '</tr>';
        }elsif($i%2 > 0){
            $msg .= '<tr><td>' . $_ . '</tr>';
        }else{
        $msg .= '<tr bgcolor="#dddddd"><td>' . $_ . '</tr>';
        }
    }
    $msg .= '</table>'. "\n";
    $msg .= '</body></html>';

    require Net::SMTP;
    #Create a new object with 'new'.
    my $smtp = Net::SMTP->new("smtp.uu.se");
    #Send the MAIL command to the server.
    $smtp->mail($sender);
    #Send the server the 'Mail To' address.
    $smtp->to($mail);
    #Start the message.
    $smtp->data();
    #Send the message.
    $smtp->datasend("From: $sender\n");
    $smtp->datasend("To: $mail\n");
    $smtp->datasend("Subject: [Sisyphus] [Quick Report]" . basename($rfPath) . "\n");
    $smtp->datasend("MIME-Version: 1.0\n");
    $smtp->datasend("Content-Type: text/html; charset=us-ascii\n");
    $smtp->datasend("\n");
    $smtp->datasend("$msg\n\n");
    #End the message.
    $smtp->dataend();
    #Close the connection to your server.
    $smtp->quit();
}


sub sortLaneFrac{
    my @a = split ':', $_[0];
    my @b = split ':', $_[1];
    return($a[0]<=>$b[0]);
}