src/c++/lib/grm/GraphAligner.cpp
/**
 * Align a read to the graph and update the graph_* fields.
 *
 * We will align both the forward and the reverse strand and return
 * the alignment which is unique, or with the better score, or default
 * to the forward strand.
 *
 * @param read read structure
 * @param alignment_flags flags for alignment (see above)
 */
Align a read to the graph and update the graph_* fields.
We will align both the forward and the reverse strand and return
the alignment which is unique, or with the better score, or default
to the forward strand.

@param read read structure
@param alignment_flags flags for alignment
Align a read to the graph and update the graph_ * fields . We will align both the forward and the reverse strand and return the alignment which is unique or with the better score or default to the forward strand . @param read read structure @param alignment_flags flags for alignment
void GraphAligner::alignRead(Read& read, unsigned int alignment_flags) const
{
    GraphAlignerImpl::p_gssw_graph_mapping gm_fwd_strand{ nullptr, safe_gssw_graph_mapping_destroy };
    GraphAlignerImpl::p_gssw_graph_mapping gm_reverse_strand{ nullptr, safe_gssw_graph_mapping_destroy };
    GraphAlignerImpl::p_gssw_graph_mapping rgm_fwd_strand{ nullptr, safe_gssw_graph_mapping_destroy };
    GraphAlignerImpl::p_gssw_graph_mapping rgm_reverse_strand{ nullptr, safe_gssw_graph_mapping_destroy };
    const std::string rev_cmp_bases = reverseComplement(read.bases());
    const bool fwd_strand_multi_align
        = _impl->alignString(_impl->graph_.get(), _impl->nodes_, _impl->node_map_, read.bases(), gm_fwd_strand);
    const bool reverse_strand_multi_align = (alignment_flags & AF_BOTH_STRANDS) != 0u
        ? _impl->alignString(_impl->graph_.get(), _impl->nodes_, _impl->node_map_, rev_cmp_bases, gm_reverse_strand)
        : false;
    bool rfwd_strand_multi_align = false;
    bool rreverse_strand_multi_align = false;
    if ((alignment_flags & AF_REVERSE_GRAPH) != 0u)
    {
        string bases_rev = read.bases();
        std::reverse(bases_rev.begin(), bases_rev.end());
        rfwd_strand_multi_align = _impl->alignString(
            _impl->graph_reversed_.get(), _impl->nodes_reversed_, _impl->node_map_reversed_, bases_rev, rgm_fwd_strand);
        rreverse_strand_multi_align = (alignment_flags & AF_BOTH_STRANDS) != 0u
            ? _impl->alignString(
                  _impl->graph_reversed_.get(), _impl->nodes_reversed_, _impl->node_map_reversed_,
                  reverseComplement(bases_rev), rgm_reverse_strand)
            : false;
    }
    bool fwd_strand_is_unique = (!fwd_strand_multi_align) && (!rfwd_strand_multi_align);
    bool reverse_strand_is_unique = (!reverse_strand_multi_align) && (!rreverse_strand_multi_align);
    bool return_reverse = false;
    if ((!fwd_strand_is_unique) && reverse_strand_is_unique && gm_reverse_strand)
    {
        return_reverse = true;
    }
    else if (fwd_strand_is_unique && (!reverse_strand_is_unique))
    {
        return_reverse = false;
    }
    else if (gm_reverse_strand)
    {
        return_reverse = gm_fwd_strand->score < gm_reverse_strand->score;
    }
    const bool resulting_strand_is_reverse = read.is_reverse_strand() != return_reverse;
    read.set_is_graph_reverse_strand(resulting_strand_is_reverse);
    auto make_node_list = [this](gssw_graph_mapping* gm) -> std::list<int> {
        std::list<int> result;
        auto const& g = gm->cigar;
        auto nc = g.elements;
        for (uint32_t i = 0; i < g.length; ++i, ++nc)
        {
            result.push_back(_impl->node_map_[nc->node->id]);
        }
        return result;
    };
    std::list<int> nodes_passed;
    if (return_reverse)
    {
        read.set_bases(rev_cmp_bases);
        std::string rev_quals = read.quals();
        std::reverse(rev_quals.begin(), rev_quals.end());
        read.set_quals(rev_quals);
        read.set_graph_pos(gm_reverse_strand->position);
        read.set_graph_alignment_score(gm_reverse_strand->score);
        read.set_is_graph_alignment_unique(reverse_strand_is_unique);
        read.set_graph_mapq(reverse_strand_is_unique ? 60 : 0);
        if (alignment_flags & AF_CIGAR)
        {
            read.set_graph_cigar(_impl->extractCigar(gm_reverse_strand.get(), _impl->node_map_));
        }
        nodes_passed = make_node_list(gm_reverse_strand.get());
    }
    else
    {
        read.set_graph_pos(gm_fwd_strand->position);
        read.set_graph_alignment_score(gm_fwd_strand->score);
        read.set_is_graph_alignment_unique(fwd_strand_is_unique);
        read.set_graph_mapq(fwd_strand_is_unique ? 60 : 0);
        if (alignment_flags & AF_CIGAR)
        {
            read.set_graph_cigar(_impl->extractCigar(gm_fwd_strand.get(), _impl->node_map_));
        }
        nodes_passed = make_node_list(gm_fwd_strand.get());
    }
    LOG()->trace("GraphAligner::alignRead read.cigar={}", read.graph_cigar());
}
void GraphAligner :: alignRead ( Read & read , unsigned int alignment_flags ) const { GraphAlignerImpl :: p_gssw_graph_mapping gm_fwd_strand { nullptr , safe_gssw_graph_mapping_destroy } ; GraphAlignerImpl :: p_gssw_graph_mapping gm_reverse_strand { nullptr , safe_gssw_graph_mapping_destroy } ; GraphAlignerImpl :: p_gssw_graph_mapping rgm_fwd_strand { nullptr , safe_gssw_graph_mapping_destroy } ; GraphAlignerImpl :: p_gssw_graph_mapping rgm_reverse_strand { nullptr , safe_gssw_graph_mapping_destroy } ; const std :: string rev_cmp_bases = reverseComplement ( read . bases ( ) ) ; const bool fwd_strand_multi_align = _impl -> alignString ( _impl -> graph_ . get ( ) , _impl -> nodes_ , _impl -> node_map_ , read . bases ( ) , gm_fwd_strand ) ; const bool reverse_strand_multi_align = ( alignment_flags & AF_BOTH_STRANDS ) != 0u ? _impl -> alignString ( _impl -> graph_ . get ( ) , _impl -> nodes_ , _impl -> node_map_ , rev_cmp_bases , gm_reverse_strand ) : false ; bool rfwd_strand_multi_align = false ; bool rreverse_strand_multi_align = false ; if ( ( alignment_flags & AF_REVERSE_GRAPH ) != 0u ) { string bases_rev = read . bases ( ) ; std :: reverse ( bases_rev . begin ( ) , bases_rev . end ( ) ) ; rfwd_strand_multi_align = _impl -> alignString ( _impl -> graph_reversed_ . get ( ) , _impl -> nodes_reversed_ , _impl -> node_map_reversed_ , bases_rev , rgm_fwd_strand ) ; rreverse_strand_multi_align = ( alignment_flags & AF_BOTH_STRANDS ) != 0u ? _impl -> alignString ( _impl -> graph_reversed_ . get ( ) , _impl -> nodes_reversed_ , _impl -> node_map_reversed_ , reverseComplement ( bases_rev ) , rgm_reverse_strand ) : false ; } bool fwd_strand_is_unique = ( ! fwd_strand_multi_align ) && ( ! rfwd_strand_multi_align ) ; bool reverse_strand_is_unique = ( ! reverse_strand_multi_align ) && ( ! rreverse_strand_multi_align ) ; bool return_reverse = false ; if ( ( ! fwd_strand_is_unique ) && reverse_strand_is_unique && gm_reverse_strand ) { return_reverse = true ; } else if ( fwd_strand_is_unique && ( ! reverse_strand_is_unique ) ) { return_reverse = false ; } else if ( gm_reverse_strand ) { return_reverse = gm_fwd_strand -> score < gm_reverse_strand -> score ; } const bool resulting_strand_is_reverse = read . is_reverse_strand ( ) != return_reverse ; read . set_is_graph_reverse_strand ( resulting_strand_is_reverse ) ; auto make_node_list = [ this ] ( gssw_graph_mapping * gm ) -> std :: list < int > { std :: list < int > result ; auto const & g = gm -> cigar ; auto nc = g . elements ; for ( uint32_t i = 0 ; i < g . length ; ++ i , ++ nc ) { result . push_back ( _impl -> node_map_ [ nc -> node -> id ] ) ; } return result ; } ; std :: list < int > nodes_passed ; if ( return_reverse ) { read . set_bases ( rev_cmp_bases ) ; std :: string rev_quals = read . quals ( ) ; std :: reverse ( rev_quals . begin ( ) , rev_quals . end ( ) ) ; read . set_quals ( rev_quals ) ; read . set_graph_pos ( gm_reverse_strand -> position ) ; read . set_graph_alignment_score ( gm_reverse_strand -> score ) ; read . set_is_graph_alignment_unique ( reverse_strand_is_unique ) ; read . set_graph_mapq ( reverse_strand_is_unique ? 60 : 0 ) ; if ( alignment_flags & AF_CIGAR ) { read . set_graph_cigar ( _impl -> extractCigar ( gm_reverse_strand . get ( ) , _impl -> node_map_ ) ) ; } nodes_passed = make_node_list ( gm_reverse_strand . get ( ) ) ; } else { read . set_graph_pos ( gm_fwd_strand -> position ) ; read . set_graph_alignment_score ( gm_fwd_strand -> score ) ; read . set_is_graph_alignment_unique ( fwd_strand_is_unique ) ; read . set_graph_mapq ( fwd_strand_is_unique ? 60 : 0 ) ; if ( alignment_flags & AF_CIGAR ) { read . set_graph_cigar ( _impl -> extractCigar ( gm_fwd_strand . get ( ) , _impl -> node_map_ ) ) ; } nodes_passed = make_node_list ( gm_fwd_strand . get ( ) ) ; } LOG ( ) -> trace ( " " , read . graph_cigar ( ) ) ; }
Align a read to the graph and update the graph_* fields.
Align a read to the graph and update the graph_ * fields .
// prefer unique alignment. if not unique, or if both unique, prefer alignment with better score