New upstream version 5.0.1 (1cc2e64b) · Commits · Debian Med / transdecoder

.gitmodules

+3 −0

Original line number	Diff line number	Diff line
		[submodule "TransDecoder.lrgTests"]
		path = TransDecoder.lrgTests
		url = https://github.com/TransDecoder/TransDecoder.lrgTests.git
		[submodule "TransDecoder.wiki"]
		path = TransDecoder.wiki
		url = https://github.com/TransDecoder/TransDecoder.wiki.git

Changelog.txt

+44 −0

Original line number	Diff line number	Diff line
		## v5.0.0 August 26, 2017

		-algorithm updates: frame[0] score > 0 and max for first 3 reading frames (instead of all 6), and orf with highest frame[0] score is chosen allowing for minimal overlap among selected predictions.
		-option --single_best_only provides the single longest of the selected orfs per contig.
		-long orfs unlikely to appear in random sequence are automatically selected as candidates with this minimal long orf length set dynamically according to GC content.
		-orf score and blast or pfam info is propagated to gff3 output



		## v4.1.0

		-single best orf now selected by default. If more than the single best orf is wanted, use the --all_good_orfs parameter.
		-start codon refinement is now done by default. To turn it off and get the original behavior of extending to the longest orf position, use parameter: --no_refine_starts
		-cdhit has been removed and replaced by our own fast method for removing redundancies.
		-selection of coding regions is strictly governed by Markov-based likelihood scores across reading frames. No auto-retention of long orfs by default, but can be activated by parameter: --retain_long_orfs_length


		all v4 releases pre-v4.1 were fairly quickly retracted due to bugs and insufficient benchmarking


		## v3.0.2 release Oct 31, 2016

		minor bugfix release - when checking for required utilities to be installed, doesn't require ^/ in path



		## v3.0.0 release April 26, 2016

		TransDecoder v3.0.0 includes the following changes:

		TransDecoder.LongOrfs now includes parameter '--gene_trans_map ' as a way to retain the gene identifier information. In the case of Cufflinks and Trinity, the gene identifiers will automatically be recognized and retained. For PASA and other inputs, it is necessary to provide the gene-to-transcript identifier mappings in order to generate isoform-clustered output files (gff3).

		TransDecoder.Predict now includes flag ' --single_best_orf ' to retain only the single 'best' ORF per transcript. ORFs are prioritized according to homology information (if given the blast and pfam results) and by sequence length, with longer ORFs preferred.

		Codon phase information is now included in the GFF3 output files.

		The .mRNA files that were generated by default for genome-free TransDecoder runs are now deprecated, but of course the .cds and .pep files are provided.

		The sample data sets include examples for running TransDecoder in a few different contexts, including starting from Trinity, PASA, or Cufflinks data.

		More useful logging information is provided to it's clearer as to how many orfs are being retained and which are being eliminated along the way.



		## 2016-03-11 v2.1 release
		-added cpu parameter to TransDecoder.predict
		-retaining gene identifier information from cufflinks output

Makefile

+3 −4

Original line number	Diff line number	Diff line
		SHELL := /bin/bash

		all:
		cd ./transdecoder_plugins/ && $(MAKE) all
		@echo Nothing to build. Run \'make test\' to run example executions on each of the sample data sets provided.

		clean:
		cd ./transdecoder_plugins/ && $(MAKE) clean
		cd ./sample_data && ./cleanme.pl
		cd ./sample_data && make clean

		test:
		cd ./sample_data/ && ./runMe.sh
		cd ./sample_data/ && make test

PerlLib/DelimParser.pm

0 → 100644

+278 −0

