File Coverage

blib/lib/FASTX/Seq.pm

Criterion	Covered	Total	%
statement	197	209	94.2
branch	73	100	73.0
condition	34	54	62.9
subroutine	34	36	94.4
pod	29	29	100.0
total	367	428	85.7

package FASTX::Seq;
#ABSTRACT: A class for representing a sequence for FASTX::Reader

use 5.012;
use warnings;
use Carp qw(confess);
use Data::Dumper;
$Data::Dumper::Sortkeys = 1;
use File::Basename;

$FASTX::Seq::VERSION           = $FASTX::Reader::VERSION;
$FASTX::Seq::DEFAULT_QUALITY   = 'I';
$FASTX::Seq::DEFAULT_LINE_LEN  = 0;
$FASTX::Seq::DEFAULT_OFFSET    = 33;
$FASTX::Seq::QUIET             = $ENV{'FASTX_QUIET'} // 0;
require Exporter;
our @ISA = qw(Exporter);



sub new {
    my $class = shift @_;
    my ($seq, $name, $comment, $qual, $offset, $line_len, $default_quality);
    if (not defined $_[0]) {
        confess "ERROR FASTX::Seq: Sequence missing, record cannot be created.\n";
    }
    # Descriptive instantiation with parameters -param => value
    if (substr($_[0], 0, 1) eq '-') {
        my %data = @_;
        # Try parsing
        for my $i (keys %data) {
            if ($i =~ /^-(seq|sequence)/) {
                $seq = $data{$i};
            } elsif ($i =~ /^-(name|id)/) {
                $name = $data{$i};
            } elsif ($i =~ /^-(comment|desc|description)/) {
                $comment = $data{$i};
            } elsif ($i =~ /^-qual(ity)?/) {
                $qual = $data{$i};
            } elsif ($i =~ /^-offset/) {
                $offset = $data{$i};
            } elsif ($i =~ /^-line_len/) {
                $line_len = $data{$i};
            } elsif ($i =~ /^-default_quality/) {
                $default_quality = $data{$i};
            } else {
                confess "ERROR FASTX::Seq: Unknown parameter $i\n";
            }
        }
    } elsif (not defined $seq) {
        # Positional instantiation
        # check number of arguments
        confess "ERROR FASTX::Seq: Wrong number of arguments\n" if (scalar(@_) < 1 || scalar(@_) > 4);
        ($seq, $name, $comment, $qual) = @_;
    }

    # Required NOT empty
    if (not defined $seq) {
        confess "ERROR FASTX::Seq: Sequence missing, record cannot be created\n";
    }

    # If quality is provided, must match sequence length
    if ( defined $seq && defined $qual
      && (length($seq) != length($qual))) {
        confess "Sequence/quality length mismatch";
    }

    my $self = bless {}, $class;


    $self->{name} = $name   // undef;
    $self->{seq}  = $seq;
    $self->{comment} = $comment // undef;
    $self->{qual} = $qual // undef;

    # Store defaults
    $self->{default_quality} = $default_quality // $FASTX::Seq::DEFAULT_QUALITY;
    $self->{line_len} = $line_len // $FASTX::Seq::DEFAULT_LINE_LEN;
    $self->{offset} = $offset // $FASTX::Seq::DEFAULT_OFFSET;

    return $self;

}


sub copy {
    my ($self) = @_;
    my $copy = __PACKAGE__->new($self->seq, $self->name, $self->comment, $self->qual);
    $copy->{default_quality} = $self->default_quality;
    $copy->line_len = $self->line_len;
    $copy->offset = $self->offset;
    return $copy;
}


sub seq : lvalue {
    # Update sequence
    my ($self, $new_val) = @_;
    $self->{seq} = $new_val if (defined $new_val);
    return $self->{seq};
}



sub name : lvalue {
    # Update name
    my ($self, $new_val) = @_;
    $self->{seq} = $new_val if (defined $new_val);
    return $self->{name};
}



sub qual : lvalue {
    # Update quality
    my ($self, $new_val) = @_;
    $self->{qual} = $new_val if (defined $new_val);
    return $self->{qual};
}


sub comment : lvalue {
    # Update comment
    my ($self, $new_val) = @_;
    $self->{comment} = $new_val if (defined $new_val);
    return $self->{comment};
}


sub offset : lvalue {
    # Update offset
    my ($self, $new_val) = @_;
    confess "ERROR FASTX::Seq: offset must be a positive integer" if (defined $new_val && $new_val !~ /^\d+$/);
    $self->{offset} = $new_val if (defined $new_val);
    return $self->{offset};
}


sub line_len : lvalue {
    # Update line_len
    my ($self, $new_val) = @_;
    confess "ERROR FASTX::Seq: line_len must be a positive integer" if (defined $new_val && $new_val !~ /^\d+$/);
    $self->{line_len} = $new_val if (defined $new_val);
    return $self->{line_len};
}


