File Coverage

blib/lib/FASTX/Abi.pm

Criterion	Covered	Total	%
statement	169	175	96.5
branch	54	66	81.8
condition	8	15	53.3
subroutine	20	20	100.0
pod	5	5	100.0
total	256	281	91.1

package FASTX::Abi;
use 5.016;
use warnings;
use Carp qw(confess);
use Bio::Trace::ABIF;
use Data::Dumper;
use File::Basename;
use FASTX::sw 'align';
$FASTX::Abi::VERSION = '1.0.0';

use constant DEFAULT_MATRIX => { 'wildcard_match'  => 0,
                                 'match'           => 1,
                                 'mismatch'        => -1,
                                 'gap'             => -2,
                                 'gap_extend'      => 0,
                                 'wildcard'        => 'N',
                                 };

use constant SCORE => 0;
use constant EVENT => 1;
use constant EXTEND => 0;
use constant GAP_SRC => 1;
use constant GAP_TGT => 2;



#ABSTRACT: Read Sanger trace file (chromatograms) in FASTQ format. For traces called with I option, the ambiguities will be split into two sequences to allow usage from NGS tools that usually do not understand IUPAC ambiguities.

our @valid_new_attributes = (
  'filename',   # *REQUIRED* input trace filepath
  'trim_ends',  # bool (default: 1) trim low quality ends
  'wnd',        # int (default: 16) sliding window for quality trim
  'min_qual',   # int (default: 22) threshold for low quality calls
  'bad_bases',  # int (default: 2)  maximum number of low quality bases per window
  'keep_abi',   # bool (default: 0) import the Bio::Trace::ABIF object in FASTX::Abi (otherwise deleted after import)
);
our @valid_obj_attributes = (
  'diff',             # number of ambiguous bases
  'diff_array',       # array of ambiguous bases position
  'sequence_name',    # sequence name from filename
  'instrument',       # Instrument
  'avg_peak_spacing', # Avg Peak Spacing in chromas
  'version',          # version chromatograms
  'chromas',          # Bio::Trace::ABIF object

  'hetero',           # ambiguity
  'seq1',             # Sequence 1 (non ambiguous, allele1)
  'seq2',             # Sequence 2 (non ambiguous, allele2)
  'sequence',         # Sequence, trimmed
  'quality',          # Quality, trimmed
  'raw_sequence',     # Raw sequence
  'raw_quality',      # Raw quality
  'iso_seq',          # Sequence are equal
  'discard',          # Low quality sequence
);

our %iupac = (
 'R' => 'AG',
 'Y' => 'CT',
 'M' => 'CA',
 'K' => 'TG',
 'W' => 'TA',
 'S' => 'CG'
);


sub new {
    # Instantiate object
    my ($class, $args) = @_;

    my $self = {
        filename  => $args->{filename},   # Chromatogram file name
        trim_ends => $args->{trim_ends},  # Trim low quality ends (bool)
        min_qual  => $args->{min_qual},   # Minimum quality
        wnd       => $args->{wnd},        # Window for end trimming
        bad_bases => $args->{bad_bases},  # Number of low qual bases per $window_width
        keep_abi  => $args->{keep_abi},   # Do not destroy $self->{chromas} after use
    };

    #check valid inputs:
    for my $input (sort keys %{ $args } ) {
      if ( ! grep( /^$input$/, @valid_new_attributes ) ) {
        confess("Method new() does not accept \"$input\" attribute. Valid attributes are:\n", join(', ', @valid_new_attributes));
      }
    }

    # CHECK INPUT FILE
    # -----------------------------------
    if (not defined $self->{filename}) {
      confess("ABI file must be provided when creating new object");
    }

    if (not -e $self->{filename}) {
      confess("ABI file not found: ", $self->{filename});
    }
    my $abif;
    my $try = eval
    {
      $abif = Bio::Trace::ABIF->new();
      $abif->open_abif($self->{filename}) or confess "Error in file: ", $self->{filename};
      1;
    };

    if (not $try) {
      confess("Bio::Trace::ABIF was unable to read: ", $self->{filename});
    }
    my $object = bless $self, $class;
    $object->{chromas} = $abif;

    my @ext = ('.abi','.ab1','.ABI','.abI','.AB1','.ab');
    my ($seqname) = basename($self->{filename}, @ext);
    $object->{sequence_name} = $abif->sample_name();

    # DEFAULTS
    # -----------------------------------
    $object->{trim_ends} = 1  unless defined $object->{trim_ends};
    $object->{wnd}       = 10 unless defined $object->{wnd};
    $object->{min_qual}  = 20 unless defined $object->{min_qual};
    $object->{bad_bases} = 4  unless defined $object->{bad_bases};
    $object->{keep_abi}  = 0  unless defined $object->{keep_abi};
    $object->{discard}   = 0;

    # GET SEQUENCE FROM AB1 FILE
    # -----------------------------------
    my $seq = _get_sequence($self);
    if ($self->{keep_abi} == 0) {
      $self->{chromas} = undef;
    }

    #check valid attributes:
    for my $input (sort keys %{ $self} ) {
      # [this is a developer's safety net]
      # uncoverable condition false
      if ( ! grep( /^$input$/, @valid_new_attributes, @valid_obj_attributes ) ) {
        confess("Method new() does not accept \"$input\" attribute. Valid attributes are:\n", join(', ', @valid_new_attributes, @valid_obj_attributes));
      }
    }


    return $object;
}


