Alphabets in BioC++

Table of Contents

• Introduction
• The nucleotide alphabets
  • Construction and assignment of alphabet symbols
  • The rank of an alphabet symbol
  • The char representation of an alphabet symbol
  • Obtaining the alphabet size
  • containers over alphabets
  • Comparability
  • Example
• Other alphabets
  • The amino acid alphabet
  • Quality alphabets
  • Gap alphabet

Learning Objective:

In this tutorial we look at alphabets and you will learn how to work with nucleotides and amino acids in BioC++. We guide you through the most important properties of BioC++'s alphabets and show you the different implemented types. After completion, you will be able to use the alphabets inside of STL containers and to compare alphabet values.

Difficulty	Easy
Duration	45 min
Prerequisite tutorials	Quick Setup (using CMake)
Recommended reading	None

The links on this page mostly point straight into the API documentation which you should use as a reference. The code examples and assignments are designed to provide some practical experience with our interface as well as a code basis for your own program development.

Introduction

An alphabet is the set of symbols of which a biological sequence – or in general a text – is composed. BioC++ implements specific and optimised alphabets not only for sequences of RNA, DNA and protein components, but also for quality, secondary structure and gap annotation as well as combinations of the aforementioned.

Task

Read the section Detailed Description of the API reference page for Alphabet. This is a detailed introduction to the alphabet module and demonstrates its main advantages.

The nucleotide alphabets

Nucleotides are the components of (Deoxy)Ribonucleic acid (DNA/RNA) and contain one of the nucleobases Adenine (A), Cytosine (C), Guanine (G), Thymine (T, only DNA) and Uracil (U, only RNA). In BioC++ the alphabets bio::alphabet::dna4 and bio::alphabet::rna4 contain exactly the four respective nucleotides. The trailed number in the alphabets' name represents the number of entities the alphabet holds – we denote this number as alphabet size. For instance, the alphabet bio::alphabet::dna5 represents five entities as it contains the additional symbol 'N' to refer to an unknown nucleotide.

Construction and assignment of alphabet symbols

Let's look at some example code which demonstrates how objects of the bio::alphabet::dna4 alphabet are assigned from characters.

#include <bio/alphabet/all.hpp> // for working with alphabets directly
 
using namespace bio::alphabet::literals;
 
int main ()
{
    // Two objects of bio::alphabet::dna4 alphabet constructed with a char literal.
    bio::alphabet::dna4 ade = 'A'_dna4;
    bio::alphabet::dna4 gua = 'G'_dna4;
 
    // Two additional objects assigned explicitly from char or rank.
    bio::alphabet::dna4 cyt, thy;
    cyt.assign_char('C');
    thy.assign_rank(3);
 
    // Further code here...

    return 0;
}

We have shown three solutions for assigning variables of alphabet type.

1. Construction by character literal, i.e. appending the operator _dna4 to the respective char symbol.
   This is the handiest way as it can be also used as a temporary object.
2. Assignment by char via the global function bio::alphabet::assign_char.
   This is useful if the assignment target already exists, e.g. in a sequence vector.
3. Assignment by rank via the global function bio::alphabet::assign_rank.
   May be used when the rank is known.

The rank of an alphabet symbol

The rank of a symbol is a number in range [0..alphabet_size) where each number is paired with an alphabet symbol by a bijective function. In BioC++ the rank is always determined by the lexicographical order of the underlying characters. For instance, in bio::alphabet::dna4 the bijection is
'A'_dna4 ⟼ 0
'C'_dna4 ⟼ 1
'G'_dna4 ⟼ 2
'T'_dna4 ⟼ 3.

BioC++ provides the function bio::alphabet::to_rank for converting a symbol to its rank value as demonstrated in the following code example. Note that the data type of the rank is usually the smallest possible unsigned type that is required for storing the values of the alphabet.

    // Get the rank type of the alphabet (here uint8_t).
    using rank_type = bio::alphabet::rank_t<bio::alphabet::dna4>;
 
    // Retrieve the numerical representation (rank) of the objects.
    rank_type rank_a = ade.to_rank();   // => 0
    rank_type rank_g = gua.to_rank();   // => 2

