ex_PairwiseSequenceAlignmentΒΆ

Simple example showing how to work with PairwiseSequenceAlignment data structure, e.g. how to create such an object and hot to print it in different formats.

USAGE:

./ex_PairwiseSequenceAlignment

Keywords:

Categories:

  • core::alignment::PairwiseAlignment; core::alignment::PairwiseSequenceAlignment

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
65
66
67
68
69
70
71
72
73
74
75
76
77
#include <iostream>

#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/data/io/alignment_io.hh>
#include <core/data/sequence/Sequence.hh>
#include <utils/options/OptionParser.hh>

std::string program_info = R"(

Simple example showing how to work with PairwiseSequenceAlignment data structure, e.g. how to
create such an object and hot to print it in different formats.

USAGE:
./ex_PairwiseSequenceAlignment

)";

std::string Q52825_1 = "IDVLLGADDGSLAFVPSEFSISPGEKIVFKNNAGFPHNIVFDEDSIPSGVDASKISMSEEDLLNAKGETFEVALSNKGEYSFYCSPHQGAGMVGKVTV";
std::string P80401_2 = "VQMLNKGTDGAMVFEPGFLKIAPGDTVTFIPTDKS-HNVETFKGLIPDGV---------PDFKSKPNEQYQVKFDIPGAYVLKCTPHVGMGMVALIQV";

/** @brief Simple example showing how to work with PairwiseSequenceAlignment data structure.
 *
 * CATEGORIES: core::alignment::PairwiseAlignment; core::alignment::PairwiseSequenceAlignment;
 * 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;
  using namespace core::data::sequence; // for core::data::sequence::Sequence

  // ---------- Test for global alignment ----------
  // --- Alignment defined as path : '-' and '-' mean a gap in a template and in a query sequence, respectively; '*' is a match
  PairwiseAlignment_SP ali = std::make_shared<PairwiseAlignment>(0, 0, 0, "--****-**|**--");

  // ---------- The two sequences that will be aligned
  Sequence_SP query = std::make_shared<Sequence>("query", "ITFTALILLAVAV", 1);
  Sequence_SP tmplt = std::make_shared<Sequence>("tmplt", "FTALLLAAV", 1);

  PairwiseSequenceAlignment seq_ali(ali, query, tmplt);

  // --- Show alignment as a path
  std::cout << "Alignment path:\n" << ali->to_path() << "\n\n";

  // ---------- Print the alignment in Edinburgh format
  core::index2 identity = sum_identical(seq_ali);
  core::index2 n_gaps = seq_ali.alignment->length() - seq_ali.alignment->n_aligned();
  std::cout << "# score: " << seq_ali.alignment_score() << " length: " << seq_ali.alignment->length()
            << " n_identical: " << identity << " n_gaps: " << n_gaps << "\n";
  core::data::io::write_edinburgh(seq_ali, std::cout, 80);

  // ---------- Test for local alignment ----------
  // --- Alignment defined as path : '-' and '-' mean a gap in a template and in a query sequence, respectively; '*' is a match
  PairwiseAlignment_SP loc_ali = std::make_shared<PairwiseAlignment>(2, 0, 0.0, "****-**|**");
  PairwiseSequenceAlignment loc_seq_ali(loc_ali, query, tmplt);

  // ---------- Print the alignment in Edinburgh format
  identity = sum_identical(loc_seq_ali);
  n_gaps = loc_seq_ali.alignment->length() - loc_seq_ali.alignment->n_aligned();
  std::cout << "# score: " << loc_seq_ali.alignment_score() << " length: " << loc_seq_ali.alignment->length()
            << " n_identical: " << identity << " n_gaps: " << n_gaps << "\n";
  core::data::io::write_edinburgh(loc_seq_ali, std::cout, 80);


  PairwiseSequenceAlignment loc_seq_ali2("Q52825_1", Q52825_1, 0, "P80401_2", P80401_2, 0, 0.0);
  loc_seq_ali2.template_sequence->first_pos(28);
  loc_seq_ali2.query_sequence->first_pos(1);
  core::data::io::write_edinburgh(loc_seq_ali2, std::cout, 80);
}
../_images/file_icon.png