sub get_fastq {
  my ($self, $name, $quality_value) = @_;

  if (not defined $name) {
    $name = $self->{sequence_name};
  } elsif ($name=~/\s+/) {
    $name =~s/\s+/_/g;
  }

  my $quality = $self->{quality};
  if (defined $quality_value) {
    if ($quality_value =~/^\d+$/ and $quality_value >= 10) {
      my $q = chr(($quality_value <= 93 ? $quality_value : 93) + 33);
      $quality = $q x length($quality);
    } elsif (length($quality_value) == 1) {
      $quality = $quality_value x length($quality);
    } else {
      confess("Supplied quality is neither a valid integer or a single char: <$quality_value>\n");
    }
  }

  my $output = '';
  if ( $self->{iso_seq} ) {
    $output .= '@' . $name . "\n" .
                $self->{seq1} . "\n+\n" .
                $quality . "\n";
  } else {
    $output .= '@' . $name . "_1\n" .
                $self->{seq1} . "\n+\n" .
                $quality . "\n";
    $output .= '@' . $name . "_2\n" .
                $self->{seq2} . "\n+\n" .
                $quality . "\n";
  }
  return $output;
}


sub get_trace_info {
  my $self   = shift;
  my $data;
  $data->{instrument} = $self->{instrument};
  $data->{version}    = $self->{version};
  $data->{avg_peak_spacing} = $self->{avg_peak_spacing};

  return $data;
}


sub rc {
  my $self = shift;
  $self->{seq1} = reverse $self->{seq1};
  $self->{seq2} = reverse $self->{seq2};
  $self->{seq1} =~ tr/ACGTacgt/TGCAtgca/;
  $self->{seq2} =~ tr/ACGTacgt/TGCAtgca/;
  return $self;
}


sub merge {
  my ($self, $other) = @_;
  $other->rc();

  my $aligner = FASTX::sw->new($self->{seq1}, $other->{seq1});
  my ($top, $bars, $bottom) = $aligner->pads;
  my $consensus = "";
  my $pos1 = 0;
  my $pos2 = 0;
  for (my $pos = 0; $pos < length($top); $pos++) {
    my $base1 = substr($top,    $pos, 1);
    my $base2 = substr($bottom, $pos, 1);

    my $pos1++ if ($base1 ne '-');
    my $pos2++ if ($base2 ne '-');

    if ($base1 eq $base2) {
      $consensus .= $base1;
    } elsif ($base1 eq '-' and $base2 ne '-') {
      $consensus .= $base2;
    } elsif ($base1 ne '-' and $base2 eq '-') {
      $consensus .= $base1;
    } elsif ($base1 ne '-' and $base2 ne '-') {
      my $qual1 = substr($self->{quality}, $pos1, 1);
      my $qual2 = substr($other->{quality}, $pos2, 1);
      $consensus .= _ascii_qual($qual1) >= _ascii_qual($qual2) ? lc($base1) : lc($base2);
    }
  }
  # evaluate longest stretch of "|" in $m
  # my $longest = 0;
  # my $longest_start = 0;
  # my $longest_end = 0;
  # my $longest_str = '';
  # my $i = 0;
  # while ($i < length($m)) {
  #   my $str = substr($m, $i, 1);
  #   if ($str eq '|') {
  #     my $start = $i;
  #     my $end = $i;
  #     while ($str eq '|') {
  #       $end++;
  #       $str = substr($m, $end, 1);
  #     }
  #     if ($end - $start > $longest) {
  #       $longest = $end - $start;
  #       $longest_start = $start;
  #       $longest_end = $end;
  #       $longest_str = $str;
  #     }
  #   }
  #   $i++;
  # }

  $self->{merge} = $aligner;
  $self->{top} = $top;
  $self->{bottom} = $bottom;
  $self->{consensus} = $consensus;
  $self->{bars} = $bars;
  return $consensus;

}


