pughlab_rnaseq_pipeline.pl

#!/usr/bin/env perl
### pughlab_rnaseq_pipeline.pl #####################################################################
use AutoLoader 'AUTOLOAD';
use strict;
use warnings;
use Carp;
use POSIX qw(strftime);
use Getopt::Std;
use Getopt::Long;
use File::Basename;
use File::Path qw(make_path);
use YAML qw(LoadFile);
use List::Util 'any';
use Data::Dumper;

my $cwd = dirname(__FILE__);
require "$cwd/scripts/utilities.pl";

####################################################################################################
# version       author	  	comment
# 1.0		sprokopec       script to run PughLab RNASeq pipeline

### USAGE ##########################################################################################
# pughlab_rnaseq_pipeline.pl -c tool_config.yaml -d data.yaml
#
# where:
#	- tool_config.yaml contains tool versions and parameters, output directory, reference
#	information, etc.
#	- data_config.yaml contains sample information (YAML file containing paths to FASTQ files,
#	generated by create_fastq_yaml.pl)

### SUBROUTINES ####################################################################################

### MAIN ###########################################################################################
sub main {
	my %args = (
		tool_config	=> undef,
		data_config	=> undef,
		report		=> undef,
		cleanup		=> undef,
		cluster		=> undef,
		dry_run		=> undef,
		@_
		);

	my $tool_config = $args{tool_config};
	my $data_config = $args{data_config};

	### PREAMBLE ######################################################################################

	# load tool config
	my $tool_data = LoadFile($tool_config);
	my $date = strftime "%F", localtime;
	my $timestamp = strftime "%F_%H-%M-%S", localtime;

	# check for and/or create output directories
	my $output_directory = $tool_data->{output_dir};
	$output_directory =~ s/\/$//;
	my $log_directory = join('/', $output_directory, 'logs', 'run_RNA_pipeline_' . $timestamp);

	unless(-e $output_directory) { make_path($output_directory); }
	unless(-e $log_directory) { make_path($log_directory); }

	# start logging
	my $log_file = join('/', $log_directory, 'run_RNASeq_pipeline.log');
	open (my $log, '>', $log_file) or die "Could not open $log_file for writing.";

	print $log "---\n";
	print $log "Running PughLab RNA-Seq pipeline.\n";
	print $log "\n  Tool config used: $tool_config";
	print $log "\n  Sample config used: $data_config";
	print $log "\n  Output directory: $output_directory";
	print $log "\n---\n\n";

	my $perl = 'perl/' . $tool_data->{perl_version};

	# get optional HPC group
	my $hpc_group = defined($tool_data->{hpc_group}) ? "-A $tool_data->{hpc_group}" : undef;

	### MAIN ###########################################################################################

	my ($fastqc_run_id, $star_run_id, $gatk_run_id, $vc_run_id, $rsem_run_id);
	my ($fc_run_id, $starfus_run_id, $arriba_run_id, $mavis_run_id);
	my ($run_script, $report_run_id);

	my @job_ids;

	# prepare directory structure
	my $fastqc_directory = join('/', $output_directory, 'fastqc');
	my $fc_directory = join('/', $output_directory, 'FusionCatcher');
	my $star_directory = join('/', $output_directory, 'STAR');
	my $starfus_directory = join('/', $output_directory, 'STAR-Fusion');
	my $arriba_directory = join('/', $output_directory, 'Arriba');
	my $rsem_directory = join('/', $output_directory, 'RSEM');
	my $gatk_directory = join('/', $output_directory, 'GATK');
	my $vc_directory = join('/', $output_directory, 'HaplotypeCaller');
	my $mavis_directory = join('/', $output_directory, 'Mavis');

	# check which tools have been requested
	my %tool_set = (
		'fastqc' => defined($tool_data->{fastqc}->{run}) ? $tool_data->{fastqc}->{run} : 'N',
		'star'	=> defined($tool_data->{star}->{run}) ? $tool_data->{star}->{run} : 'N',
		'rsem'	=> defined($tool_data->{rsem}->{run}) ? $tool_data->{rsem}->{run} : 'N',
		'gatk'	=> defined($tool_data->{gatk}->{run}) ? $tool_data->{gatk}->{run} : 'N',
		'haplotype_caller' => defined($tool_data->{haplotype_caller}->{run}) ? $tool_data->{haplotype_caller}->{run} : 'N',
		'star_fusion'	=> defined($tool_data->{star_fusion}->{run}) ? $tool_data->{star_fusion}->{run} : 'N',
		'arriba'	=> defined($tool_data->{arriba}->{run}) ? $tool_data->{arriba}->{run} : 'N',
		'fusioncatcher'	=> defined($tool_data->{fusioncatcher}->{run}) ? $tool_data->{fusioncatcher}->{run} : 'N',
		'mavis'	=> defined($tool_data->{mavis}->{run}) ? $tool_data->{mavis}->{run} : 'N'
		);

	print $log Dumper \%tool_set;
	print $log "\n";

	# indicate YAML files for processed BAMs
	my $star_output_yaml = join('/', $star_directory, 'star_bam_config_' . $timestamp . '.yaml');
	my $gatk_output_yaml = join('/', $gatk_directory, 'gatk_bam_config_' . $timestamp . '.yaml');

	## run FASTQC pipeline
	unless(-e $fastqc_directory) { make_path($fastqc_directory); }

	if ('Y' eq $tool_set{'fastqc'}) {

		my $fastqc_command = join(' ',
			"perl $cwd/scripts/collect_fastqc_metrics.pl",
			"-o", $fastqc_directory,
			"-t", $tool_config,
			"-d", $data_config,
			"-c", $args{cluster}
			);

		# record command (in log directory) and then run job
		print $log "Submitting job for collect_fastqc_metrics.pl\n";
		print $log "  COMMAND: $fastqc_command\n\n";

		$run_script = write_script(
			log_dir	=> $log_directory,
			name	=> 'pughlab_rna_pipeline__run_fastqc',
			cmd	=> $fastqc_command,
			modules	=> [$perl],
			mem		=> '256M',
			max_time	=> '48:00:00',
			hpc_driver	=> $args{cluster},
			extra_args	=> [$hpc_group]
			);

		if ($args{dry_run}) {

			$fastqc_command .= " --dry-run";
			`$fastqc_command`;
			$fastqc_run_id = 'pughlab_rna_pipeline__run_fastqc';

			} else {

			$fastqc_run_id = submit_job(
				jobname		=> $log_directory,
				shell_command	=> $run_script,
				hpc_driver	=> $args{cluster},
				dry_run		=> $args{dry_run},
				log_file	=> $log
				);

			print $log ">>> FASTQC job id: $fastqc_run_id\n\n";
			push @job_ids, $fastqc_run_id;
			}
		}

	## run FusionCatcher pipeline
	unless(-e $fc_directory) { make_path($fc_directory); }

	if ('Y' eq $tool_set{'fusioncatcher'}) {

		my $fc_command = join(' ',
			"perl $cwd/scripts/fusioncatcher.pl",
			"-o", $fc_directory,
			"-t", $tool_config,
			"-d", $data_config,
			"-c", $args{cluster}
			);

		if ($args{cleanup}) {
			$fc_command .= " --remove";
			}

		# record command (in log directory) and then run job
		print $log "Submitting job for fusioncatcher.pl\n";
		print $log "  COMMAND: $fc_command\n\n";

		$run_script = write_script(
			log_dir	=> $log_directory,
			name	=> 'pughlab_rna_pipeline__run_fusioncatcher',
			cmd	=> $fc_command,
			modules	=> [$perl],
			mem		=> '256M',
			max_time	=> '5-00:00:00',
			extra_args	=> [$hpc_group],
			hpc_driver	=> $args{cluster}
			);

		if ($args{dry_run}) {

			$fc_command .= " --dry-run";
			`$fc_command`;

			} else {

			$fc_run_id = submit_job(
				jobname		=> $log_directory,
				shell_command	=> $run_script,
				hpc_driver	=> $args{cluster},
				dry_run		=> $args{dry_run},
				log_file	=> $log
				);

			print $log ">>> FUSIONCATCHER job id: $fc_run_id\n\n";
			push @job_ids, $fc_run_id;
			}
		}

	## run STAR-alignment pipeline
	unless(-e $star_directory) { make_path($star_directory); }

	if ('Y' eq $tool_set{'star'}) {

		my $star_command = join(' ',
			"perl $cwd/scripts/star.pl",
			"-o", $star_directory,
			"-t", $tool_config,
			"-d", $data_config,
			"-b", $star_output_yaml,
			"-c", $args{cluster}
			);

		if ($args{cleanup}) {
			$star_command .= " --remove";
			}

		# record command (in log directory) and then run job
		print $log "Submitting job for star.pl\n";
		print $log "  COMMAND: $star_command\n\n";

		$run_script = write_script(
			log_dir	=> $log_directory,
			name	=> 'pughlab_rna_pipeline__run_star',
			cmd	=> $star_command,
			modules	=> [$perl],
			mem		=> '256M',
			max_time	=> '5-00:00:00',
			extra_args	=> [$hpc_group],
			hpc_driver	=> $args{cluster}
			);

		if ($args{dry_run}) {

			$star_command .= " --dry-run";
			`$star_command`;

			} else {

			$star_run_id = submit_job(
				jobname		=> $log_directory,
				shell_command	=> $run_script,
				hpc_driver	=> $args{cluster},
				dry_run		=> $args{dry_run},
				log_file	=> $log
				);

			print $log ">>> STAR job id: $star_run_id\n\n";
			push @job_ids, $star_run_id;
			}
		}

	## run STAR-Fusion pipeline
	unless(-e $starfus_directory) { make_path($starfus_directory); }

	if ('Y' eq $tool_set{'star_fusion'}) {

		my $starfus_command = join(' ',
			"perl $cwd/scripts/star_fusion.pl",
			"-o", $starfus_directory,
			"-t", $tool_config,
			"-d", $star_output_yaml,
			"-c", $args{cluster}
			);

		if ($args{cleanup}) {
			$starfus_command .= " --remove";
			}

		# record command (in log directory) and then run job
		print $log "Submitting job for star_fusion.pl\n";
		print $log "  COMMAND: $starfus_command\n\n";

		$run_script = write_script(
			log_dir	=> $log_directory,
			name	=> 'pughlab_rna_pipeline__run_star_fusion',
			cmd	=> $starfus_command,
			modules	=> [$perl],
			dependencies	=> $star_run_id,
			mem		=> '256M',
			max_time	=> '5-00:00:00',
			extra_args	=> [$hpc_group],
			hpc_driver	=> $args{cluster}
			);

		if ($args{dry_run}) {

			$starfus_command .= " --dry-run";
			`$starfus_command`;

			} else {

			$starfus_run_id = submit_job(
				jobname		=> $log_directory,
				shell_command	=> $run_script,
				hpc_driver	=> $args{cluster},
				dry_run		=> $args{dry_run},
				log_file	=> $log
				);

			print $log ">>> STAR-FUSION job id: $starfus_run_id\n\n";
			push @job_ids, $starfus_run_id;
			}
		}

	## run Arriba pipeline
	if ('Y' eq $tool_set{'arriba'}) {

		unless(-e $arriba_directory) { make_path($arriba_directory); }

		my $arriba_command = join(' ',
			"perl $cwd/scripts/arriba.pl",
			"-o", $arriba_directory,
			"-t", $tool_config,
			"-d", $data_config,
			"-c", $args{cluster}
			);

		if ($args{cleanup}) {
			$arriba_command .= " --remove";
			}

		# record command (in log directory) and then run job
		print $log "Submitting job for arriba.pl\n";
		print $log "  COMMAND: $arriba_command\n\n";

		$run_script = write_script(
			log_dir	=> $log_directory,
			name	=> 'pughlab_rna_pipeline__run_arriba',
			cmd	=> $arriba_command,
			modules	=> [$perl],
			mem		=> '256M',
			max_time	=> '5-00:00:00',
			extra_args	=> [$hpc_group],
			hpc_driver	=> $args{cluster}
			);

		if ($args{dry_run}) {

			$arriba_command .= " --dry-run";
			`$arriba_command`;

			} else {

			$arriba_run_id = submit_job(
				jobname		=> $log_directory,
				shell_command	=> $run_script,
				hpc_driver	=> $args{cluster},
				dry_run		=> $args{dry_run},
				log_file	=> $log
				);

			print $log ">>> ARRIBA job id: $arriba_run_id\n\n";
			push @job_ids, $arriba_run_id;
			}
		}

	## run Mavis pipeline
	unless(-e $mavis_directory) { make_path($mavis_directory); }

	if ('Y' eq $tool_set{'mavis'}) {

		my $mavis_command = join(' ',
			"perl $cwd/scripts/mavis.pl",
			"-o", $mavis_directory,
			"-t", $tool_config,
			"-r", $star_output_yaml,
			"-c", $args{cluster}
			);

		# because mavis.pl will search provided directories and run based on what it finds
		# (Manta/Delly/NovoBreak/Pindel/SViCT resuts), this can only be run AFTER these
		# respective jobs finish
		my @depends;
		if (defined($fc_run_id)) {
			$mavis_command .= " --fusioncatcher $fc_directory";
			push @depends, $fc_run_id;
			}
		if (defined($starfus_run_id)) {
			$mavis_command .= " --starfus $starfus_directory";
			push @depends, $starfus_run_id;
			}
		if (defined($arriba_run_id)) {
			$mavis_command .= " --arriba $arriba_directory";
			push @depends, $arriba_run_id;
			}

		if ($args{cleanup}) {
			$mavis_command .= " --remove";
			}

		# record command (in log directory) and then run job
		print $log "Submitting job for mavis.pl\n";
		print $log "  COMMAND: $mavis_command\n\n";

		$run_script = write_script(
			log_dir	=> $log_directory,
			name	=> 'pughlab_rna_pipeline__run_mavis',
			cmd	=> $mavis_command,
			modules	=> [$perl],
			dependencies	=> join(':', @depends),
			mem		=> '256M',
			max_time	=> '5-00:00:00',
			extra_args	=> [$hpc_group],
			hpc_driver	=> $args{cluster}
			);

		if ($args{dry_run}) {
			$mavis_command .= 'pughlab_rna_pipeline__run_mavis';
			} else {
			$mavis_run_id = submit_job(
				jobname		=> $log_directory,
				shell_command	=> $run_script,
				hpc_driver	=> $args{cluster},
				dry_run		=> $args{dry_run},
				log_file	=> $log
				);
			print $log ">>> MAVIS job id: $mavis_run_id\n\n";
			push @job_ids, $mavis_run_id;
			}
		}

	## run RSEM pipeline
	unless(-e $rsem_directory) { make_path($rsem_directory); }

	if ('Y' eq $tool_set{'rsem'}) {

		my $rsem_command = join(' ',
			"perl $cwd/scripts/rsem.pl",
			"-o", $rsem_directory,
			"-t", $tool_config,
			"-d", $star_output_yaml,
			"-c", $args{cluster}
			);

		if ($args{cleanup}) {
			$rsem_command .= " --remove";
			}

		# record command (in log directory) and then run job
		print $log "Submitting job for rsem.pl\n";
		print $log "  COMMAND: $rsem_command\n\n";

		$run_script = write_script(
			log_dir	=> $log_directory,
			name	=> 'pughlab_rna_pipeline__run_rsem',
			cmd	=> $rsem_command,
			modules	=> [$perl],
			dependencies	=> $star_run_id,
			mem		=> '256M',
			max_time	=> '5-00:00:00',
			extra_args	=> [$hpc_group],
			hpc_driver	=> $args{cluster}
			);

		if ($args{dry_run}) {

			$rsem_command .= " --dry-run";
			`$rsem_command`;

			} else {

			$rsem_run_id = submit_job(
				jobname		=> $log_directory,
				shell_command	=> $run_script,
				hpc_driver	=> $args{cluster},
				dry_run		=> $args{dry_run},
				log_file	=> $log
				);

			print $log ">>> RSEM job id: $rsem_run_id\n\n";
			push @job_ids, $rsem_run_id;
			}
		}

	## run GATK indel realignment/recalibration pipeline
	unless(-e $gatk_directory) { make_path($gatk_directory); }

	if ('Y' eq $tool_set{'gatk'}) {

		my $gatk_command = join(' ',
			"perl $cwd/scripts/gatk.pl",
			"--rna",
			"-o", $gatk_directory,
			"-t", $tool_config,
			"-d", $star_output_yaml,
			"-b", $gatk_output_yaml,
			"-c", $args{cluster}
			);

		if ($args{cleanup}) {
			$gatk_command .= " --remove";
			}

		# record command (in log directory) and then run job
		print $log "Submitting job for gatk.pl\n";
		print $log "  COMMAND: $gatk_command\n\n";

		$run_script = write_script(
			log_dir	=> $log_directory,
			name	=> 'pughlab_rna_pipeline__run_gatk',
			cmd	=> $gatk_command,
			modules	=> [$perl],
			dependencies	=> $star_run_id,
			mem		=> '256M',
			max_time	=> '5-00:00:00',
			extra_args	=> [$hpc_group],
			hpc_driver	=> $args{cluster}
			);

		if ($args{dry_run}) {

			$gatk_command .= " --dry-run";
			`$gatk_command`;

			} else {

			$gatk_run_id = submit_job(
				jobname		=> $log_directory,
				shell_command	=> $run_script,
				hpc_driver	=> $args{cluster},
				dry_run		=> $args{dry_run},
				log_file	=> $log
				);

			print $log ">>> GATK job id: $gatk_run_id\n\n";
			push @job_ids, $gatk_run_id;
			}
		}

	## run GATK's HaplotypeCaller pipeline
	unless(-e $vc_directory) { make_path($vc_directory); }

	if ('Y' eq $tool_set{'haplotype_caller'}) {

		my $vc_command = join(' ',
			"perl $cwd/scripts/haplotype_caller.pl",
			"--rna",
			"-o", $vc_directory,
			"-t", $tool_config,
			"-d", $gatk_output_yaml,
			"-c", $args{cluster}
			);

		if ($args{cleanup}) {
			$vc_command .= " --remove";
			}

		# record command (in log directory) and then run job
		print $log "Submitting job for haplotype_caller.pl\n";
		print $log "  COMMAND: $vc_command\n\n";

		$run_script = write_script(
			log_dir	=> $log_directory,
			name	=> 'pughlab_rna_pipeline__run_haplotypecaller',
			cmd	=> $vc_command,
			modules	=> [$perl],
			dependencies	=> $gatk_run_id,
			mem		=> '256M',
			max_time	=> '5-00:00:00',
			extra_args	=> [$hpc_group],
			hpc_driver	=> $args{cluster}
			);

		if ($args{dry_run}) {

			$vc_command .= " --dry-run";
			`$vc_command`;

			} else {

			$vc_run_id = submit_job(
				jobname		=> $log_directory,
				shell_command	=> $run_script,
				hpc_driver	=> $args{cluster},
				dry_run		=> $args{dry_run},
				log_file	=> $log
				);

			print $log ">>> HaplotypeCaller job id: $vc_run_id\n\n";
			push @job_ids, $vc_run_id;
			}
		}

	########################################################################################
	# Create a final report for the project.
	########################################################################################
	if ($args{report}) {

		my $report_command = join(' ',
			"perl $cwd/scripts/pughlab_pipeline_auto_report.pl",
			"-t", $tool_config,
			"-c", $args{cluster},
			"-d", $date
			);
	
		# record command (in log directory) and then run job
		print $log "Submitting job for pughlab_pipeline_auto_report.pl\n";
		print $log "  COMMAND: $report_command\n\n";

		$run_script = write_script(
			log_dir	=> $log_directory,
			name	=> 'pughlab_rna_pipeline__run_report',
			cmd	=> $report_command,
			modules	=> [$perl],
			dependencies	=> join(':', @job_ids),
			mem		=> '256M',
			max_time	=> '24:00:00',
			extra_args	=> [$hpc_group],
			hpc_driver	=> $args{cluster}
			);

		if ($args{dry_run}) {

			$report_command .= " --dry-run";

			} else {

			$report_run_id = submit_job(
				jobname		=> $log_directory,
				shell_command	=> $run_script,
				hpc_driver	=> $args{cluster},
				dry_run		=> $args{dry_run},
				log_file	=> $log
				);

			print $log ">>> Report job id: $report_run_id\n\n";
			}
		}

	# finish up
	print $log "\nProgramming terminated successfully.\n\n";
	close $log;

	}

