File Coverage

blib/lib/Bio/WebService/LANL/SequenceLocator.pm

Criterion	Covered	Total	%
statement	29	193	15.0
branch	0	78	0.0
condition	0	36	0.0
subroutine	10	26	38.4
pod	1	6	16.6
total	40	339	11.8

use strictures 1;
use utf8;
use 5.018;

=head1 NAME

Bio::WebService::LANL::SequenceLocator - Locate sequences within HIV using LANL's web tool

=head1 SYNOPSIS

    use Bio::WebService::LANL::SequenceLocator;

    my $locator = Bio::WebService::LANL::SequenceLocator->new(
        agent_string => 'Your Organization - [e-mail]',
    );
    my @sequences = $locator->find([
        "agcaatcagatggtcagccaaaattgccctatagtgcagaacatccaggggcaagtggtacatcaggccatatcacctagaactttaaatgca",
    ]);

See L below.

=head1 DESCRIPTION

This library provides simple programmatic access to
L
web tool and is also used to power
L
for the same tool (via L).

Nearly all of the information output by LANL's sequence locator is parsed and
provided by this library, though the results do vary slightly depending on the
base type of the query sequence.  Multiple query sequences can be located at
the same time and results will be returned for all.

Results are extracted from both tab-delimited files provided by LANL as well as
the HTML itself.

=head1 EXAMPLE RESULTS

    # Using @sequences from the SYNOPSIS above
    use JSON;
    print encode_json(\@sequences);

    __END__
    [
       {
          "query" : "sequence_1",
          "query_sequence" : "AGCAATCAGATGGTCAGCCAAAATTGCCCTATAGTGCAGAACATCCAGGGGCAAGTGGTACATCAGGCCATATCACCTAGAACTTTAAATGCA",
          "base_type" : "nucleotide",
          "reverse_complement" : "0",
          "alignment" : "\n Query AGCAATCAGA TGGTCAGCCA AAATTGCCCT ATAGTGCAGA ACATCCAGGG  50\n       ::::::::    ::::::::: ::::: :::: :::::::::: :::::::::: \n  HXB2 AGCAATCA-- -GGTCAGCCA AAATTACCCT ATAGTGCAGA ACATCCAGGG  1208\n\n Query GCAAGTGGTA CATCAGGCCA TATCACCTAG AACTTTAAAT GCA  93\n       :::: ::::: :::::::::: :::::::::: :::::::::: ::: \n  HXB2 GCAAATGGTA CATCAGGCCA TATCACCTAG AACTTTAAAT GCA  1251\n\n  ",
          "hxb2_sequence" : "AGCAATCA---GGTCAGCCAAAATTACCCTATAGTGCAGAACATCCAGGGGCAAATGGTACATCAGGCCATATCACCTAGAACTTTAAATGCA",
          "similarity_to_hxb2" : "94.6",
          "start" : "373",
          "end" : "462",
          "genome_start" : "1162",
          "genome_end" : "1251",
          "polyprotein" : "Gag",
          "region_names" : [
             "Gag",
             "p17",
             "p24"
          ],
          "regions" : [
             {
                "cds" : "Gag",
                "aa_from_protein_start" : [ "125", "154" ],
                "na_from_cds_start" : [ "373", "462" ],
                "na_from_hxb2_start" : [ "1162", "1251" ],
                "na_from_query_start" : [ "1", "93" ],
                "protein_translation" : "SNQMVSQNCPIVQNIQGQVVHQAISPRTLNA"
             },
             {
                "cds" : "p17",
                "aa_from_protein_start" : [ "125", "132" ],
                "na_from_cds_start" : [ "373", "396" ],
                "na_from_hxb2_start" : [ "1162", "1185" ],
                "na_from_query_start" : [ "1", "27" ],
                "protein_translation" : "SNQMVSQNC"
             },
             {
                "cds" : "p24",
                "aa_from_protein_start" : [ "1", "22" ],
                "na_from_cds_start" : [ "1", "66" ],
                "na_from_hxb2_start" : [ "1186", "1251" ],
                "na_from_query_start" : [ "28", "93" ],
                "protein_translation" : "PIVQNIQGQVVHQAISPRTLNA"
             }
          ]
       }
    ]

=cut

