Loading...
Searching...
No Matches
CIGAR Conversion
Sequence Alignment

The CIGAR Conversion submodule contains utility functions to convert a CIGAR to an alignment or vice versa. More...

Collaboration diagram for CIGAR Conversion:

## Classes

struct | seqan3::cigar_clipped_bases
| Helper struct to specialise soft and hard clipping when using seqan3::cigar_from_alignment. More...

## Functions

template<typename reference_type, typename sequence_type>
auto | seqan3::alignment_from_cigar (std::string const &cigar_string, reference_type const &reference, uint32_t const zero_based_reference_start_position, sequence_type const &query)
template<typename reference_type, typename sequence_type>
auto | seqan3::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)
| Construct an alignment from a CIGAR string and the corresponding sequences.
template<typename alignment_type>
auto | seqan3::cigar_from_alignment (alignment_type const &alignment, cigar_clipped_bases const &clipped_bases={}, bool const extended_cigar=false)
| Creates a CIGAR string (SAM format) given a seqan3::detail::pairwise_alignment represented by two seqan3::aligned_sequences.

## Detailed Description

The CIGAR Conversion submodule contains utility functions to convert a CIGAR to an alignment or vice versa.

# Quick background on the CIGAR string

The CIGAR string is a compact representation of an aligned read against a reference and was introduced by the SAM format. The SAM format stores the result of mapping short/long read sequences from a sequencing experiment (e.g., Illumina/Nanopore) against a reference (e.g., hg38).

# Creating an alignment from a CIGAR String

To create an alignment from a CIGAR String, you can do the following:

// SPDX-FileCopyrightText: 2006-2025 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2025 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
#include <seqan3/alignment/cigar_conversion/alignment_from_cigar.hpp>
#include <seqan3/core/debug_stream.hpp>
#include <seqan3/io/sam_file/all.hpp>
using namespace seqan3::literals;
auto sam_file_raw = R"(@HD VN:1.6
@SQ SN:ref LN:34
read1 41 ref 1 61 1S1M1D1M1I ref 10 300 ACGT !##$ AS:i:2 NM:i:7
read2 42 ref 2 62 1H7M1D1M1S2H ref 10 300 AGGCTGNAG !##$&'()* xy:B:S,3,4,5
read3 43 ref 3 63 1S1M1P1M1I1M1I1D1M1S ref 10 300 GGAGTATA !!*+,-./
)";
int main()
{
// The reference sequence might be read from a different file.
seqan3::dna5_vector reference = "ACTGATCGAGAGGATCTAGAGGAGATCGTAGGAC"_dna5;
seqan3::sam_file_input fin{std::istringstream{sam_file_raw}, seqan3::format_sam{}};
// You will probably read it from a file, e.g., like this:
// seqan3::sam_file_input fin{"test.sam"};
for (auto && rec : fin)
{
auto alignment =
alignment_from_cigar(rec.cigar_sequence(), reference, rec.reference_position().value(), rec.sequence());
seqan3::debug_stream << alignment << '\n';
}
// prints:
// (ACT-,C-GT)
// (CTGATCGAG,AGGCTGN-A)
// (T-G-A-TC,G-AGTA-T)
}
alignment_from_cigar.hpp
Provides the function seqan3::alignment_from_cigar.
std::istringstream
seqan3::format_sam
The SAM format (tag).
Definition format_sam.hpp:105
seqan3::sam_file_input
A class for reading SAM files, both SAM and its binary representation BAM are supported.
Definition sam_file/input.hpp:236
debug_stream.hpp
Provides seqan3::debug_stream and related types.
seqan3::alignment_from_cigar
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)
Construct an alignment from a CIGAR string and the corresponding sequences.
Definition alignment_from_cigar.hpp:81
seqan3::debug_stream
debug_stream_type debug_stream
A global instance of seqan3::debug_stream_type.
Definition debug_stream.hpp:38
seqan3::field::alignment
@ alignment
The (pairwise) alignment stored in an object that models seqan3::detail::pairwise_alignment.
Definition record.hpp:78
all.hpp
Meta-header for the IO / SAM File submodule .
seqan3::literals
The SeqAn namespace for literals.