sub _get_sequence {
    my $self   = shift;
    my $abif = $self->{chromas};

    $self->{raw_sequence}  = $abif->sequence();

    # Get quality values
    my @qv       = $abif->quality_values();
    # Encode quality in FASTQ chars
    my @fqv      = map {chr(int(($_<=93? $_ : 93)*4/6) + 33)} @qv;

    # FASTQ
    my $q = join('', @fqv);


    $self->{raw_quality} = $q;

    $self->{sequence} = $self->{raw_sequence};
    $self->{quality}  = $self->{raw_quality};

    # Trim
    if ($self->{trim_ends}) {
        #The Sequencing Analysis program determines the clear range of the sequence by trimming bases from the 5' to 3'
        #ends until fewer than 4 bases out of 20 have a quality value less than 20.
        #You can change these parameters by explicitly passing arguments to this method
        #(the default values are $window_width = 20, $bad_bases_threshold = 4, $quality_threshold = 20).
        # Note that Sequencing Analysis counts the bases starting from one, so you have to add one to the return values to get consistent results.

        my ($begin_pos, $end_pos) = $abif->clear_range(
                                    $self->{wnd},
                                    $self->{bad_bases},
                                    $self->{min_qual},

                                   );

# This can be tested with low quality chromatograms
# *TODO* to ask for some bad trace

# uncoverable branch false
# uncoverable condition left
# uncoverable condition right

             if ($begin_pos>0 and $end_pos>0) {
                my $hi_qual_length = $end_pos-$begin_pos+1;
                $self->{sequence} = substr($self->{sequence}, $begin_pos, $hi_qual_length);
                $self->{quality}  = substr($self->{quality} , $begin_pos, $hi_qual_length);
             } else {
                $self->{discard} = 1;
             }
    }

    # Check hetero bases
    if ($self->{sequence}!~/[ACGT][RYMKWS]+[ACGT]/i) {
        $self->{hetero} = 0;
    } else {
        $self->{hetero} = 1;
    }

    # Check
    $self->{diff_array} = ();
    $self->{diff} = 0;
    my $seq1 = '';
    my $seq2 = '';
    for (my $i = 0; $i{sequence}); $i++) {
      my $q0 = substr($self->{quality}, $i, 1);
      my $s0 = substr($self->{sequence}, $i,1);

      # Ambiguity detected:
      if ($iupac{$s0}) {
        my ($base1, $base2) = split //, $iupac{$s0};
        $seq1.=$base1;
        $seq2.=$base2;
        $self->{diff}++;
        push(@{ $self->{diff_array} }, $i);
      } else {
        $seq1.=$s0;
        $seq2.=$s0;

      }
    }
    $self->{seq1} = $seq1;
    $self->{seq2} = $seq2;

    if ($seq1 eq $seq2) {
       $self->{iso_seq} = 1
     } else {
       $self->{iso_seq} = 0;
     }


     $self->{instrument} = $self->{chromas}->official_instrument_name();
     $self->{version}    = $self->{chromas}->abif_version();
     $self->{avg_peak_spacing} = $self->{chromas}->avg_peak_spacing();

}

sub _ascii_qual {
  my ($qual, $offset) = @_;
  $offset = 33 unless defined $offset;
  return ord($qual) - $offset;
}
1;

__END__