Merge branch 'release/v3.0.0'

README.md

@@ -3,7 +3,7 @@ AscatNGS
 AscatNGS contains the Cancer Genome Projects workflow implementation of the ASCAT copy number
 algorithm for paired end sequencing.
 
-For details of the underlying algorithm please see the [ASCAT](http://heim.ifi.uio.no/bioinf/Projects/ASCAT/) site.
+For details of the underlying algorithm please see the [ASCAT](https://www.crick.ac.uk/peter-van-loo/software/ASCAT) site.  The GitHub repoitory can be found [here](https://github.com/Crick-CancerGenomics/ascat).
 
 ---
 

perl/MANIFEST

@@ -24,18 +24,18 @@ docs/pod_html/ascatFailedCnCsv.html
 docs/pod_html/index.html
 docs/pod_html/Sanger/CGP/Ascat.html
 docs/pod_html/utilities/ascatSnpPanelGenerator.html
-docs/reports_html/blib-lib-Sanger-CGP-Ascat-Implement-pm--branch.html
-docs/reports_html/blib-lib-Sanger-CGP-Ascat-Implement-pm--subroutine.html
-docs/reports_html/blib-lib-Sanger-CGP-Ascat-Implement-pm.html
-docs/reports_html/blib-lib-Sanger-CGP-Ascat-pm--subroutine.html
-docs/reports_html/blib-lib-Sanger-CGP-Ascat-pm.html
 docs/reports_html/common.js
 docs/reports_html/cover.14
 docs/reports_html/cover.css
 docs/reports_html/coverage.html
 docs/reports_html/css.js
 docs/reports_html/digests
 docs/reports_html/index.html
+docs/reports_html/lib-Sanger-CGP-Ascat-Implement-pm--branch.html
+docs/reports_html/lib-Sanger-CGP-Ascat-Implement-pm--subroutine.html
+docs/reports_html/lib-Sanger-CGP-Ascat-Implement-pm.html
+docs/reports_html/lib-Sanger-CGP-Ascat-pm--subroutine.html
+docs/reports_html/lib-Sanger-CGP-Ascat-pm.html
 docs/reports_html/standardista-table-sorting.js
 docs/reports_text/coverage.txt
 lib/Sanger/CGP/Ascat.pm

perl/bin/ascat.pl

@@ -84,6 +84,7 @@ sub setup {
   my %opts;
   pod2usage(-msg  => "\nERROR: Option must be defined.\n", -verbose => 1,  -output => \*STDERR) if(scalar @ARGV == 0);
   $opts{'cmd'} = join " ", $0, @ARGV;
+  warn "Executing: $opts{cmd}\n";
   GetOptions( 'h|help' => \$opts{'h'},
               'm|man' => \$opts{'m'},
               'v|version' => \$opts{'v'},
@@ -95,8 +96,8 @@ sub setup {
               'c|cpus=i' => \$opts{'threads'},
               'x|limit=i' => \$opts{'limit'},
               'q|minbasequal=i' => \$opts{'minbasequal'},
-              's|snp_loci=s' => \$opts{'snp_loci'},
-              'sp|snp_pos=s' => \$opts{'snp_pos'},
+              's|snp_loci=s' => \$opts{'snp_loci'}, # still here so it doesn't break if defined
+              'sp|snp_pos=s' => \$opts{'snp_pos'}, # still here so it doesn't break if defined
               'sg|snp_gc=s' => \$opts{'snp_gc'},
               'g|gender=s' => \$opts{'gender'},
               'gc|genderChr=s' => \$opts{'genderChr'},
@@ -112,8 +113,8 @@ sub setup {
               'nc|noclean' => \$opts{'noclean'},
   ) or pod2usage(2);
 
-  pod2usage(-message => Sanger::CGP::Ascat::license, -verbose => 1) if(defined $opts{'h'});
-  pod2usage(-verbose => 2) if(defined $opts{'m'});
+  pod2usage(-verbose => 1, -exitval => 0) if(defined $opts{'h'});
+  pod2usage(-verbose => 2, -exitval => 0) if(defined $opts{'m'});
 
   if($opts{'v'}){
     print Sanger::CGP::Ascat->VERSION."\n";
@@ -129,22 +130,33 @@ sub setup {
   for(keys %opts) { $defined++ if(defined $opts{$_}); }
   pod2usage(-msg  => "\nERROR: Options must be defined.\n", -verbose => 1,  -output => \*STDERR) unless($defined);
 
+  for my $item (qw(tumour normal snp_gc reference outdir protocol gender)) {
+    pod2usage(-msg  => "\nERROR: Option '-$item' must be defined.\n", -verbose => 1,  -output => \*STDERR) unless(defined $opts{$item});
+  }
+
   if((defined($opts{'purity'}) && !defined($opts{'ploidy'})) || (!defined($opts{'purity'}) && defined($opts{'ploidy'}))){
     pod2usage(-msg  => "\nERROR: If one of purity/ploidy are defined, both should be defined.\n", -verbose => 1,  -output => \*STDERR);
   }
 
-  for my $item(qw(tumour normal snp_loci snp_pos snp_gc reference outdir locus)) {
+  for my $item(qw(tumour normal snp_gc reference outdir locus)) {
     $opts{$item} = File::Spec->rel2abs( $opts{$item} ) if(defined $opts{$item});
   }
 
   PCAP::Cli::file_for_reading('tumour', $opts{'tumour'});
   PCAP::Cli::file_for_reading('normal', $opts{'normal'});
-  PCAP::Cli::file_for_reading('snp_loci', $opts{'snp_loci'});
-  PCAP::Cli::file_for_reading('snp_pos', $opts{'snp_pos'});
   PCAP::Cli::file_for_reading('snp_gc', $opts{'snp_gc'});
   PCAP::Cli::file_for_reading('reference', $opts{'reference'});
   PCAP::Cli::out_dir_check('outdir', $opts{'outdir'});
 
+  if(defined $opts{'snp_loci'}) {
+    warn qq{NOTE: '-snp_loci' file is no longer required\n};
+    delete $opts{'snp_loci'};
+  }
+  if(defined $opts{'snp_pos'}) {
+    warn qq{NOTE: '-snp_pos' file is no longer required\n};
+    delete $opts{'snp_pos'};
+  }
+
   my $final_logs = File::Spec->catdir($opts{'outdir'}, 'logs');
   if(-e $final_logs) {
     warn "NOTE: Presence of '$final_logs' directory suggests successful complete analysis, please delete to rerun\n";
@@ -179,7 +191,7 @@ sub setup {
     }
   }
   elsif(exists $opts{'index'}) {
-    die "ERROR: -index cannot be defined without -process\n";
+    pod2usage(-msg  => "\nERROR: -index cannot be defined without -process.\n", -verbose => 1,  -output => \*STDERR);
   }
 
   $opts{'threads'} = 1 unless(defined $opts{'threads'});
@@ -194,20 +206,20 @@ sub setup {
   $opts{'tmp'} = $tmpdir;
 
   if(defined $opts{'gender'}){
-    pod2usage(-message => 'unknown gender value: '.$opts{'gender'}, -verbose => 1) unless(first {$_ eq $opts{'gender'}} @VALID_GENDERS);
+    pod2usage(-message => "\nERROR: Unknown gender value: $opts{gender}\n", -verbose => 1) unless(first {$_ eq $opts{'gender'}} @VALID_GENDERS);
     if($opts{'gender'} eq 'L') {
       my ($is_male, $gender_chr) = Sanger::CGP::Ascat::Implement::determine_gender(\%opts);
       $opts{'genderChr'} = $gender_chr;
       $opts{'genderIsMale'} = $gender_chr;
       $opts{'gender'} = $is_male eq 'N' ? 'XX' : 'XY';
     }
     else {
-      pod2usage(-message => 'genderChr must be set when gender is XX/XY', -verbose => 1) if(!defined $opts{'genderChr'});
-      pod2usage(-message => 'gender must be XX, XY or L', -verbose => 1)if($opts{'gender'} !~ m/^X[XY]$/);
+      pod2usage(-message => "\nERROR: genderChr must be set when gender is XX/XY\n", -verbose => 1) if(!defined $opts{'genderChr'});
+      pod2usage(-message => "\nERROR: gender must be XX, XY or L\n", -verbose => 1)if($opts{'gender'} !~ m/^X[XY]$/);
       $opts{'genderIsMale'} = $opts{'gender'} eq 'XX' ? 'N' : 'Y';
     }
   } else {
-    pod2usage(-message => 'gender not set', -verbose => 1);
+    pod2usage(-message => "\nERROR: gender not set\n", -verbose => 1);
   }
 
   return \%opts;
@@ -234,9 +246,8 @@ =head1 SYNOPSIS
     -tumour       -t    Tumour BAM/CRAM file
     -normal       -n    Normal BAM/CRAM file
     -reference    -r    Reference fasta
-    -snp_loci     -s    Snp locus file
-    -snp_pos      -sp   Snp position file
     -snp_gc       -sg   Snp GC correction file
+    -protocol     -pr   Sequencing protocol (e.g. WGS, WXS)
     -gender       -g    Sample gender (XX, XY, L)
                           For XX/XY see '-gc'
                           When 'L' see '-l'
@@ -251,7 +262,6 @@ =head1 SYNOPSIS
     -genderChr    -gc   Specify the 'Male' sex chromosome: Y,chrY...
     -species      -rs   Reference species [BAM HEADER]
     -assembly     -ra   Reference assembly [BAM HEADER]
-    -protocol     -pr   Sequencing protocol (e.g. WGS, WXS)
     -platform     -pl   Seqeuncing platform [BAM HEADER]
     -minbasequal  -q    Minimum base quality required before allele is used. [20]
     -cpus         -c    Number of cores to use. [1]

perl/lib/Sanger/CGP/Ascat.pm

@@ -26,7 +26,7 @@ use strict;
 use Const::Fast qw(const);
 use base 'Exporter';
 
-our $VERSION = '2.1.3';
+our $VERSION = '3.0.0';
 our @EXPORT = qw($VERSION);
 
 const my $LICENSE =>

perl/lib/Sanger/CGP/Ascat/Implement.pm

@@ -35,6 +35,7 @@ use File::Temp qw(tempfile);
 use File::Copy qw(copy move);
 use Capture::Tiny qw(capture);
 use FindBin qw($Bin);
+use List::Util qw(first);
 
 use File::ShareDir qw(module_dir);
 
@@ -53,7 +54,7 @@ const my $ALLELE_COUNT_PARA => ' -b %s -o %s -l %s -c %s -r %s ';
 
 const my $GREP_ALLELE_COUNTS => q{grep -v '^#' %s >> %s};
 
-const my @ASCAT_RESULT_FILES => qw( %s.aberrationreliability.png
+const my @ASCAT_RESULT_FILES => qw(
                                     %s.ASCATprofile.png
                                     %s.ASPCF.png
                                     %s.germline.png
@@ -100,19 +101,25 @@ sub allele_count {
 
     my $ac_out = File::Spec->catdir($tmp, 'allele_count');
     make_path($ac_out) unless(-e $ac_out);
+    my $loci_files = File::Spec->catdir($tmp, 'allele_count', 'loci_files');
+    make_path($loci_files) unless(-e $loci_files);
 
     my $chr = $seqs[$seq_idx-1];
     my $sname = sanitised_sample_from_bam($options->{$samp_type});
     my $alleleCountOut = File::Spec->catfile($ac_out,sprintf '%s.%s.allct', $sname, $chr);
 
+    # first we need a loci file, generate from the gc file:
+    my $loci_file = "$loci_files/$sname.$chr";
+    my $gc_to_loci = qq{cut -f 2,3 $options->{snp_gc} | grep -vP '^Chr\\tPosition' | grep -P '^$chr\t' > $loci_file};
+
     my $allc_exe = _which('alleleCounter');
     my $allc_lib = dirname($allc_exe);
 
     my $command = $allc_exe;
-    $command .= sprintf $ALLELE_COUNT_PARA, $options->{$samp_type}, $alleleCountOut, $options->{'snp_loci'}, $chr, $options->{'reference'};
+    $command .= sprintf $ALLELE_COUNT_PARA, $options->{$samp_type}, $alleleCountOut, $loci_file, $chr, $options->{'reference'};
     $command .= '-m '.$options->{'minbasequal'} if exists $options->{'minbasequal'};
 
-    PCAP::Threaded::external_process_handler(File::Spec->catdir($tmp, 'logs'), $command, $index);
+    PCAP::Threaded::external_process_handler(File::Spec->catdir($tmp, 'logs'), [$gc_to_loci, $command], $index);
 
     PCAP::Threaded::touch_success(File::Spec->catdir($tmp, 'progress'), $index);
   }
@@ -149,7 +156,12 @@ sub ascat {
     PCAP::Threaded::touch_success(File::Spec->catdir($tmp, 'progress'), 'merge_counts_wt', 0);
   }
 
-  my $core_chrs = snpLociChrs($options);
+  my $snp_pos = File::Spec->catfile($tmp,'SnpPositions.tsv');
+  my $gc_to_snppos = qq{cut -f 1-3 $options->{snp_gc} > $snp_pos};
+
+  my @chr_set = snpLociChrs($options);
+  my $core_chrs = @chr_set;
+  $core_chrs++ unless(first {$_ eq 'Y'  || $_ eq 'chrY'} @chr_set); # incase Y is not in the SNP set
 
   my $command = "cd $ascat_out; "._which('Rscript');
 
@@ -161,7 +173,7 @@ sub ascat {
 
   $command .= " $ascat_exe";
   $command .= " $ascat_path";
-  $command .= ' '.$options->{'snp_pos'};
+  $command .= ' '.$snp_pos;
   $command .= ' '.$options->{'snp_gc'};
   $command .= ' '.$options->{'tumour_name'};
   $command .= ' '.$tumcountfile;
@@ -176,7 +188,7 @@ sub ascat {
     $command .= ' '.$options->{'ploidy'};
   }
 
-  PCAP::Threaded::external_process_handler(File::Spec->catdir($tmp, 'logs'), $command, 0);
+  PCAP::Threaded::external_process_handler(File::Spec->catdir($tmp, 'logs'), [$gc_to_snppos, $command], 0);
 
   PCAP::Threaded::touch_success(File::Spec->catdir($tmp, 'progress'), 0);
 
@@ -337,13 +349,20 @@ sub determine_gender {
   return ($norm_gender, $gender_chr);
 }
 
+=head
+
+Updated to run from SnpGcPositions as this file is the only one technically needed
+
+=cut
+
 sub snpLociChrs {
   my $options = shift;
   my %chr_set;
   my @chrs;
-  open my $IN, '<', $options->{'snp_loci'};
+  open my $IN, '<', $options->{'snp_gc'};
   while(my $line = <$IN>) {
-    my ($chr, undef) = split /\t/, $line;
+    next if($line =~ m/Chr\tPosition/);
+    my (undef, $chr, undef) = split /\t/, $line;
     unless(exists $chr_set{$chr}) {
       $chr_set{$chr} = 1;
       push @chrs, $chr;