Original line number	Diff line number	Diff line
		#!/usr/bin/env perl

		# classes for DelimParser::Reader and DelimParser::Writer

		package DelimParser;
		use strict;
		use warnings;
		use Carp;

		####
		sub new {
		my ($packagename, $fh, $delimiter) = @_;

		unless ($fh && $delimiter) {
		confess "Error, need filehandle and delimiter params";
		}


		my $self = { _delim => $delimiter,
		_fh => $fh,

		# set below in _init()
		_column_headers => [],
		};


		bless ($self, $packagename);

		return($self);
		}


		####
		sub get_fh {
		my $self = shift;
		return($self->{_fh});
		}

		####
		sub get_delim {
		my $self = shift;
		return($self->{_delim});
		}

		####
		sub get_column_headers {
		my $self = shift;
		return(@{$self->{_column_headers}});
		}

		####
		sub set_column_headers {
		my $self = shift;
		my (@columns) = @_;

		$self->{_column_headers} = \@columns;

		return;
		}

		####
		sub get_num_columns {
		my $self = shift;
		return(length($self->get_column_headers()));
		}


		###
		sub reconstruct_header_line {
		my $self = shift;
		my @column_headers = $self->get_column_headers();

		my $header_line = join("\t", @column_headers);
		return($header_line);
		}

		###
		sub reconstruct_line_from_row {
		my $self = shift;
		my $row_href = shift;
		unless ($row_href && ref $row_href) {
		confess "Error, must set row_href as param";
		}

		my @column_headers = $self->get_column_headers();

		my @vals;
		foreach my $col_header (@column_headers) {
		my $val = $row_href->{$col_header};
		push (@vals, $val);
		}

		my $row_text = join("\t", @vals);

		return($row_text);

		}


		##################################################
		package DelimParser::Reader;
		use strict;
		use warnings;
		use Carp;
		use Data::Dumper;

		our @ISA;
		push (@ISA, 'DelimParser');

		sub new {
		my ($packagename) = shift;
		my $self = $packagename->DelimParser::new(@_);

		$self->_init();

		return($self);
		}


		####
		sub _init {
		my $self = shift;

		my $fh = $self->get_fh();
		my $delim = $self->get_delim();

		my $header_row = <$fh>;
		chomp $header_row;

		unless ($header_row) {
		confess "Error, no header row read.";
		}

		my @fields = split(/$delim/, $header_row);

		$self->set_column_headers(@fields);


		return;
		}

		####
		sub get_row {
		my $self = shift;

		my $fh = $self->get_fh();
		my $line = <$fh>;
		unless ($line) {
		return(undef); # eof
		}

		my $delim = $self->get_delim();
		my @fields = split(/$delim/, $line);
		chomp $fields[$#fields]; ## it's important that this is done after the delimiter splitting in case the last field is actually empty.


		my @column_headers = $self->get_column_headers();

		my $num_col = scalar (@column_headers);
		my $num_fields = scalar(@fields);

		if ($num_col != $num_fields) {
		confess "Error, line: [$line] " . Dumper(\@fields) . " is lacking $num_col fields: " . Dumper(\@column_headers);
		}

		my %dict;
		foreach my $colname (@column_headers) {
		my $field = shift @fields;
		$dict{$colname} = $field;
		}

		return(\%dict);
		}


		##################################################

		package DelimParser::Writer;
		use strict;
		use warnings;
		use Carp;

		our @ISA;
		push (@ISA, 'DelimParser');

		sub new {
		my ($packagename) = shift;
		my ($ofh, $delim, $column_fields_aref, $FLAGS) = @_;

		## FLAGS can be:
		# NO_WRITE_HEADER\|...

		unless (ref $column_fields_aref eq 'ARRAY') {
		confess "Error, need constructor params: ofh, delim, column_fields_aref";
		}

		my $self = $packagename->DelimParser::new($ofh, $delim);

		$self->_initialize($column_fields_aref, $FLAGS);

		return($self);
		}


		####
		sub _initialize {
		my $self = shift;
		my $column_fields_aref = shift;
		my $FLAGS = shift;

		unless (ref $column_fields_aref eq 'ARRAY') {
		confess "Error, require column_fields_aref as param";
		}


		my $ofh = $self->get_fh();
		my $delim = $self->get_delim();



		$self->set_column_headers(@$column_fields_aref);

		unless (defined($FLAGS) && $FLAGS =~ /NO_WRITE_HEADER/) {
		my $output_line = join($delim, @$column_fields_aref);
		print $ofh "$output_line\n";
		}


		return;
		}


		####
		sub write_row {
		my $self = shift;
		my $dict_href = shift;

		unless (ref $dict_href eq "HASH") {
		confess "Error, need dict_href as param";
		}

		my $num_dict_fields = scalar(keys %$dict_href);

		my @column_headers = $self->get_column_headers();


		my $delim = $self->get_delim();

		my @out_fields;
		for my $column_header (@column_headers) {
		my $field = $dict_href->{$column_header};
		unless (defined $field) {
		confess "Error, missing value for required column field: $column_header";
		}
		if ($field =~ /$delim/) {
		# don't allow any delimiters to contaminate the field value, otherwise it'll introduce offsets.
		$field =~ s/$delim/ /g;
		}
		# also avoid newlines, which will also break the output formatting.
		if ($field =~ /\n/) {
		$field =~ s/\n/ /g;
		}

		push (@out_fields, $field);
		}

		my $outline = join("\t", @out_fields);

		my $ofh = $self->get_fh();

		print $ofh "$outline\n";

		return;
		}


		1; #EOM

PerlLib/PWM.pm

0 → 100644

+413 −0

Original line number	Diff line number	Diff line
		package PWM;

		use strict;
		use warnings;
		use Carp;

		my $PSEUDOCOUNT = 0.1;

		srand(1); # reproducibility

		###
		sub new {
		my $packagename = shift;

		my $self = { pos_freqs => [],
		pos_probs => [],
		_pwm_built_flag => 0,
		};

		bless($self, $packagename);

		return($self);
		}


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

		return($self->{_pwm_built_flag});

		}

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

		if (@{$self->{pos_freqs}}) {
		return(1);
		}
		else {
		return(0);
		}
		}



		sub add_feature_seq_to_pwm {
		my ($self, $feature_seq) = @_;

		$feature_seq = uc $feature_seq;

		my $pwm_len = $self->get_pwm_length();
		if ($pwm_len && length($feature_seq) != $pwm_len) {
		confess "Error, pwm_len: $pwm_len and feature_seq len: " . length($feature_seq) . " are unmatched.";
		}
		if ($feature_seq =~ /[^GATC]/) {
		print STDERR "Error, feature_seq: $feature_seq contains non-GATC chars... skipping\n";
		return;
		}

		my @chars = split(//, $feature_seq);
		for (my $i = 0; $i <= $#chars; $i++) {
		my $char = $chars[$i];

		$self->{pos_freqs}->[$i]->{$char}++;
		}

		return;
		}


		sub remove_feature_seq_from_pwm {
		my ($self, $feature_seq) = @_;

		if (! $self->has_pos_freqs()) {
		confess "Error, pwm obj doesn't have position frequencies set";
		}

		$feature_seq = uc $feature_seq;

		my $pwm_len = $self->get_pwm_length();
		if ($pwm_len && length($feature_seq) != $pwm_len) {
		confess "Error, pwm_len: $pwm_len and feature_seq len: " . length($feature_seq) . " are unmatched.";
		}
		if ($feature_seq =~ /[^GATC]/) {
		print STDERR "Warning, feature_seq: $feature_seq contains non-GATC chars... skipping.\n";
		return;
		}

		my @chars = split(//, $feature_seq);
		for (my $i = 0; $i <= $#chars; $i++) {
		my $char = $chars[$i];
		$self->{pos_freqs}->[$i]->{$char}--;
		}

		return;
		}



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

		my $pwm_length;

		if ($self->is_pwm_built()) {
		$pwm_length = $#{$self->{pos_probs}} + 1;
		}
		else {
		$pwm_length = $#{$self->{pos_freqs}} + 1;
		}

		return($pwm_length);
		}


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

		if (! $self->is_pwm_built()) {
		confess "Error, pwm not built yet";
		}

		my $probs_aref = $self->{pos_probs};

		my @chars = qw(G A T C);
		my $feature_seq = "";

		my $pwm_length = $self->get_pwm_length();
		for (my $i = 0; $i < $pwm_length; $i++) {
		my $sum_prob = 0;
		my $selected_flag = 0;
		#print "rand val: $rand_val\n";
		my $tot_prob = 0;
		for (my $j = 0; $j <= $#chars; $j++) {
		my $char = $chars[$j];
		my $p = $probs_aref->[$i]->{$char};
		$tot_prob += $p;
		}
		my $rand_val = rand($tot_prob); # tot_prob should be 1 or very close...


		for (my $j = 0; $j <= $#chars; $j++) {
		my $char = $chars[$j];
		my $p = $probs_aref->[$i]->{$char};
		#print "char: $char, p: $p\n";
		$sum_prob += $p;
		if ($j == $#chars) {
		$sum_prob = 1;
		}

		if ($rand_val <= $sum_prob) {
		# choose char
		$feature_seq .= $char;
		$selected_flag = 1;
		last;
		}
		}
		if (! $selected_flag) {
		croak "Error, didn't select a random char for feature seq";
		}
		}

		return($feature_seq);

		}



		####
		sub write_pwm_file {
		my ($self, $filename) = @_;

		unless ($self->is_pwm_built()) {
		croak "Error, pwm needs to be built first before writing to file";
		}

		open (my $ofh, ">$filename") or croak "Error, cannot write to file: $filename";

		my $pos_probs_aref = $self->{pos_probs};
		my $pwm_length = $self->get_pwm_length();

		my @chars = qw(G A T C);
		print $ofh join("\t", "pos", @chars) . "\n";
		for (my $i = 0; $i < $pwm_length; $i++) {

		my @vals = ($i);
		foreach my $char (@chars) {
		my $prob = $pos_probs_aref->[$i]->{$char} \|\| 0;
		push (@vals, $prob);
		}

		print $ofh join("\t", @vals) . "\n";
		}

		close $ofh;

		return;
		}



		####
		sub build_pwm {
		my ($self) = @_;

		my $pos_freqs_aref = $self->{pos_freqs};

		for (my $i = 0; $i <= $#{$pos_freqs_aref}; $i++) {
		my @chars = keys %{$pos_freqs_aref->[$i]};

		my $sum = 0;
		foreach my $char (@chars) {
		my $val = $pos_freqs_aref->[$i]->{$char};
		$sum += $val;
		}

		# now convert to relative freqs
		@chars = qw(G A T C); # the complete set of chars we care about.
		foreach my $char (@chars) {
		my $val = $pos_freqs_aref->[$i]->{$char} \|\| 0;
		my $prob = sprintf("%.6f", ($val + $PSEUDOCOUNT) / ($sum + 4 * $PSEUDOCOUNT) );
		$self->{pos_probs}->[$i]->{$char} = $prob;
		}

		}

		$self->{_pwm_built_flag} = 1;

		return;
		}


		sub score_pwm_using_base_freqs {
		my ($self, $target_sequence, $base_freqs_href, %options ) = @_;

		### Options can include:
		###
		### mask => [ coordA, coordB, coordC ], # ignored pwm positions in scoring
		### pwm_range => [pwm_idx_start, pwm_idx_end], # start and end are inclusive
		###

		for my $key (keys %options) {
		if (! grep {$_ eq $key} ('mask', 'pwm_range')) {
		confess "Error, option $key is not recognized";
		}
		}

		#print STDERR "target_seq: [$target_sequence]\n";

		$target_sequence = uc $target_sequence;
		unless ($target_sequence =~ /^[GATC]+$/) {
		# can only score GATC-containing sequences.
		return("NA");
		}

		if (! $self->is_pwm_built()) {
		croak("pwm not built yet!");
		}

		my $pwm_length = $self->get_pwm_length();

		if (length($target_sequence) != $pwm_length) {
		croak "Error, len(target_sequence)=" . length($target_sequence) . " and pwm length = $pwm_length";
		}

		my %mask;
		if (my $mask_positions_aref = $options{'mask'}) {
		%mask = map { + $_ => 1 } @$mask_positions_aref;
		}

		my $motif_score = 0;

		my @seqarray = split(//, $target_sequence);

		my ($pwm_start, $pwm_end) = (0, $pwm_length-1);
		if (my $pwm_range_aref = $options{'pwm_range'}) {
		($pwm_start, $pwm_end) = @$pwm_range_aref;
		}

		for (my $i = $pwm_start; $i <= $pwm_end; $i++) {

		if ($mask{$i}) {
		next;
		}

		my $char = $seqarray[$i];
		my $prob = $self->{pos_probs}->[$i]->{$char};

		unless ($prob) {
		return("NA");
		}

		my $prob_rand = $base_freqs_href->{$char};
		unless ($prob_rand) {
		die "Error, no non-zero probability value specified for char [$char] ";
		}

		my $loglikelihood = log($prob/$prob_rand);
		$motif_score += $loglikelihood;

		}

		return($motif_score);

		}


		sub score_plus_minus_pwm {
		my ($self, $target_sequence, $pwm_minus, %options) = @_;

		### Options can include:
		###
		### mask => [ coordA, coordB, coordC ], # ignored pwm positions in scoring
		### pwm_range => [pwm_idx_start, pwm_idx_end], # start and end are inclusive
		###

		for my $key (keys %options) {
		if (! grep {$_ eq $key} ('mask', 'pwm_range')) {
		confess "Error, option $key is not recognized";
		}
		}

		#print STDERR "target_seq: [$target_sequence]\n";

		$target_sequence = uc $target_sequence;
		unless ($target_sequence =~ /^[GATC]+$/) {
		# can only score GATC-containing sequences.
		return("NA");
		}

		if (! $self->is_pwm_built()) {
		croak("pwm not built yet!");
		}

		my $pwm_length = $self->get_pwm_length();

		if (length($target_sequence) != $pwm_length) {
		croak "Error, len(target_sequence)=" . length($target_sequence) . " and pwm length = $pwm_length";
		}

		my %mask;
		if (my $mask_positions_aref = $options{'mask'}) {
		%mask = map { + $_ => 1 } @$mask_positions_aref;
		}

		my $motif_score = 0;

		my @seqarray = split(//, $target_sequence);

		my ($pwm_start, $pwm_end) = (0, $pwm_length-1);
		if (my $pwm_range_aref = $options{'pwm_range'}) {
		($pwm_start, $pwm_end) = @$pwm_range_aref;
		}

		for (my $i = $pwm_start; $i <= $pwm_end; $i++) {

		if ($mask{$i}) {
		#print STDERR "masking $i\n";
		next;
		}

		my $char = $seqarray[$i];
		my $prob = $self->{pos_probs}->[$i]->{$char};

		unless ($prob) {
		return("NA");
		}

		my $prob_rand = $pwm_minus->{pos_probs}->[$i]->{$char};
		unless ($prob_rand) {
		die "Error, no non-zero probability value specified for char [$char] ";
		}

		my $loglikelihood = log($prob/$prob_rand);
		$motif_score += $loglikelihood;

		}

		return($motif_score);

		}

		####
		sub load_pwm_from_file {
		my ($self, $pwm_file) = @_;

		open(my $fh, $pwm_file) or confess "Error, cannot open file: $pwm_file";

		my $header = <$fh>;
		chomp $header;
		my @chars = split(/\t/, $header);
		shift @chars;

		while (<$fh>) {
		chomp;
		my @x = split(/\t/);
		my $idx_pos = shift @x;
		for (my $i = 0; $i <= $#chars; $i++) {
		my $char = $chars[$i];
		my $prob = $x[$i];

		$self->{pos_probs}->[$idx_pos]->{$char} = $prob;
		}
		}
		close $fh;

		$self->{_pwm_built_flag} = 1;

		return;
		}


		1; #EOM