# Creating a CIGAR string from an alignment

To create a CIGAR String from an alignment, you can do the following:

// SPDX-FileCopyrightText: 2006-2025 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2025 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
#include <seqan3/alignment/aligned_sequence/aligned_sequence_concept.hpp>
#include <seqan3/alignment/cigar_conversion/cigar_from_alignment.hpp>
#include <seqan3/alignment/decorator/gap_decorator.hpp>
#include <seqan3/alphabet/gap/gapped.hpp>
#include <seqan3/alphabet/nucleotide/dna5.hpp>
#include <seqan3/core/debug_stream.hpp>
#include <seqan3/utility/views/slice.hpp>
int main()
{
using namespace seqan3::literals;
seqan3::dna5_vector reference = "ATGGCGTAGAGCTTCCCCCCCCCCCCCCCCC"_dna5;
seqan3::dna5_vector read = "ATGCCCCGTTGCTT"_dna5; // length 14
// Let's say, we want to ignore the last 2 bases of the query because the quality is low.
// We thus only align the first 12 bases, the last two will be soft-clipped bases in the CIGAR string.
seqan3::gap_decorator aligned_reference{reference | seqan3::views::slice(0, 12)};
seqan3::gap_decorator aligned_query{read | seqan3::views::slice(0, 12)};
// insert gaps
seqan3::insert_gap(aligned_reference, aligned_reference.begin() + 4, 2);
seqan3::insert_gap(aligned_query, aligned_query.begin() + 11, 2);
auto cigar_sequence =
seqan3::cigar_from_alignment(std::tie(aligned_reference, aligned_query),
{.hard_front = 1, .hard_back = 0, .soft_front = 0, .soft_back = 2});
seqan3::debug_stream << cigar_sequence << '\n'; // prints [1H,4M,2I,5M,2D,1M,2S]
}
aligned_sequence_concept.hpp
Includes the aligned_sequence and the related insert_gap and erase_gap functions to enable stl contai...
cigar_from_alignment.hpp
Provides the function seqan3::cigar_from_alignment and a helper struct seqan3::cigar_clipped_bases.
seqan3::gap_decorator
A gap decorator allows the annotation of sequences with gap symbols while leaving the underlying sequ...
Definition gap_decorator.hpp:78
seqan3::gap_decorator::begin
const_iterator begin() const noexcept
Returns an iterator to the first element of the container.
Definition gap_decorator.hpp:349
dna5.hpp
Provides seqan3::dna5, container aliases and string literals.
gap_decorator.hpp
Provides seqan3::gap_decorator.
gapped.hpp
Provides seqan3::gapped.
seqan3::insert_gap
std::ranges::iterator_t< range_type > insert_gap(range_type &rng, std::ranges::iterator_t< range_type > const pos_it, typename range_type::size_type const size=1)
An implementation of seqan3::writable_aligned_sequence::insert_gap for ranges with the corresponding ...
Definition aligned_sequence_concept.hpp:409
seqan3::cigar_from_alignment
auto cigar_from_alignment(alignment_type const &alignment, cigar_clipped_bases const &clipped_bases={}, bool const extended_cigar=false)
Creates a CIGAR string (SAM format) given a seqan3::detail::pairwise_alignment represented by two seq...
Definition cigar_from_alignment.hpp:111
seqan3::views::slice
constexpr auto slice
A view adaptor that returns a half-open interval on the underlying range.
Definition slice.hpp:137
slice.hpp
Provides seqan3::views::slice.
std::tie
T tie(T... args)

See also
  seqan3::alignment_from_cigar
  seqan3::cigar_from_alignment

## Function Documentation