package Bio::WebService::LANL::SequenceLocator;

use Moo;
use Data::Dumper;
use HTML::LinkExtor;
use HTML::TableExtract;
use HTML::TokeParser;
use HTTP::Request::Common;
use List::AllUtils qw< pairwise part min max >;

our $VERSION = 20170324;

=head1 METHODS

=head2 new

Returns a new instance of this class.  An optional parameter C
should be provided to identify yourself to LANL out of politeness.  See the
L for an example.

=cut

has agent_string => (
    is      => 'ro',
    lazy    => 1,
    builder => sub { '' },
);

has agent => (
    is      => 'ro',
    lazy    => 1,
    builder => sub {
        require LWP::UserAgent;
        my $self  = shift;
        my $agent = LWP::UserAgent->new(
            agent => join(" ", __PACKAGE__ . "/$VERSION", $self->agent_string),
        );
        $agent->env_proxy;
        return $agent;
    },
);

has lanl_base => (
    is      => 'ro',
    lazy    => 1,
    builder => sub { 'https://www.hiv.lanl.gov' },
);

has lanl_endpoint => (
    is      => 'ro',
    lazy    => 1,
    builder => sub { shift->lanl_base . '/cgi-bin/LOCATE/locate.cgi' },
);

has _bogus_slug => (
    is      => 'ro',
    default => sub { 'BOGUS_SEQ_SO_TABULAR_FILES_ARE_LINKED_IN_OUTPUT' },
);

sub _request {
    my $self = shift;
    my $req  = shift;
    my $response = $self->agent->request($req);

    if (not $response->is_success) {
        warn sprintf "Request failed: %s %s -> %s\n",
            $req->method, $req->uri, $response->status_line;
        return;
    }

    return $response->decoded_content;
}

=head2 find

Takes an array ref of sequence strings.  Sequences may be in amino acids or
nucleotides and mixed freely.  Sequences should not be in FASTA format.

If sequence bases are not clearly nucleotides or clearly amino acids, LANL
seems to default to nucleotides.  This can be an issue for some sequences since
the full alphabet for nucleotides overlaps with the alphabet for amino acids.
To overcome this problem, you may specify C<< base => 'nucleotide' >>
or C<< base => 'amino acid' >> after the array ref of sequences.  This forces
every sequence to be interpreted as nucleotides or amino acids, so you cannot
mix base types in your sequences if you use this option.  C, C, and
C are accepted aliases for C.  C, C, C,
and C are accepted aliases for C.

Returns a list of hashrefs when called in list context, otherwise returns an
arrayref of hashrefs.

See L for the structure of the data returned.

=cut

sub find {
    my ($self, $sequences, %args) = @_;

    my $content = $self->submit_sequences($sequences, %args)
        or return;

    return $self->parse_html($content);
}

sub submit_sequences {
    my ($self, $sequences, %args) = @_;

    if (defined $args{base}) {
        my $base = lc $args{base};
        if ($base =~ /^n(uc(leotides?)?)?$/i) {
            $args{base} = 1;
        } elsif ($base =~ /^(a(mino( acids?)?)?|aa)$/i) {
            $args{base} = 0;
        } else {
            warn "Unknown base type <$args{base}>, ignoring";
            delete $args{base};
        }
    }

    # Submit multiple sequences at once using FASTA
    my $fasta = join "\n", map {
        ("> sequence_$_", $sequences->[$_ - 1])
    } 1 .. @$sequences;

    # LANL only presents the parseable table.txt we want if there's more
    # than a single sequence.  We always add it so we can reliably skip it.
    $fasta .= "\n> " . $self->_bogus_slug . "\n";

    return $self->_request(
        POST $self->lanl_endpoint,
        Content_Type => 'form-data',
        Content      => [
            organism            => 'HIV',
            DoReverseComplement => 0,
            seq_input           => $fasta,
            (defined $args{base}
                ? ( base => $args{base} )
                : ()),
        ],
    );
}

