bamqc

#!/usr/bin/perl
use warnings;
use strict;
use FindBin qw($RealBin);
use Getopt::Long;

#####################################################################################
# Copyright Copyright 2010-16 Simon Andrews                                         #
#                                                                                   #
#    This file is part of BamQC.                                                    #
#                                                                                   #
#    BamQC is free software; you can redistribute it and/or modify                  #
#    it under the terms of the GNU General Public License as published by           #
#    the Free Software Foundation; either version 3 of the License, or              #
#    (at your option) any later version.                                            #
#                                                                                   #
#    BamQC is distributed in the hope that it will be useful,                       #
#    but WITHOUT ANY WARRANTY; without even the implied warranty of                 #
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                  #
#    GNU General Public License for more details.                                   #
#                                                                                   #
#    You should have received a copy of the GNU General Public License              #
#    along with BamQC; if not, write to the Free Software                           #
#    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA     #
#####################################################################################


# Check to see if they've mistakenly downloaded the source distribution
# since several people have made this mistake

if (-e "$RealBin/uk/ac/babraham/BamQC/BamQCApplication.java") {
	die "This is the source distribution of BamQC.  You need to get the compiled version if you want to run the program\n";
}

my $delimiter = ':';

if ($^O =~ /Win/) {
	$delimiter = ';';
}

if ($ENV{CLASSPATH}) {
	$ENV{CLASSPATH} .= "$delimiter$RealBin$delimiter$RealBin/sam-1.103.jar$delimiter$RealBin/jbzip2-0.9.1.jar$delimiter$RealBin/hamcrest-core-1.3.jar$delimiter$RealBin/junit-4.11.jar$delimiter$RealBin/log4j-1.2.16.jar$delimiter$RealBin/commons-lang3-3.4.jar$delimiter$RealBin/commons-math3-3.5.jar";
}
else {
	$ENV{CLASSPATH} = "$RealBin$delimiter$RealBin/sam-1.103.jar$delimiter$RealBin/jbzip2-0.9.1.jar$delimiter$RealBin/hamcrest-core-1.3.jar$delimiter$RealBin/junit-4.11.jar$delimiter$RealBin/log4j-1.2.16.jar$delimiter$RealBin/commons-lang3-3.4.jar$delimiter$RealBin/commons-math3-3.5.jar";
}

# "$delimiter$RealBin$delimiter$RealBin/sam-1.103.jar$delimiter$RealBin/jbzip2-0.9.1.jar$delimiter$RealBin/hamcrest-core-1.3.jar$delimiter$RealBin/junit-4.11.jar$delimiter$RealBin/log4j-1.2.16.jar$delimiter$RealBin/commons-lang3-3.4.jar$delimiter$RealBin/commons-math3-3.5.jar



my @java_args;
my @files;

# We now need to scan the command line for switches which we're going
# to pass on to the main java program.



my $version;
my $help;
my $outdir;
my $unzip;
my $gff;
my $genome;
my $species;
my $assembly;
my $genome_pattern = 'unassigned';
my $saved;
my $limits;
my $threads;
my $quiet;
my $nogroup;
my $expgroup;
my $temp_directory;
my $java_bin = 'java';

my $result = GetOptions('v|version' => \$version,
			'h|help' => \$help,
			'b|available:s' => \$genome_pattern,
			'e|saved' => \$saved,
			'q|quiet' => \$quiet,
			'nogroup' => \$nogroup,
			'expgroup' => \$expgroup,
			'o|outdir=s' => \$outdir,
			'extract!' => \$unzip,
			't|threads=i' => \$threads,
			'f|gff=s' => \$gff,
			'g|genome=s' => \$genome,
			's|species=s' => \$species,
			'a|assembly=s' => \$assembly,
			'l|limits=s' => \$limits,
			'd|dir=s' => \$temp_directory,
			'j|java=s' => \$java_bin,
);

# Check the simple stuff first

if ($help) {
	# Just print the help and exit
	print while(<DATA>);
	exit;
}

if ($version) {
	push @java_args ,"-Dbamqc.show_version=true";
}

if ($genome_pattern ne 'unassigned') {
	push @java_args ,"-Dbamqc.available_genomes=true";
	# For some reason, if a user uses the filter '*', this is converted 
	# to 'bamqc' by Perl. The following test is used to reverse this back.
	if($genome_pattern eq 'bamqc') {
		push @java_args ,"-Dbamqc.genome_pattern=*";
	} else {
		push @java_args ,"-Dbamqc.genome_pattern=$genome_pattern";
	}
}