### GETOPTS AND DEFAULT VALUES #####################################################################
# declare variables
my ($tool_config, $data_config, $summarize);
my $hpc_driver = 'slurm';
my ($remove_junk, $dry_run);

my $help;

# read in command line arguments
GetOptions(
	'h|help'	=> \$help,
	't|tool=s'	=> \$tool_config,
	'd|data=s'	=> \$data_config,
	'c|cluster=s'	=> \$hpc_driver,
	'remove'	=> \$remove_junk,
	'summarize'	=> \$summarize,
	'dry-run'	=> \$dry_run
	 );

if ($help) {
	my $help_msg = join("\n",
		"Options:",
		"\t--help|-h\tPrint this help message",
		"\t--data|-d\t<string> data config (yaml format)",
		"\t--tool|-t\t<string> tool config (yaml format)",
		"\t--cluster|-c\t<string> cluster scheduler (default: slurm)",
		"\t--remove\t<boolean> should intermediates be removed? (default: false)",
		"\t--summarize\t<boolean> should output be summarized and a final report be generated? (default: false)",
		"\t--dry-run\t<boolean> should jobs be submitted? (default: false)"
		);

	print "$help_msg\n";
	exit;
	}

if (!defined($tool_config)) { die("No tool config file defined; please provide -t | --tool (ie, tool_config.yaml)"); }
if (!defined($data_config)) { die("No data config file defined; please provide -d | --data (ie, sample_config.yaml)"); }

# check for compatible HPC driver; if not found, change dry_run to Y
my @compatible_drivers = qw(slurm);
if ( (!any { /$hpc_driver/ } @compatible_drivers ) && (!$dry_run) ) {
	print "Unrecognized HPC driver requested: setting dry_run to true -- jobs will not be submitted but commands will be written to file.\n";
	$dry_run = 1;
	}

main(
	tool_config	=> $tool_config,
	data_config	=> $data_config,
	cluster		=> $hpc_driver,
	cleanup		=> $remove_junk,
	report		=> $summarize,
	dry_run		=> $dry_run
	);