alignment_from_cigar.hpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <vector>
 
#include <seqan3/alignment/aligned_sequence/aligned_sequence_concept.hpp>
#include <seqan3/alignment/decorator/gap_decorator.hpp>
#include <seqan3/alphabet/cigar/cigar.hpp>
#include <seqan3/io/sam_file/detail/cigar.hpp>
#include <seqan3/utility/views/slice.hpp>
 
namespace seqan3
{
 
template <typename reference_type, typename sequence_type>
inline auto alignment_from_cigar(std::vector<cigar> const & cigar_vector,
                                 reference_type const & reference,
                                 uint32_t const zero_based_reference_start_position,
                                 sequence_type const & query)
{
    if (cigar_vector.empty())
        throw std::logic_error{"An empty CIGAR is not a valid alignment representation."};
 
    // compute the length of the aligned region in the reference sequence
    // -------------------------------------------------------------------------
    // this requires a first stream over the cigar vector.
    uint32_t reference_length{0};
    uint32_t query_length{0};
 
    for (auto [cigar_count, cigar_operation] : cigar_vector)
    {
        if (('M'_cigar_operation == cigar_operation) || ('='_cigar_operation == cigar_operation)
            || ('X'_cigar_operation == cigar_operation))
        {
            reference_length += cigar_count;
            query_length += cigar_count;
        }
        else if ('D'_cigar_operation == cigar_operation)
        {
            reference_length += cigar_count;
        }
        else if ('I'_cigar_operation == cigar_operation)
        {
            query_length += cigar_count;
        }
    }
 
    if (static_cast<size_t>(zero_based_reference_start_position + reference_length) > std::ranges::size(reference))
        throw std::logic_error{"The CIGAR string indicates a reference length of at least "
                               + std::to_string(zero_based_reference_start_position + reference_length)
                               + ", but the supplied reference sequence is only of size"
                               + std::to_string(std::ranges::size(reference)) + "."};
 
    // Get soft clipping at the start and end of the CIGAR string
    // -------------------------------------------------------------------------
    uint32_t soft_clipping_start{0};
    uint32_t soft_clipping_end{0};
 
    // Checks whether the given index in the cigar vector is a soft clip.
    auto soft_clipping_at = [&cigar_vector](size_t const index)
    {
        return cigar_vector[index] == 'S'_cigar_operation;
    };
    // Checks whether the given index in the cigar vector is a hard clip.
    auto hard_clipping_at = [&](size_t const index)
    {
        return cigar_vector[index] == 'H'_cigar_operation;
    };
    // Checks whether the given cigar vector has at least min_size many elements.
    auto vector_size_at_least = [&](size_t const min_size)
    {
        return cigar_vector.size() >= min_size;
    };
    // Returns the cigar count of the ith cigar element in the given cigar vector.
    auto cigar_count_at = [&](size_t const index)
    {
        return get<0>(cigar_vector[index]);
    };
 
    // check for soft clipping at the first two positions
    // cigar is non-empty, checked at the very beginning.
    if (soft_clipping_at(0))
        soft_clipping_start = cigar_count_at(0);
    else if (vector_size_at_least(2) && hard_clipping_at(0) && soft_clipping_at(1))
        soft_clipping_start = cigar_count_at(1);
 
    // Check for soft clipping at the last two positions to validate the CIGAR string.
    // Even if the two following arithmetics overflow, they are protected by the corresponding if expressions below.
    auto const last_index = cigar_vector.size() - 1;
    auto const second_last_index = last_index - 1;
 
    if (vector_size_at_least(2) && soft_clipping_at(last_index))
        soft_clipping_end = cigar_count_at(last_index);
    else if (vector_size_at_least(3) && hard_clipping_at(last_index) && soft_clipping_at(second_last_index))
        soft_clipping_end = cigar_count_at(second_last_index);
 
    if (soft_clipping_start + query_length + soft_clipping_end != std::ranges::size(query))
        throw std::logic_error{"The CIGAR string indicates a query/read sequence length of "
                               + std::to_string(soft_clipping_start + query_length + soft_clipping_end)
                               + ", but the supplied query/read sequence is of size"
                               + std::to_string(std::ranges::size(query)) + "."};
 
    // Assign the sequence to the alignment (a tuple of 2 gap decorators)
    // -------------------------------------------------------------------------
    using gapped_reference_type = gap_decorator<decltype(reference | views::slice(0, 0))>;
    using gapped_sequence_type = gap_decorator<decltype(query | views::slice(0, 0))>;
    using alignment_type = std::tuple<gapped_reference_type, gapped_sequence_type>;
 
    alignment_type alignment{};
 
    assign_unaligned(get<0>(alignment),
                     reference
                         | views::slice(zero_based_reference_start_position,
                                        zero_based_reference_start_position + reference_length));
    // query_length already accounts for soft clipping at begin and end
    assign_unaligned(get<1>(alignment), query | views::slice(soft_clipping_start, soft_clipping_start + query_length));
 
    // Insert gaps into the alignment based on the cigar vector
    // -------------------------------------------------------------------------
    using std::get;
    auto current_ref_pos = std::ranges::begin(get<0>(alignment));
    auto current_read_pos = std::ranges::begin(get<1>(alignment));
 
    for (auto [cigar_count, cigar_operation] : cigar_vector)
    {
        if (('M'_cigar_operation == cigar_operation) || ('='_cigar_operation == cigar_operation)
            || ('X'_cigar_operation == cigar_operation))
        {
            std::ranges::advance(current_ref_pos, cigar_count);
            std::ranges::advance(current_read_pos, cigar_count);
        }
        else if (('D'_cigar_operation == cigar_operation) || ('N'_cigar_operation == cigar_operation))
        {
            // insert gaps into query
            current_read_pos = get<1>(alignment).insert_gap(current_read_pos, cigar_count);
            ++current_read_pos;
            std::ranges::advance(current_ref_pos, cigar_count);
        }
        else if (('I'_cigar_operation == cigar_operation)) // Insert gaps into ref
        {
            current_ref_pos = get<0>(alignment).insert_gap(current_ref_pos, cigar_count);
            ++current_ref_pos;
            std::ranges::advance(current_read_pos, cigar_count);
        }
        else if (('P'_cigar_operation == cigar_operation)) // skip padding
        {
            current_ref_pos = get<0>(alignment).insert_gap(current_ref_pos, cigar_count);
            ++current_ref_pos;
 
            current_read_pos = get<1>(alignment).insert_gap(current_read_pos, cigar_count);
            ++current_read_pos;
        }
        // S and H are ignored as they are handled by cropping the sequence
    }
 
    return alignment;
}
 
template <typename reference_type, typename sequence_type>
inline auto alignment_from_cigar(std::string const & cigar_string,
                                 reference_type const & reference,
                                 uint32_t const zero_based_reference_start_position,
                                 sequence_type const & query)
{
    alignment_from_cigar(detail::parse_cigar(cigar_string), reference, zero_based_reference_start_position, query);
}
 
} // namespace seqan3