The char representation of an alphabet symbol

Our alphabets also have a character representation because it is more intuitive to work with them than using the rank. Each alphabet symbol is represented by its respective character whenever possible (A ⟼ 'A'). Analogously to the rank, BioC++ provides the function bio::alphabet::to_char for converting a symbol to its character representation.

    // Get the character type of the alphabet (here char).
    using char_type = bio::alphabet::char_t<bio::alphabet::dna4>;
 
    // Retrieve the character representation.
    char_type char_a = ade.to_char();   // => 'A'
    char_type char_g = gua.to_char();   // => 'G'

Above you have seen that you can assign an alphabet symbol from a character with bio::alphabet::from_char. In contrast to the rank interface, this assignment is not a bijection because the whole spectrum of available chars is mapped to values inside the alphabet. For instance, assigning to bio::alphabet::dna4 from any character other than C, G or T results in the value 'A'_dna4 and assigning from any character except A, C, G or T to bio::alphabet::dna5 results in the value 'N'_dna5. You can avoid the implicit conversion by using bio::alphabet::assign_char_strict which throws bio::alphabet::invalid_char_assignment on invalid characters.

    // Assign from character with value check.
    bio::alphabet::assign_char_strictly_to('C', cyt);
 
    // bio::alphabet::assign_char_strictly_to('X', thy); // would throw bio::alphabet::invalid_char_assignment

You can test the validity of a character by calling bio::alphabet::char_is_valid_for. It returns true if the character is valid and false otherwise.

Obtaining the alphabet size

You can retrieve the alphabet size by accessing the class member variable alphabet_size which is implemented in most bio::alphabet::alphabet instances.

    // Get the alphabet size as class member of the alphabet.
    uint8_t const size1 = bio::alphabet::dna4::alphabet_size;        // => 4

containers over alphabets

In BioC++ you can use the STL containers to model e.g. sequences, sets or mappings with our alphabets. The following example shows some exemplary contexts for their use. For sequences we recommend the std::vector with one of BioC++'s alphabet types. Please note how easily a sequence can be created via the string literal.

    using namespace bio::alphabet::literals;
 
    // Examples of different container types with SeqAn's alphabets.
    std::vector<bio::alphabet::dna5> dna_sequence{"GATTANAG"_dna5};
    std::pair<bio::alphabet::gapped<bio::alphabet::dna4>, bio::alphabet::gapped<bio::alphabet::dna4>> alignment_column{bio::alphabet::gap{}, thy};
    std::set<bio::alphabet::dna4> pyrimidines{'C'_dna4, 'T'_dna4};

Comparability

All alphabets in BioC++ are regular and comparable. This means that you can use the <, <=, > and >= operators to compare the values based on the rank. For instance, this is useful for sorting a text over the alphabet. Regular implies that the equality and inequality of two values can be tested with == and !=.

    // Equality and comparison of bio::alphabet::dna4 symbols.
    bool eq  = (cyt == 'C'_dna4); // true
    bool neq = (thy != 'C'_dna4); // true
    bool geq = (cyt >= 'C'_dna4); // true
    bool gt  = (thy >  'C'_dna4); // true
    bool seq = (cyt <= 'C'_dna4); // true
    bool st  = (ade <  'C'_dna4); // true
 
    // Sort a vector of symbols.
    std::vector<bio::alphabet::dna4> some_nucl{"GTA"_dna4};
    std::sort(some_nucl.begin(), some_nucl.end()); // some_nucl: "AGT"

Example

To wrap up this section on the nucleotide alphabet, the following assignment lets you practise the use of a BioC++ alphabet and its related functions. It will also show you a handy advantage of using a vector over an alphabet instead of using std::string: The rank representation can be used straight as an array index (opposed to using a map with logarithmic access times, for example).

Assignment 1: GC content of a sequence

An important property of DNA and RNA sequences is the GC content, which is the percentage of nucleobases that are either Guanine or Cytosine. Given the nucleotide counts , , , the GC content is calculated as

Write a program that