sub parse_html {
    my ($self, $content) = @_;

    # Fetch and parse the two tables provided as links which removes the need
    # to parse all of the HTML.
    my @results = $self->parse_tsv($content);

    # Now parse the table data from the HTML
    my @tables = $self->parse_tables($content);

    # Extract the alignments, parsing the HTML a third time!
    my @alignments = $self->parse_alignments($content);

    unless (@results and @tables and @alignments) {
        warn "Didn't find all three of TSV, tables, and alignments!\n";
        warn "TSV:             ", scalar @results, "\n";
        warn "HTML tables:     ", scalar @tables, "\n";
        warn "HTML alignments: ", scalar @alignments, "\n";
        warn "Content:\n$content\n", "=" x 80, "\n";
        return;
    }

    unless (@results == @tables and @results == @alignments) {
        warn "Tab-delimited results count doesn't match parsed HTML result count.  Bug!\n";
        warn "TSV:             ", scalar @results, "\n";
        warn "HTML tables:     ", scalar @tables, "\n";
        warn "HTML alignments: ", scalar @alignments, "\n";
        warn "Content:\n$content\n", "=" x 80, "\n";
        return;
    }

    @results = pairwise {
        my $new = {
            %$a,
            base_type       => $b->{base_type},
            regions         => $b->{rows},
            region_names    => [ map { $_->{cds} } @{$b->{rows}} ],
        };
        delete $new->{$_} for qw(protein protein_start protein_end);
        $new;
    } @results, @tables;

    @results = pairwise { +{ %$b, %$a } } @results, @alignments;

    # Fill in genome start/end for amino acid sequences
    for my $r (@results) {
        next unless $r->{base_type} eq 'amino acid';

        if ($r->{genome_start} or $r->{genome_end}) {
            warn "Amino acid sequence with genome start/end already?!",
                 " query <$r->{query_sequence}>";
            next;
        }

        $r->{genome_start} = min map { $_->{na_from_hxb2_start}[0] } @{$r->{regions}};
        $r->{genome_end}   = max map { $_->{na_from_hxb2_start}[1] } @{$r->{regions}};
    }

    return wantarray ? @results : \@results;
}

sub parse_tsv {
    my ($self, $content) = @_;
    my @results;
    my %urls;

    my $extract = HTML::LinkExtor->new(
        sub {
            my ($tag, %attr) = @_;
            return unless $tag eq 'a' and $attr{href};
            return unless $attr{href} =~ m{/(table|simple_results)\.txt$};
            $urls{$1} = $attr{href};
        },
        $self->lanl_base,
    );
    $extract->parse($content);

    for my $table_name (qw(table simple_results)) {
        next unless $urls{$table_name};
        my $table = $self->_request(GET $urls{$table_name})
            or next;

        my (@these_results, %seen);
        my @lines  = split "\n", $table;
        my @fields = map {
            s/^SeqName$/query/;         # standard key
            s/(?<=[a-z])(?=[A-Z])/_/g;  # undo CamelCase
            s/ +/_/g;                   # no spaces
            y/A-Z/a-z/;                 # normalize to lowercase
            # Account for the same field twice in the same data table
            if ($seen{$_}++) {
                $_ = /^(start|end)$/
                    ? "protein_$_"
                    : join "_", $_, $seen{$_};
            }
            $_;
        } split "\t", shift @lines;

        for (@lines) {
            my @values = split "\t";
            my %data;
            @data{@fields} = @values;

            next if $data{query} eq $self->_bogus_slug;

            $data{query_sequence} =~ s/\s+//g
                if $data{query_sequence};
            push @these_results, \%data;
        }

        # Merge with existing results, if any
        @results = @results
                 ? pairwise { +{ %$a, %$b } } @results, @these_results
                 : @these_results;
    }

    return @results;
}

