view release on metacpan or  search on metacpan

lib/CLIPSeqTools/App/reads_long_gaps_size_distribution.pm  view on Meta::CPAN

		-excludes => 'validate_args',
	},
	"CLIPSeqTools::Role::Option::Plot" => {
		-alias    => { validate_args => '_validate_args_for_plot' },
		-excludes => 'validate_args',
	},
	"CLIPSeqTools::Role::Option::OutputPrefix" => {
		-alias    => { validate_args => '_validate_args_for_output_prefix' },
		-excludes => 'validate_args',
	};


#######################################################################
########################   Interface Methods   ########################
#######################################################################
sub validate_args {
	my ($self) = @_;

	$self->_validate_args_for_library;
	$self->_validate_args_for_plot;
	$self->_validate_args_for_output_prefix;
}

sub run {
	my ($self) = @_;

	warn "Starting analysis: reads_long_gaps_size_distribution\n";

	warn "Validating arguments\n" if $self->verbose;
	$self->validate_args();

	warn "Creating reads collection\n" if $self->verbose;
	my $reads_collection = $self->reads_collection;
	$reads_collection->schema->storage->debug(1) if $self->verbose > 1;
	my $total_cn = $reads_collection->total_copy_number;

	warn "Measuring long gaps\n" if $self->verbose;
	my %gap_size_count;
	$reads_collection->foreach_record_do( sub {
		my ($rec) = @_;

		my $cigar = $rec->cigar;
		if ($cigar !~ /N/) {
			$gap_size_count{0} += $rec->copy_number;
		}
		else {
			while ($cigar =~ /(\d+)N/g) {
				my $gap_length = $1;
				$gap_size_count{$gap_length} += $rec->copy_number;
			}
		}

		return 0;
	});

	warn "Creating output path\n" if $self->verbose;
	$self->make_path_for_output_prefix();

	warn "Printing results\n" if $self->verbose;
	open (my $OUT1, '>', $self->o_prefix.'reads_long_gaps_size_distribution.tab');
	say $OUT1 join("\t", 'gap_size', 'count', 'percent');
	say $OUT1 join("\t", $_, $gap_size_count{$_}, ($gap_size_count{$_}/$total_cn)*100) for sort {$a <=> $b} keys %gap_size_count;
	close $OUT1;

	if ($self->plot) {
		warn "Creating plot\n" if $self->verbose;
		CLIPSeqTools::PlotApp->initialize_command_class(
			'CLIPSeqTools::PlotApp::reads_long_gaps_size_distribution',
			file     => $self->o_prefix.'reads_long_gaps_size_distribution.tab',
			o_prefix => $self->o_prefix
		)->run();
	}
}

1;