ex_alignment_ioΒΆ

Unit test which reads alignment in Edinburgh format or calculates a global sequence alignment for two predefined sequences. It saves output alignment in Edinburgh format.

USAGE:

./ex_alignment_io [alignment]

EXAMPLE:

./ex_alignment_io example.edinb

Keywords:

Categories:

  • core/data/io/alignment_io.hh

Input files:

Output files:

Program source:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#include <iostream>

#include <core/data/io/alignment_io.hh>
#include <core/BioShellEnvironment.hh>
#include <core/alignment/NWAligner.hh>
#include <core/alignment/on_alignment_computations.hh>
#include <core/alignment/PairwiseAlignment.hh>
#include <core/alignment/PairwiseSequenceAlignment.hh>
#include <core/alignment/scoring/SimilarityMatrix.hh>
#include <core/alignment/scoring/SimilarityMatrixScore.hh>
#include <utils/options/OptionParser.hh>
#include <utils/exit.hh>

std::string program_info = R"(

Unit test which reads alignment in Edinburgh format or calculates a global sequence alignment
for two predefined sequences. It saves output alignment in Edinburgh format.

USAGE:
  ./ex_alignment_io [alignment]

EXAMPLE:
  ./ex_alignment_io example.edinb

)";

/** @brief Read alignment in Edinburgh format or calculate a new one from given sequences; write Edinburgh.
 *
 * CATEGORIES: core/data/io/alignment_io.hh
 * KEYWORDS:   sequence alignment
 */
int main(const int argc, const char* argv[]) {

  if ((argc > 1) && utils::options::call_for_help(argv[1]))
    utils::exit_OK_with_message(program_info);

  using namespace core::alignment;
  using namespace core::alignment::scoring;

  if (argc > 1) { // If there was an input alignment file given, read it!
    std::vector<PairwiseSequenceAlignment_SP> alignments;
    auto ali = core::data::io::read_edinburgh(argv[1], alignments);
    for (const PairwiseSequenceAlignment_SP & ali : alignments)
      core::data::io::write_edinburgh(*ali, std::cout, 80);
  } else { // otherwise, align the two sequences defined below
    // ---------- The two sequences that will be aligned
    std::string query = "MKGWLFLVIAIVGEVIATSALKSSEGFTKLAPSAVVIIGYGIAFYFLSLVLKSIPVGVAYAVWSGLGVVIITAIAWLLHGQKLDAWGFVGMGLI";
    std::string tmplt = "MIYLYLLCAIFAEVVATSLLKSTEGFTRLWPTVGCLVGYGIAFALLALSISHGMQTDVAYALWSAIGTAAIVLVAVLFLGSPISVMKVVGVGLI";
    // ---------- Gap penalties
    short int open = -10;
    short int extend = -2;
    // ---------- load BLOSUM matrix from bioshell's library; the directory must be defined as a shell variable
    const std::shared_ptr<SimilarityMatrix<short int>> b62_matrix = SimilarityMatrix<short int>::from_ncbi_file("alignments/BLOSUM62");
    const SimilarityMatrixScore<short int> b62_score(query, tmplt, *b62_matrix);
    NWAligner<short int, SimilarityMatrixScore<short int>> global(std::max(query.length(), tmplt.length()));
    // ---------- Compute and backtrace the alignment
    global.align(open, extend, b62_score);
    const PairwiseAlignment_SP ali = global.backtrace();
    // ---------- Convert the abstract alignment to a pairwise sequence alignment object
   core::data::sequence::Sequence query_seq("Q7B1Y7_SALEN",query);
   core::data::sequence::Sequence tmplt_seq("MMR_MYCTU",tmplt);
   core::data::io::write_edinburgh(query_seq, *ali, tmplt_seq, std::cout, 80);
  }
}
../_images/file_icon.png