sub default_quality : lvalue {
    # Update default_quality
    my ($self, $new_val) = @_;
    confess "ERROR FASTX::Seq: default_quality must be a single character" if (defined $new_val && length($new_val) != 1);
    $self->{default_quality} = $new_val if (defined $new_val);
    return $self->{default_quality};
}


sub len {
    # Update comment
    my ($self) = @_;
    return length($self->{seq});
}





sub rev {
    # Update comment
    my ($self) = @_;
    $self->{seq} = reverse($self->{seq});
    $self->{qual} = reverse($self->{qual}) if (defined reverse($self->{qual}));
    return $self;
}



sub rc {
    # Update comment
    my ($self) = @_;
    $self->rev();
    if ($self->{seq} =~ /U/i) {
        $self->{seq} =~ tr/ACGURYSWKMBDHVacguryswkmbdhv/UGCAYRSWMKVHDBugcayrswmkvhdb/;
    } else {
        $self->{seq} =~ tr/ACGTRYSWKMBDHVacgtryswkmbdhv/TGCAYRSWMKVHDBtgcayrswmkvhdb/;
    }
    #$self->{qual} = reverse($self->{qual}) if (defined reverse($self->{qual}));
    return $self;
}


sub slice {
    # Update comment
    my ($self, $from, $len) = @_;
    my $new_seq;
    my $new_qual;
    if (defined($len)) {
      $new_seq = substr($self->{seq}, $from, $len);
      $new_qual = defined($self->{qual}) ? substr($self->{qual}, $from, $len) : undef;
    } else {
        $new_seq = substr($self->{seq}, $from);
        $new_qual = defined($self->{qual}) ? substr($self->{qual}, $from) : undef;
    }
    return __PACKAGE__->new($new_seq, $self->{name}, $self->{comment}, $new_qual);
}


sub char2qual {
    my ($self, $quality_encoded, $offset) = @_;
    # Check quality_encoded is a single character
    confess "Quality encoded character must be a single character" if (length($quality_encoded) != 1);
    $offset = defined $offset ? $offset : $FASTX::Seq::DEFAULT_OFFSET;
    return unpack("C*", $quality_encoded) - $offset;
}


sub qual2char {
    my ($self, $quality_integer, $offset) = @_;
    confess "Quality integer must be an integer value" if ($quality_integer !~ /^\d+$/);
    $offset = defined $offset ? $offset : $FASTX::Seq::DEFAULT_OFFSET;
    return chr($quality_integer + $offset);
}



sub qualities {
    my ($self) = @_;
    my @qualities;
    if (defined $self->{qual}) {
        for (my $i = 0; $i < length($self->{qual}); $i++) {
            push @qualities, $self->char2qual(substr($self->{qual}, $i, 1), $self->{offset});
        }
    }
    return @qualities;
}


sub min_qual {
    my ($self) = @_;
    my @qualities = $self->qualities();
    # Calculate minimum quality
    @qualities = sort {$a <=> $b} @qualities;
    return $qualities[0];
}



sub max_qual {
    my ($self) = @_;
    my @qualities = $self->qualities();
    # Calculate minimum quality
    @qualities = sort {$b <=> $a} @qualities;
    return $qualities[0];
}


sub trim_after  : lvalue {
    my ($self, $qual_value) = @_;
    confess "Quality integer must be an integer value" if ($qual_value !~ /^\d+$/);
    # Detect first base with quality lower or equal than $qual_value
    my $i = 0;
    my @qualities = $self->qualities();
    for ($i = 0; $i < scalar(@qualities); $i++) {
        if ($qualities[$i] <= $qual_value) {
            last;
        }
    }
    # Trim sequence and quality
    my $slice = $self->slice(0, $i);
    $self->{seq} = $slice->{seq};
    $self->{qual} = $slice->{qual};
    return $self;
}



sub trim_until  : lvalue {
    my ($self, $qual_value) = @_;
    confess "Quality integer must be an integer value" if ($qual_value !~ /^\d+$/);
    # Detect first base with quality lower or equal than $qual_value
    my $i = 0;
    my @qualities = $self->qualities();
    for ($i = 0; $i < scalar(@qualities); $i++) {
        if ($qualities[$i] >= $qual_value) {
            last;
        }
    }
    # Trim sequence and quality
    my $slice = $self->slice($i);
    $self->{seq} = $slice->{seq};
    $self->{qual} = $slice->{qual};

    return $self;
}

sub _kmer2num {
    my ($kmer) = @_;

    my %baseVal = ('T' => 0, 'C' => 1, 'A' => 2, 'G' => 3, 'U' => 0);
    my $klen = length($kmer);
    my $num = 0;

    for my $i (0..($klen - 1)) {
        if (exists $baseVal{substr($kmer, $i, 1)}) {
            my $p = 4 ** ($klen - 1 - $i);
            $num += $p * $baseVal{substr($kmer, $i, 1)};
        } else {
            $num = -1;
            last;
        }
    }

    return $num;
}


