Mercurial > repos > dereeper > roary_plots
diff Roary/lib/Bio/Roary/QC/Report.pm @ 0:c47a5f61bc9f draft
Uploaded
author | dereeper |
---|---|
date | Fri, 14 May 2021 20:27:06 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Roary/lib/Bio/Roary/QC/Report.pm Fri May 14 20:27:06 2021 +0000 @@ -0,0 +1,241 @@ +package Bio::Roary::QC::Report; + +# ABSTRACT: generate a report based on kraken output + +=head1 SYNOPSIS + +=cut + +use Moose; +use File::Temp; +use File::Path 'rmtree'; +use Cwd; +use File::Basename; +with 'Bio::Roary::JobRunner::Role'; + +has 'input_files' => ( is => 'ro', isa => 'ArrayRef', required => 1 ); +has 'kraken_exec' => ( is => 'ro', isa => 'Str', default => 'kraken' ); +has 'kraken_report_exec' => ( is => 'ro', isa => 'Str', default => 'kraken-report' ); +has 'kraken_db' => ( is => 'ro', isa => 'Str', required => 1 ); +has 'outfile' => ( is => 'rw', isa => 'Str', default => 'qc_report.csv' ); +has '_kraken_data' => ( is => 'rw', isa => 'ArrayRef', lazy_build => 1 ); +has '_header' => ( is => 'rw', isa => 'Str', lazy_build => 1 ); +has 'kraken_memory' => ( is => 'rw', isa => 'Int', default => 2000 ); + +has '_tmp_directory_obj' => ( is => 'rw', lazy_build => 1 ); +has '_tmp_directory' => ( is => 'rw', lazy_build => 1, isa => 'Str', ); + + +sub _nuc_fasta_filename +{ + my ($self, $gff) = @_; + + my $prefix = basename( $gff, ".gff" ); + my $outfile = $self->_tmp_directory . "/$prefix.fna"; + return $outfile; +} + +sub _extract_nuc_fasta_cmd { + my ($self, $gff) = @_; + my $outfile = $self->_nuc_fasta_filename($gff); + my $cmd = "sed -n '/##FASTA/,//p' $gff | grep -v \'##FASTA\' > $outfile"; + + return $cmd; +} + +sub _extract_nuc_files_from_all_gffs +{ + my ($self) = @_; + my @nuc_files; + my @commands_to_run; + for my $input_file(@{$self->input_files}) + { + push(@nuc_files,$self->_nuc_fasta_filename($input_file)); + push(@commands_to_run,$self->_extract_nuc_fasta_cmd($input_file)); + } + my $kraken_runner_obj = $self->_job_runner_class->new( + commands_to_run => \@commands_to_run, + memory_in_mb => $self->kraken_memory, + verbose => $self->verbose, + cpus => $self->cpus + ); + $kraken_runner_obj->run(); + return \@nuc_files; +} + +sub _kraken_cmd { + my ( $self, $a, $kraken_output ) = @_; + + my $kcmd = $self->kraken_exec . + " --fasta-input ". + " --preload ". + " --db " . $self->kraken_db . + " --output $kraken_output $a > /dev/null 2>&1"; + return $kcmd; +} + +sub _kraken_report_cmd { + my ( $self, $k, $report_output ) = @_; + + my $krcmd = $self->kraken_report_exec . + " --db " . $self->kraken_db . + " $k > $report_output"; + return $krcmd; +} + +sub _kraken_output_filename +{ + my ( $self, $assembly ) = @_; + my $kraken_output = $assembly; + $kraken_output =~ s/fna$/kraken/; + return $kraken_output; +} + +sub _run_kraken_on_nuc_files +{ + my ( $self, $nuc_files ) = @_; + my @kraken_output_files; + my @commands_to_run; + for my $nuc_file(@{$nuc_files}) + { + my $kraken_output = $self->_kraken_output_filename($nuc_file); + push(@kraken_output_files, $kraken_output ); + push(@commands_to_run, $self->_kraken_cmd( $nuc_file, $kraken_output )); + } + + my $kraken_runner_obj = $self->_job_runner_class->new( + commands_to_run => \@commands_to_run, + memory_in_mb => $self->kraken_memory, + verbose => $self->verbose, + cpus => $self->cpus + ); + $kraken_runner_obj->run(); + + for my $filename(@{$nuc_files}) + { + unlink($filename); + } + + return \@kraken_output_files; +} + +sub _kraken_report_output_filename +{ + my ( $self, $assembly ) = @_; + return $assembly.".report"; +} + +sub _run_kraken_report_on_kraken_files +{ + my ( $self, $kraken_files ) = @_; + + my @kraken_report_output_files; + my @commands_to_run; + for my $nuc_file(@{$kraken_files}) + { + my $kraken_output = $self->_kraken_report_output_filename($nuc_file); + push(@kraken_report_output_files, $kraken_output ); + push(@commands_to_run, $self->_kraken_report_cmd( $nuc_file, $kraken_output )); + } + + my $kraken_runner_obj = $self->_job_runner_class->new( + commands_to_run => \@commands_to_run, + memory_in_mb => $self->kraken_memory, + verbose => $self->verbose, + cpus => $self->cpus + ); + $kraken_runner_obj->run(); + for my $filename(@{$kraken_files}) + { + unlink($filename); + } + return \@kraken_report_output_files; +} + +sub _build__kraken_data { + my $self = shift; + my $nuc_files = $self->_extract_nuc_files_from_all_gffs(); + my $kraken_files = $self->_run_kraken_on_nuc_files($nuc_files); + my $kraken_report_files = $self->_run_kraken_report_on_kraken_files( $kraken_files ); + + return $self->_parse_kraken_reports($kraken_report_files); +} + +sub _parse_kraken_reports +{ + my ( $self, $kraken_report_files ) = @_; + + my @report_rows; + for my $kraken_report(@{$kraken_report_files}) + { + push(@report_rows, $self->_parse_kraken_report($kraken_report)); + } + + for my $kraken_report(@{$kraken_report_files}) + { + unlink($kraken_report); + } + + return \@report_rows; +} + +sub _parse_kraken_report { + my ( $self, $kraken_report ) = @_; + + # parse report + open( my $report_fh, '<', $kraken_report ); + + my $sample_name = $kraken_report; + $sample_name =~ s/.report$//; + $sample_name =~ s/.kraken$//; + my($sample_base_name, $dirs, $suffix) = fileparse($sample_name); + + my ( $top_genus, $top_species ); + while ( <$report_fh> ){ + my @parts = split( "\t" ); + chomp @parts; + + $top_genus = $parts[5] if ( (! defined $top_genus) && $parts[3] eq 'G' ); + $top_species = $parts[5] if ( (! defined $top_species) && $parts[3] eq 'S' ); + + last if (defined $top_genus && defined $top_species); + } + close($report_fh); + + $top_genus ||= "not_found"; + $top_genus =~ s/^\s+//g; + $top_species ||= "not_found"; + $top_species =~ s/^\s+//g; + + return [ $sample_base_name, $top_genus, $top_species ]; +} + + +sub _build__header { + return join( ',', ( 'Sample', 'Genus', 'Species' ) ); +} + +sub _build__tmp_directory_obj { + return File::Temp->newdir(DIR => getcwd, CLEANUP => 1 ); +} + +sub _build__tmp_directory { + my $self = shift; + return $self->_tmp_directory_obj->dirname(); +} + +sub report { + my $self = shift; + + open( OUTFILE, '>', $self->outfile ); + print OUTFILE $self->_header . "\n"; + for my $line ( @{ $self->_kraken_data } ){ + print OUTFILE join( ',', @{ $line } ) . "\n"; + } + close OUTFILE; +} + + +__PACKAGE__->meta->make_immutable; +no Moose; +1;