1. reads a sequence as command line argument into a vector of bio::alphabet::dna5,
2. counts the number of occurrences for each nucleotide in an array of size alphabet size and
3. calculates the GC content.

The bio::alphabet::dna5 type ensures that invalid characters in the input sequence are converted to 'N'. Note that these characters should not influence the GC content.

Pass the sequences CATTACAG (3/8 = 37.5%) and ANNAGAT (1/5 = 20%) to your program and check if your results are correct.

Solution

#include <array>     // std::array
#include <string>    // std::string
#include <vector>    // std::vector
 
#include <bio/alphabet/all.hpp>
#include <bio/alphabet/fmt.hpp>
#include <bio/ranges/views/all.hpp>       // optional: use views to convert the input string to a dna5 sequence
 
using namespace bio::alphabet::literals;
 
/* call with single argument: SEQUENCE  */
int main (int argc, char * argv[])
{
    if (argc != 2)
    {
        // add error handling here
        return 0;
    }
    std::string input = argv[1];
 
    // Option1: Convert the input to a dna5 sequence by copy-conversion
    std::vector<bio::alphabet::dna5> sequence;
    sequence.resize(input.size());
    for (size_t i = 0; i < input.size(); ++i)
        bio::alphabet::assign_char_to(input[i], sequence[i]);
 
    // Option1: use views for the conversion. Views will be introduced in the next chapter.
    //auto sequence = input | bio::views::char_to<bio::alphabet::dna5>;
 
    // Initialise an array with count values for dna5 symbols.
    std::array<size_t, bio::alphabet::dna5::alphabet_size> count{}; // default initialised with zeroes
 
    // Increase the symbol count according to the sequence.
    for (bio::alphabet::dna5 symbol : sequence)
        ++count[symbol.to_rank()];
 
    // Calculate the GC content: (#G + #C) / (#A + #T + #G + #C).
    size_t gc = count['C'_dna5.to_rank()] + count['G'_dna5.to_rank()];
    size_t atgc = input.size() - count['N'_dna5.to_rank()];
    float gc_content = 1.0f * gc / atgc;
 
    fmt::print("The GC content of {} is {}%\n", sequence, 100 * gc_content);
 
    return 0;
}

Other alphabets

Until now, we have focused on alphabets for nucleotides to introduce the properties of BioC++'s alphabet on a specific example. BioC++ implements, however, many more alphabets. In this section, we want to give you an overview of the other existing alphabets.

The amino acid alphabet

Proteins consist of one or more long chains of amino acids. The so-called primary structure of a protein is expressed as sequences over an amino acid alphabet. The bio::alphabet::aa27 alphabet contains the standard one-letter code of the 20 canonical amino acids as well as the two proteinogenic amino acids, a termination symbol and some wildcard characters. For details read the Aminoacid page.

Quality alphabets

The alphabets for quality are sequence annotations since they describe additional properties of the respective sequence. The values are produced by sequencing machines and represent the probability that a nucleobase was recorded incorrectly. The characters are most commonly found in FASTQ files. See Quality for details.

You can build an Alphabet Tuple Composite with a nucleotide and quality alphabet that stores both information together. For the use cases just described we offer a pre-defined composite (bio::alphabet::qualified). See our API documentation for a detailed description of each.

Gap alphabet

The bio::alphabet::gap alphabet is the smallest alphabet in BioC++, consisting only of the gap character. It is most often used in an Alphabet Variant with a nucleotide or amino acid alphabet to represent gapped sequences, e.g. in alignments. To create a gapped alphabet simply use bio::alphabet::gapped<> with the alphabet type you want to refine.

    // Assign a gap symbol to a gapped RNA alphabet.
    bio::alphabet::gapped<bio::alphabet::rna5> sym = bio::alphabet::gap{};                         // => -
 
    using namespace bio::alphabet::literals;
    // Each bio::alphabet::rna5 symbol is still valid.
    sym = 'U'_rna5;                                   // => U
 
    // The alphabet size is six (AUGCN-).
    uint8_t const size2 = bio::alphabet::gapped<bio::alphabet::rna5>::alphabet_size;   // => 6