cigar_from_alignment.hpp

Go to the documentation of this file.

// -----------------------------------------------------------------------------------------------------
// Copyright (c) 2006-2022, Knut Reinert & Freie Universität Berlin
// Copyright (c) 2016-2022, Knut Reinert & MPI für molekulare Genetik
// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License
// shipped with this file and also available at: https://github.com/seqan/seqan3/blob/master/LICENSE.md
// -----------------------------------------------------------------------------------------------------
 
#pragma once
 
#include <vector>
 
#include <seqan3/alphabet/cigar/cigar.hpp>
#include <seqan3/alphabet/gap/gap.hpp>
#include <seqan3/utility/views/zip.hpp>
 
namespace seqan3
{
 
struct cigar_clipped_bases
{
    uint32_t hard_front{}; 
    uint32_t hard_back{};  
    uint32_t soft_front{}; 
    uint32_t soft_back{};  
};
 
template <typename alignment_type>
inline auto cigar_from_alignment(alignment_type const & alignment,
                                 cigar_clipped_bases const & clipped_bases = {},
                                 bool const extended_cigar = false)
{
    static_assert((tuple_like<std::remove_cvref_t<alignment_type>>
                   && std::tuple_size_v<std::remove_cvref_t<alignment_type>> == 2),
                  "The alignment must be a std::pair or std::tuple of size 2.");
 
    static_assert(
        (std::equality_comparable_with<gap, std::ranges::range_reference_t<decltype(std::get<0>(alignment))>>
         && std::equality_comparable_with<gap, std::ranges::range_reference_t<decltype(std::get<1>(alignment))>>),
        "The alignment must contain two ranges whose value_type is comparable to seqan3::gap.");
 
    /*brief: Compares two seqan3::aligned_sequence values and returns their cigar operation.
    * param  reference_char      The aligned character of the reference to compare.
    * param  query_char          The aligned character of the query to compare.
    * param  extended_cigar      Whether to print the extended cigar alphabet or not. See cigar operation.
    * returns A seqan3::cigar::operation representing the alignment operation between the two values.
    *
    * The resulting CIGAR operation is based on the query character (`query_char`).
    *
    * ### Example:
    *
    * The following alignment column shows the reference char ('C') on top and a
    * gap for the query char at the bottom.
    * ```
    * ... C ...
    *     |
    * ... - ...
    * ```
    * In this case, `to_cigar_op` will return
    * 'D' since the query char is "deleted".
    *
    * The next alignment column shows the reference char ('C') on top and a
    * query char ('G') at the bottom.
    * ```
    * ... C ...
    *     |
    * ... G ...
    * ```
    * In this case, `to_cigar_op` will return
    * 'M', for the basic CIGAR the two bases are aligned, while
    * in the extended CIGAR alphabet (`extended_cigar` = `true`) the function
    * will return an 'X' since the bases are aligned but are not
    * equal.
    * \sa seqan3::aligned_sequence
    */
    auto to_cigar_op = [extended_cigar](auto const reference_char, auto const query_char)
    {
        // note that N is not considered because it is equivalent to D but has a special meaning:
        // SAM spec: "For mRNA-to-genome alignment, an N operation represents an intron. For other types of alignments,
        //            the interpretation of N is not defined."
        // as we cannot know the meaning, the user has to change D -> N themself
        constexpr std::array<char, 6> operators{'M', 'D', 'I', 'P', 'X', '='}; // contains the possible cigar operators.
 
        // no gaps               -> 00 -> 0
        // gap only in query     -> 01 -> 1
        // gap only in reference -> 10 -> 2
        // both gaps             -> 11 -> 3
        uint8_t key = (static_cast<uint8_t>(reference_char == gap{}) << 1) | static_cast<uint8_t>(query_char == gap{});
 
        // mismatch -> 100 -> 4
        // match    -> 101 -> 5
        if (extended_cigar && (key == 0)) // in extended format, refine the substitution operator to match/mismatch.
            key |= ((1 << 2) | static_cast<uint8_t>(query_char == reference_char)); // maps to [4, 5].
 
        return assign_char_to(operators[key], cigar::operation{});
    };
 
    using std::get;
 
    auto & ref_seq = get<0>(alignment);
    auto & query_seq = get<1>(alignment);
 
    if (ref_seq.size() != query_seq.size())
        throw std::logic_error{"The aligned sequences (including gaps) must have the same length."};
 
    if (std::ranges::empty(ref_seq)) // only check ref_seq because query_seq was checked to have to same size
        throw std::logic_error{"The aligned sequences may not be empty."};
 
    std::vector<cigar> result{};
 
    // Add (H)ard-clipping at the start of the query
    if (clipped_bases.hard_front)
        result.emplace_back(clipped_bases.hard_front, 'H'_cigar_operation);
 
    // Add (S)oft-clipping at the start of the query
    if (clipped_bases.soft_front)
        result.emplace_back(clipped_bases.soft_front, 'S'_cigar_operation);
 
    // Create cigar string from alignment
    // -------------------------------------------------------------------------
    // initialize first operation and count value:
    cigar::operation operation{to_cigar_op(ref_seq[0], query_seq[0])};
    uint32_t count{0};
 
    // go through alignment columns
    for (auto && [reference_char, query_char] : views::zip(ref_seq, query_seq))
    {
        cigar::operation next_op = to_cigar_op(reference_char, query_char);
 
        if (operation == next_op)
        {
            ++count;
        }
        else
        {
            result.emplace_back(count, operation);
            operation = next_op;
            count = 1;
        }
    }
 
    // append last cigar element
    result.emplace_back(count, operation);
 
    // Add (S)oft-clipping at the end of the query
    if (clipped_bases.soft_back)
        result.emplace_back(clipped_bases.soft_back, 'S'_cigar_operation);
 
    // Add (H)ard-clipping at the end of the query
    if (clipped_bases.hard_back)
        result.emplace_back(clipped_bases.hard_back, 'H'_cigar_operation);
 
    return result;
}
 
} // namespace seqan3