sub parse_tables {
    my ($self, $content) = @_;
    my @tables;

    my %columns_for = (
        'amino acid'    => [
            "CDS"                                               => "cds",
            "AA position relative to protein start in HXB2"     => "aa_from_protein_start",
            "AA position relative to query sequence start"      => "aa_from_query_start",
            "AA position relative to polyprotein start in HXB2" => "aa_from_polyprotein_start",
            "NA position relative to CDS start in HXB2"         => "aa_from_cds_start",
            "NA position relative to HXB2 genome start"         => "na_from_hxb2_start",
        ],
        'nucleotide'    => [
            "CDS"                                                   => "cds",
            "Nucleotide position relative to CDS start in HXB2"     => "na_from_cds_start",
            "Nucleotide position relative to query sequence start"  => "na_from_query_start",
            "Nucleotide position relative to HXB2 genome start"     => "na_from_hxb2_start",
            "Amino Acid position relative to protein start in HXB2" => "aa_from_protein_start",
        ],
    );

    for my $base_type (sort keys %columns_for) {
        my ($their_cols, $our_cols) = part {
            state $i = 0;
            $i++ % 2
        } @{ $columns_for{$base_type} };

        my $extract = HTML::TableExtract->new( headers => $their_cols );
        $extract->parse($content);

        # Examine all matching tables
        for my $table ($extract->tables) {
            my %table = (
                coords      => [$table->coords],
                base_type   => $base_type,
                columns     => $our_cols,
                rows        => [],
            );
            for my $row ($table->rows) {
                @$row = map { defined $_ ? s/^\s+|\s*$//gr : $_ } @$row;

                # An empty row with only a sequence string in the first column.
                if (    $row->[0]
                    and $row->[0] =~ /^[A-Za-z]+$/
                    and not grep { defined and length } @$row[1 .. scalar @$row - 1])
                {
                    $table{rows}->[-1]{protein_translation} = $row->[0];
                    next;
                }

                # Not all rows are data, some are informational sentences.
                next if grep { not defined } @$row;

                my %row;
                @row{@$our_cols} =
                    map { ($_ and $_ eq "NA")       ? undef     : $_ }
                    map { ($_ and /(\d+) → (\d+)/)  ? [$1, $2]  : $_ }
                        @$row;

                push @{$table{rows}}, \%row;
            }
            push @tables, \%table
                if @{$table{rows}};
        }
    }

    # Sort by depth, then within each depth by count
    @tables = sort {
        $a->{coords}[0] <=> $b->{coords}[0]
     or $a->{coords}[1] <=> $b->{coords}[1]
    } @tables;

    if (@tables > 1) {
        unless (    $tables[-1]->{rows}[0]{na_from_query_start} eq "1 →"
                and $tables[-1]->{rows}[0]{protein_translation} eq "X") {
            warn "Last table appears to be real!?  It should be the bogus table of the bogus sequence.";
            warn "Table is ", Dumper($tables[-1]), "\n";
            return;
        } else {
            pop @tables;
        }
    }

    return @tables;
}

sub parse_alignments {
    my ($self, $content) = @_;
    my @alignments;

    my $doc = HTML::TokeParser->new(
        \$content,
        unbroken_text => 1,
    );

    my $expect_alignment = 0;

    while (my $tag = $doc->get_tag("b", "pre")) {
        my $name = lc $tag->[0];
        my $text = $doc->get_text;
        next unless defined $text;

        # s are preceeded by a bold header, which we use as an indicator
        if ($name eq 'b') {
            $expect_alignment = $text =~ /Alignment\s+of\s+the\s+query\s+sequence\s+to\s+HXB2/i;
        } elsif ($name eq 'pre') {
            if ($text =~ /^\s*Query\b/m and $text =~ /^\s*HXB2\b/m) {
                push @alignments, $text;
                warn "Not expecting alignment, but found one‽"
                    unless $expect_alignment;
            }
            elsif ($text =~ /^\s+$/ and $expect_alignment) {
                push @alignments, undef;    # We appear to have found an unaligned sequence.
            }
            $expect_alignment = 0;
        }
    }

    if (defined $alignments[-1]) {
        warn "Last alignment is non-null!  It should be the empty alignment of the bogus sequence.";
        warn "Alignment is <$alignments[-1]>\n";
        return;
    } else {
        pop @alignments;
    }

    my @results;
    for (@alignments) {
        my @hxb2;
        if (defined) {
            push @hxb2, $1 =~ s/\s+//gr
                while /^\s*HXB2\b\s+(.+?)(?:\s+\d+|\s*)$/gm;
        }
        push @results, {
            alignment       => $_,
            hxb2_sequence   => @hxb2 ? join("", @hxb2) : undef,
        };
    }

    return @results;
}

=head1 AUTHOR

Thomas Sibley [e-mail]

=head1 COPYRIGHT

Copyright 2014 by the Mullins Lab, Department of Microbiology, University of
Washington.

=head1 LICENSE

Licensed under the same terms as Perl 5 itself.

=cut

42;