if ($saved) {
	push @java_args ,"-Dbamqc.saved_genomes=true";
}


# Now parse any additional options
if ($outdir) {
	unless(-e $outdir and -d $outdir) {
		die "Specified output directory '$outdir' does not exist\n";
	}
	
	push @java_args ,"-Dbamqc.output_dir=$outdir";
}

our $contaminant;

if ($contaminant)  {
	unless (-e $contaminant and -r $contaminant) {
		die "Contaminant file '$contaminant' did not exist, or could not be read\n";
	}
	push @java_args ,"-Dbamqc.contaminant_file=$contaminant";
}

our $adapter;

if ($adapter)  {
	unless (-e $adapter and -r $adapter) {
		die "Adapters file '$adapter' did not exist, or could not be read\n";
	}
	push @java_args ,"-Dbamqc.adapter_file=$adapter";
}

if ($gff)  {
	unless (-e $gff and -r $gff) {
		die "GFF file '$gff' did not exist, or could not be read\n";
	}
	push @java_args ,"-Dbamqc.gff_file=$gff";
}

if ($genome)  {
# We will download this later if this does not exist. Let's not check anything here.
#	unless (-e $genome and -r $genome) {
#		die "Genome '$genome' did not exist, or could not be read\n";
#	}
	push @java_args ,"-Dbamqc.genome=$genome";
}

if ($genome)  {
	push @java_args ,"-Dbamqc.genome=$genome";
}

if ($species and $assembly)  {
	push @java_args ,"-Dbamqc.species=$species";
	push @java_args ,"-Dbamqc.assembly=$assembly";
} elsif ($species or $assembly)  {
	die "If used, the parameters -s (--species) and -a (--assembly) must be specified together\n";  
}

if ($limits)  {
	unless (-e $limits and -r $limits) {
		die "Limits file '$limits' did not exist, or could not be read\n";
	}
	push @java_args ,"-Dbamqc.limits_file=$limits";
}

if ($temp_directory) {
	unless (-e $temp_directory and -d $temp_directory and -w $temp_directory) {
		die "Temp directory '$temp_directory' doesn't exist, or can't be written to\n";
	}
	push @java_args, "-Djava.io.tmpdir=$temp_directory";

}

if ($threads) {
	if ($threads < 1) {
		die "Number of threads must be a positive integer";
	}
	
	push @java_args ,"-Dbamqc.threads=$threads";
	my $memory = 250 * $threads;
	unshift @java_args,"-Xmx${memory}m";
}
else {
	unshift @java_args,'-Xmx250m';
}

if ($quiet) {
	push @java_args ,"-Dbamqc.quiet=true";	
}

if ($nogroup) {
	push @java_args ,"-Dbamqc.nogroup=true";	
}

if (defined $unzip) {
	
	if ($unzip) {
		$unzip = 'true';
	}
	else {
		$unzip = 'false';
	}
	
	push @java_args,"-Dbamqc.unzip=$unzip";	
}


if ($java_bin ne 'java') {
#	$java_bin =~ s/\\/\//g;
	
	unless (-e $java_bin) {
		die "Couldn't find java interpreter at '$java_bin'";
	}
	
#	if ($java_bin =~ / /) {
#		$java_bin = "\"$java_bin\"";
#	}	
}


# We've found that on systems with large numbers of CPUs, and especially on systems
# where IO contention can occur, that the program can consume large amounts of CPU
# during parallel garbage collection.  Since interactive performance isn't really
# important when running non-interactively we'll limit the GC to using only a single
# thread if running with a list of filenames

if (@ARGV) {
	push @java_args,"-XX:ParallelGCThreads=1";
}


foreach (@ARGV) {
  if (/^\-D/) {
    push @java_args,$_;
  }
  else {
    push @files,$_;
  }
}

# This is set internally as well, but on some JREs it doesn't
# pick up the internally set value properly, so we'll set it
# outside as well which should work.
if (@files or $version or $help or $genome_pattern ne 'unassigned' or $saved) {
	push @java_args, "-Djava.awt.headless=true";
} else {
    # Antialiasing by default on text  
	push @java_args, "-Dawt.useSystemAAFontSettings=gasp";
}

#warn "Running 'exec $java_bin,@java_args, \"uk.ac.babraham.BamQC.BamQCApplication\", @files, CLASSPATH is $ENV{CLASSPATH}";


