Global Alignments

Example Program

Global Alignments

Computing an optimal global alignment between two sequences.

A tutorial about global alignments.

1	//#include <iostream>
2
3	#include <seqan/basic.h>
4	#include <seqan/file.h>
5	#include <seqan/align.h>
6
7	int main()
8	{
9	using namespace seqan;
10	typedef Value<Gaps<Dna5String, ArrayGaps> >::Type TValue;
11	using namespace seqan;
12

Two DNA sequences that shall be aligned.

13	typedef String<Dna> TSequence;
14	TSequence seq1 = "atcgaatgcgga";
15	TSequence seq2 = "actcgttgca";

Scoring objects are used to define a scoring scheme. In this case, affine gap costs with match = 0, mismatch = -1, gapextend = -1 and gapopen = -2.

16	Score<int> scoringScheme(0, -1, -1, -2);

Example 1: We use Align to align the two sequences. Since we do not specify an algorithm tag when we call globalAlignment, a suitable algorithm (Gotoh) is automatically choosen.

17	Align<TSequence, ArrayGaps> align;
18	resize(rows(align), 2);
19	assignSource(row(align, 0), seq1);
20	assignSource(row(align, 1), seq2);
21
22	int score = globalAlignment(align, scoringScheme);
23	std::cout << "Score = " << score << std::endl;
24	std::cout << align << std::endl;

Example 2: We now choose explicitely the algorithm MyersHirschberg. Since this algorithm always works on Levenshtein distance, we do not need score.

25	score = globalAlignment(align, MyersHirschberg());
26	std::cout << "Score = " << score << std::endl;
27	std::cout << align << std::endl;

Example 3: We now do the same as in case 1, but now we use an Alignment Graph for storing the alignment. Here we use Gotoh's algorithm.

28	typedef StringSet<TSequence, Dependent<> > TStringSet;
29	typedef Graph<Alignment<TStringSet, void> > TAlignmentGraph;
30
31	TStringSet string_set;
32	appendValue(string_set, seq1);
33	appendValue(string_set, seq2);
34	TAlignmentGraph alignment_graph(string_set);
35
36	score = globalAlignment(alignment_graph, scoringScheme, Gotoh());
37	std::cout << "Score = " << score << std::endl;
38	std::cout << alignment_graph << std::endl;
39	return 0;
40	}

See

Alignment Graph, Align, Score, globalAlignment

SeqAn - Sequence Analysis Library - www.seqan.de