## alignment_from_cigar() [1/2]

template<typename reference_type, typename sequence_type>

auto seqan3::alignment_from_cigar | ( | std::string const & | cigar_string,
|  | reference_type const & | reference,
|  | uint32_t const | zero_based_reference_start_position,
|  | sequence_type const & | query )

inline

This is an overloaded member function, provided for convenience. It differs from the above function only in what argument(s) it accepts.

## alignment_from_cigar() [2/2]

template<typename reference_type, typename sequence_type>

auto seqan3::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 )

inline

Construct an alignment from a CIGAR string and the corresponding sequences.

Template Parameters

  reference_type | The type of the reference sequence for a SAM record.
  sequence_type | The type of the read sequence for a SAM record.

Parameters

  [in] | cigar_vector | The CIGAR information to convert to an alignment.
  [in] | reference | The reference sequence to which the query was aligned to, the alignment being represented by cigar_vector.
  [in] | zero_based_reference_start_position | The start position of the alignment in the reference sequence. The position is zero-based. When using our seqan3::sam_file_input, note that the seqan3::sam_file_input::record_type::reference_position() is always zero-based.
  [in] | query | The query or read sequence of the alignment represented by cigar_vector.

Returns
  An alignment represented by a std::tuple of size 2 holding 2 seqan3::gap_decorators. At position 0 is the aligned reference sequence and at position 1 the aligned read sequence.

### Quick background on the CIGAR string

The CIGAR string is a compact representation of an aligned read against a reference and was introduced by the SAM format. The SAM format stores the result of mapping short/long read sequences from a sequencing experiment (e.g., Illumina/Nanopore) against a reference (e.g., hg38).

### Conversion to an alignment

You can reconstruct a full alignment from a CIGAR string, if you have the respective sequences at hand:

// SPDX-FileCopyrightText: 2006-2025 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2025 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
#include <seqan3/alignment/cigar_conversion/alignment_from_cigar.hpp>
#include <seqan3/alphabet/cigar/cigar.hpp>
#include <seqan3/alphabet/nucleotide/dna5.hpp>
#include <seqan3/core/debug_stream.hpp>
using namespace seqan3::literals;
int main()
{
// CIGAR string = 2M1D2M
std::vector<seqan3::cigar> cigar_vector{{2, 'M'_cigar_operation},
{1, 'D'_cigar_operation},
{2, 'M'_cigar_operation}};
uint32_t reference_start_position{0}; // The read is aligned at the start of the reference.
seqan3::dna5_vector reference = "ACTGATCGAGAGGATCTAGAGGAGATCGTAGGAC"_dna5;
seqan3::dna5_vector query = "ACGA"_dna5;
auto alignment = alignment_from_cigar(cigar_vector, reference, reference_start_position, query);
seqan3::debug_stream << alignment << '\n'; // prints (ACTGA,AC-GA)
}
cigar.hpp
Provides the seqan3::cigar alphabet.
std::vector

### Quick explanation of the alignment representation

In seqan3, an alignment is represented by a std::tuple of size 2 that holds 2 seqan3::aligned_sequences.

The data structure that we use most often to model seqan3::aligned_sequence is the seqan3::gap_decorator. It is a lightweight data structure that only holds a view on the sequence (no copy is made) and on top can hold seqan3::gaps.

In the above example, the read sequence ACGT is aligned to the reference with one gap: AC-GA where - represents a gap. In the CIGAR string, the gap in the query/read is represented by 1D.

The full alignment consist of two aligned sequences (read and reference). In the above example, the alignment

position 01234
reference ACTGA
read AC-GA

is represented by a tuple of the aligned reference at the first position (std::get<0>) and the aligned read at the second position (std::get<1>): (ACTGA,AC-GA).

### IO Example

A more realistic example is extracting the information directly from a SAM file:

// SPDX-FileCopyrightText: 2006-2025 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2025 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
#include <seqan3/alignment/cigar_conversion/alignment_from_cigar.hpp>
#include <seqan3/core/debug_stream.hpp>
#include <seqan3/io/sam_file/all.hpp>
using namespace seqan3::literals;
auto sam_file_raw = R"(@HD VN:1.6
@SQ SN:ref LN:34
read1 41 ref 1 61 1S1M1D1M1I ref 10 300 ACGT !##$ AS:i:2 NM:i:7
read2 42 ref 2 62 1H7M1D1M1S2H ref 10 300 AGGCTGNAG !##$&'()* xy:B:S,3,4,5
read3 43 ref 3 63 1S1M1P1M1I1M1I1D1M1S ref 10 300 GGAGTATA !!*+,-./
)";
int main()
{
// The reference sequence might be read from a different file.
seqan3::dna5_vector reference = "ACTGATCGAGAGGATCTAGAGGAGATCGTAGGAC"_dna5;
seqan3::sam_file_input fin{std::istringstream{sam_file_raw}, seqan3::format_sam{}};
// You will probably read it from a file, e.g., like this:
// seqan3::sam_file_input fin{"test.sam"};
for (auto && rec : fin)
{
auto alignment =
alignment_from_cigar(rec.cigar_sequence(), reference, rec.reference_position().value(), rec.sequence());
seqan3::debug_stream << alignment << '\n';
}
// prints:
// (ACT-,C-GT)
// (CTGATCGAG,AGGCTGN-A)
// (T-G-A-TC,G-AGTA-T)
}

See also
  seqan3::sam_file_input
  seqan3::cigar

## cigar_from_alignment()

template<typename alignment_type>

auto seqan3::cigar_from_alignment | ( | alignment_type const & | alignment,
|  | cigar_clipped_bases const & | clipped_bases = {},
|  | bool const | extended_cigar = false )

inline

Creates a CIGAR string (SAM format) given a seqan3::detail::pairwise_alignment represented by two seqan3::aligned_sequences.

Template Parameters

  alignment_type | Must model seqan3::detail::pairwise_alignment.

Parameters

  alignment | The alignment, represented by a pair of aligned sequences, to be transformed into CIGAR vector based on the second (query) sequence.
  clipped_bases | Provides information on whether the query sequence was cropped (hard clipping) before the alignment or whether part of the query sequence does not take part (soft clipping) in the alignment.
  extended_cigar | Whether to print the extended CIGAR alphabet or not. See seqan3::cigar::operation.

Returns
  A std::vector<seqan3::cigar> representing the alignment.

Note
  The resulting cigar_vector is based on the query sequence, which is the second sequence in the Sequence Alignment pair.

### Example:

Given the following alignment, reference sequence on top and the query sequence at the bottom:

ATGG--CGTAGAGCTT
|||X |||X| |||
ATGCCCCGTTG--CTT

In this case, the function seqan3::cigar_from_alignment will return the following CIGAR vector:

[('M',4),('I',2),('M',5),('D',2),('M',3)]

The extended CIGAR string would look like this:

[('=',3)('X',1)('I',2)('=',3)('X',1)('=',1)('D',2)('=',3)]
// SPDX-FileCopyrightText: 2006-2025 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2025 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
#include <seqan3/alignment/aligned_sequence/aligned_sequence_concept.hpp>
#include <seqan3/alignment/cigar_conversion/cigar_from_alignment.hpp>
#include <seqan3/alignment/decorator/gap_decorator.hpp>
#include <seqan3/alphabet/gap/gapped.hpp>
#include <seqan3/alphabet/nucleotide/dna5.hpp>
#include <seqan3/core/debug_stream.hpp>
#include <seqan3/utility/views/slice.hpp>
int main()
{
using namespace seqan3::literals;
seqan3::dna5_vector reference = "ATGGCGTAGAGCTTCCCCCCCCCCCCCCCCC"_dna5;
seqan3::dna5_vector read = "ATGCCCCGTTGCTT"_dna5; // length 14
// Align the full query against the first 14 bases of the reference.
seqan3::gap_decorator aligned_reference{reference | seqan3::views::slice(0, 14)};
seqan3::gap_decorator aligned_read{read};
// Insert gaps to represent the alignment:
seqan3::insert_gap(aligned_read, aligned_read.begin() + 11, 2);
seqan3::insert_gap(aligned_reference, aligned_reference.begin() + 4, 2);
seqan3::debug_stream << aligned_reference << '\n' << aligned_read << '\n';
// prints:
// ATGG--CGTAGAGCTT
// ATGCCCCGTTG--CTT
auto cigar_sequence = seqan3::cigar_from_alignment(std::tie(aligned_reference, aligned_read));
seqan3::debug_stream << cigar_sequence << '\n'; // prints [4M,2I,5M,2D,3M]
}