if ($java_bin ne 'java') {
  if ($genome_pattern ne 'unassigned' or $saved) {
      system $java_bin, @java_args, "uk.ac.babraham.BamQC.Utilities.BamQCListGenomes";
  } else {
      system $java_bin, @java_args, "uk.ac.babraham.BamQC.BamQCApplication", @files;
  }
} else {
  if ($genome_pattern ne 'unassigned' or $saved) {
      exec $java_bin, @java_args, "uk.ac.babraham.BamQC.Utilities.BamQCListGenomes";
  } else {
      exec $java_bin, @java_args, "uk.ac.babraham.BamQC.BamQCApplication", @files;
  }
}






__DATA__

        BamQC - A high throughput mapped sequence QC analysis tool

SYNOPSIS

    bamqc bamfile1 .. bamfileN | <folder containing SAM/BAM mapped files> 

    bamqc [-o output dir] [--(no)extract] [-f file.gtf] 
           bamfile1 .. bamfileN | <folder containing SAM/BAM mapped files> 
           
    bamqc [-o output dir] [--(no)extract] [-g genome_dir] 
           bamfile1 .. bamfileN | <folder containing SAM/BAM mapped files>          

    bamqc [-o output dir] [--(no)extract] [-p species -e assembly] 
           bamfile1 .. bamfileN | <folder containing SAM/BAM mapped files>                 

DESCRIPTION

    BamQC reads a set of mapped BAM files and produces from each one a quality
    control report consisting of a number of different modules, each one of 
    which will help to identify a different potential type of problem in your
    data.
    
    If no files to process are specified on the command line then the program
    will start as an interactive graphical application.  If files are provided
    on the command line then the program will run with no user interaction
    required.  In this mode it is suitable for inclusion into a standardised
    analysis pipeline.
    
    The options for the program as as follows:

    -f --gff        Use a specified annotation file as annotation set
    
    -g --genome     The directory containing species/assembly to use. If the 
                    couple species assembly does not exist, BamQC will try to 
                    download it.
    
    -s --species    The genome species to use. If the couple species assembly
                    does not exist, BamQC will try to download it.
    
    -a --assembly   The genome assembly associated to the species to use. 
                    If the couple species assembly does not exist, BamQC 
                    will try to download it.  
     
    -b --available [Pattern] List the genomes available on the Babraham Server
    				which are filtered using the string Pattern if this is provided. 
    				This string is considered as a substring in the pattern matching. 
    				For instance, the pattern=m?s will retrieve mus, mis, musculus, 
    				humus. Therefore its behaviour is equivalent to *m?s*.
    
    -e --saved      List the saved genomes
    
    -h --help       Print this help file and exit
    
    -v --version    Print the version of the program and exit
    
    -o --outdir     Create all output files in the specified output directory.
                    Please note that this directory must exist as the program
                    will not create it.  If this option is not set then the 
                    output file for each sequence file is created in the same
                    directory as the sequence file which was processed.
                                       
    --extract       If set then the zipped output file will be uncompressed in
                    the same directory after it has been created.  By default
                    this option will be set if bamqc is run in non-interactive
                    mode.
                    
    -j --java       Provides the full path to the java binary you want to use to
                    launch bamqc. If not supplied then java is assumed to be in
                    your path.
                   
    --noextract     Do not uncompress the output file after creating it.  You
                    should set this option if you do not wish to uncompress
                    the output when running in non-interactive mode.
                    
    --nogroup       Disable grouping of bases for reads >50bp. All reports will
                    show data for every base in the read.  WARNING: Using this
                    option will cause bamqc to crash and burn if you use it on
                    really long reads, and your plots may end up a ridiculous size.
                    You have been warned!
                                        
    -t --threads    Specifies the number of files which can be processed
                    simultaneously.  Each thread will be allocated 250MB of
                    memory so you shouldn't run more threads than your
                    available memory will cope with, and not more than
                    6 threads on a 32 bit machine
                                      
    -l --limits     Specifies a non-default file which contains a set of criteria
                    which will be used to determine the warn/error limits for the
                    various modules.  This file can also be used to selectively 
                    remove some modules from the output all together.  The format
                    needs to mirror the default limits.txt file found in the
                    Configuration folder.
                    
   -q --quiet       Supress all progress messages on stdout and only report errors.
   
   -d --dir         Selects a directory to be used for temporary files written when
                    generating report images. Defaults to system temp directory if
                    not specified.
                    
BUGS

    Any bugs in bamqc should be reported either to simon.andrews@babraham.ac.uk
    or in www.bioinformatics.babraham.ac.uk/bugzilla/