sub translate {
    my ($self, $code_number) = @_;

    my $record = $self->copy();
    # Default genetic code if not specified
    $code_number //= 11;

    my @code_map = (
     "FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG", #1
     "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSS**VVVVAAAADDEEGGGG",
     "FFLLSSSSYY**CCWWTTTTPPPPHHQQRRRRIIMMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
     "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
     "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSSSSVVVVAAAADDEEGGGG",
     "FFLLSSSSYYQQCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
     "",
     "",
     "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNNKSSSSVVVVAAAADDEEGGGG",
     "FFLLSSSSYY**CCCWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
     "FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG", #11 Bact, Ar and Plant Plast
     "FFLLSSSSYY**CC*WLLLSPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
     "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSSGGVVVVAAAADDEEGGGG",
     "FFLLSSSSYYY*CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNNKSSSSVVVVAAAADDEEGGGG",
     "FFLLSSSSYY*QCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG", #15
     "FFLLSSSSYY*LCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
     "", "", "", "",
     "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNNKSSSSVVVVAAAADDEEGGGG",
     "FFLLSS*SYY*LCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
     "FF*LSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
     "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSSKVVVVAAAADDEEGGGG",
     "FFLLSSSSYY**CCGWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
     "FFLLSSSSYY**CC*WLLLAPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
     "FFLLSSSSYYQQCCWWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
     "FFLLSSSSYYQQCCWWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
     "FFLLSSSSYYYYCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
     "FFLLSSSSYYEECC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG",
     "FFLLSSSSYYEECCWWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG"
    );

    my $code = $code_map[$code_number - 1];
    # Check sequence is nucleotidic or RNA
    confess "ERROR FASTX::Seq: Sequence must be nucleotidic" if ($self->{seq} !~ /^[ACGTU]+$/i);
    my @nucleotides = split //, $self->seq;
    my @translated_sequence;

    for (my $i = 0; $i < scalar(@nucleotides) - 2; $i += 3) {
        my $codon = $nucleotides[$i] . $nucleotides[$i+1] . $nucleotides[$i+2];
        my $num = _kmer2num($codon);  # You need to define this function
        if ($num != -1) {
            push @translated_sequence, substr($code, $num, 1);

        } else {
            push @translated_sequence, '-';

        }
    }

    my $translated_seq = join('', @translated_sequence);

    $record->seq($translated_seq);
    # Remove qual if present
    $record->qual(undef);
    return $record;
}


sub string {
    # Update comment
    my ($self, @args) = @_;
    if (defined $self->{qual}) {
        return $self->asfastq(@args);
    } else {
        return $self->asfasta(@args);
    }

}


sub as_string {
    # Update comment
    my ($self) = @_;
    return $self->string();
}



sub asfasta {
    # Update comment
    my ($self, $len) = @_;

    my $name  = $self->{name} // "sequence";
    my $comment = (defined $self->{comment} and length($self->{comment}) > 0) ? " " . $self->{comment} : "";
    return ">" . $name . $comment . "\n" . _split_string($self->{seq}, $len) . "\n";
}


sub asfastq {
    # Update comment
    my ($self, $user_quality) = @_;
    my $quality = $self->{qual};

    if (defined $user_quality) {
        # User requests new quality
        if (length($user_quality) == 1 ) {
            # And it's valid!
            $quality = $user_quality x length($self->{seq});
        } else {
            if (not $FASTX::Seq::QUIET) {
              say STDERR "[WARNING] FASTX::Seq->as_fastq(): Provide a _char_ as quality, not a value (", $user_quality,"): defaulting to $FASTX::Seq::DEFAULT_QUALITY";
            }
            $quality = $FASTX::Seq::DEFAULT_QUALITY x length($self->{seq});
        }
    } elsif (not defined $quality) {
        $quality = $FASTX::Seq::DEFAULT_QUALITY x length($self->{seq});
    }


    my $name  = $self->{name} // "sequence";
    my $comment = (defined $self->{comment} and length($self->{comment}) > 0) ? " " . $self->{comment} : "";

    return "@" . $name . $comment . "\n" . $self->{seq} . "\n+\n" . $quality  . "\n";
}




sub as_fasta {
    # Update comment
    my ($self, @args) = @_;
    return $self->asfasta(@args);
}


sub as_fastq {
    # Update comment
    my ($self, @args) = @_;
    return $self->asfastq(@args);
}


sub is_fasta {
    # Return true if record has no quality
    my ($self) = @_;
    return ((defined $self->{qual}) > 0 and length($self->{qual}) == length($self->{seq})) ? 0 : 1;
}


sub is_fastq {
    # Update comment
    my ($self) = @_;
    return ((defined $self->{qual}) > 0 and length($self->{qual}) == length($self->{seq})) ? 1 : 0;
}


sub _split_string {
	my ($string, $width) = @_;
	if (not defined $width or $width == 0) {
		return "$string";
	}

	my $output_string;
	for (my $i=0; $i
		my $frag = substr($string, $i, $width);
		$output_string.=$frag."\n";
	}
    chomp($output_string);
	return $output_string;
}

1;

__END__