### Soft and Hard clipping

The terms soft and hard clipping were introduced by the SAM specifications. A SAM file is only storing a semi alignment represented by the CIGAR string. The semi alignment of a query sequence is most often the result of a read mapping.

#### Hard clipping

Before aligning a query or read to the reference, some tools crop the query sequence because the quality is bad at one end (e.g., Illumina reads tend to display a bad sequence quality towards the end of the read).

To inform the user of a SAM file that query sequences were altered, hard-clipping information is appended to the CIGAR string. E.g. 100M50H indicates that a read of former length 150 has been cropped at the right end by 50 bases. The sequence in the SAM file will thus only be 100 bases long.

#### Soft clipping

In contrast to hard clipping, soft clipping indicates that the read was cropped and the respective bases do not participate in the alignment, but they are still part of the reported sequence. E.g., 100M50S indicates that a read of length 150 has been aligned without the rightmost 50 bases. The sequence in the SAM file will still be 150 bases long.

#### Adding soft and hard clipping with seqan3::cigar_from_alignment

You can add the respective clipping information by passing an instance of seqan3::cigar_clipped_bases:

// SPDX-FileCopyrightText: 2006-2025 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2025 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
#include <seqan3/alignment/aligned_sequence/aligned_sequence_concept.hpp>
#include <seqan3/alignment/cigar_conversion/cigar_from_alignment.hpp>
#include <seqan3/alignment/decorator/gap_decorator.hpp>
#include <seqan3/alphabet/gap/gapped.hpp>
#include <seqan3/alphabet/nucleotide/dna5.hpp>
#include <seqan3/core/debug_stream.hpp>
#include <seqan3/utility/views/slice.hpp>
int main()
{
using namespace seqan3::literals;
seqan3::dna5_vector reference = "ATGGCGTAGAGCTTCCCCCCCCCCCCCCCCC"_dna5;
seqan3::dna5_vector read = "ATGCCCCGTTGCTT"_dna5; // length 14
// Let's say, we want to ignore the last 2 bases of the query because the quality is low.
// We thus only align the first 12 bases, the last two will be soft-clipped bases in the CIGAR string.
seqan3::gap_decorator aligned_reference{reference | seqan3::views::slice(0, 12)};
seqan3::gap_decorator aligned_query{read | seqan3::views::slice(0, 12)};
// insert gaps
seqan3::insert_gap(aligned_reference, aligned_reference.begin() + 4, 2);
seqan3::insert_gap(aligned_query, aligned_query.begin() + 11, 2);
auto cigar_sequence =
seqan3::cigar_from_alignment(std::tie(aligned_reference, aligned_query),
{.hard_front = 1, .hard_back = 0, .soft_front = 0, .soft_back = 2});
seqan3::debug_stream << cigar_sequence << '\n'; // prints [1H,4M,2I,5M,2D,1M,2S]
}

See also
  seqan3::cigar_clipped_bases
  seqan3::sam_file_output
  seqan3::cigar
  seqan3::aligned_sequence

Hide me