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
|
#include <iostream>
#include <seqan/sequence.h>
#include <seqan/stream.h>
using namespace seqan2;
// Function to print simple alignment between two sequences with the same length
template <typename TText1, typename TText2>
void printAlign(TText1 const & genomeFragment, TText2 const & read)
{
std::cout << "Alignment " << std::endl;
std::cout << " genome : " << genomeFragment << std::endl;
std::cout << " read : " << read << std::endl;
}
int main(int, char const **)
{
// Build reads and genomes
DnaString chr1 = "TATAATATTGCTATCGCGATATCGCTAGCTAGCTACGGATTATGCGCTCTGCGATATATCGCGCTAGATGTGCAGCTCGATCGAATGCACGTGTGTGCGATCGATTAGCGTCGATCATCGATCTATATTAGCGCGCGGTATCGGACGATCATATTAGCGGTCTAGCATTTAG";
// Build List containing all reads
typedef String<DnaString> TDnaList;
TDnaList readList;
resize(readList, 4);
readList[0] = "TTGCTATCGCGATATCGCTAGCTAGCTACGGATTATGCGCTCTGCGATATATCGCGCT";
readList[1] = "TCGATTAGCGTCGATCATCGATCTATATTAGCGCGCGGTATCGGACGATCATATTAGCGGTCTAGCATT";
readList[2] = "AGCCTGCGTACGTTGCAGTGCGTGCGTAGACTGTTGCAAGCCGGGGGTTCATGTGCGCTGAAGCACACATGCACA";
readList[3] = "CGTGCACTGCTGACGTCGTGGTTGTCACATCGTCGTGCGTGCGTACTGCTGCTGACA";
// Append a second chromosome sequence fragment to chr1
DnaString chr2 = "AGCCTGCGTACGTTGCAGTGCGTGCGTAGACTGTTGCAAGCCGGGGGTTCATGTGCGCTGAAGCACACATGCACACGTCTCTGTGTTCCGACGTGTGTCACGTGCACTGCTGACGTCGTGGTTGTCACATCGTCGTGCGTGCGTACTGCTGCTGACACATGCTGCTG";
append(chr1, chr2);
// Print readlist
std::cout << " \n Read list: " << std::endl;
for (unsigned i = 0; i < length(readList); ++i)
std::cout << readList[i] << std::endl;
// Assume we have mapped the 4 reads to chr1 (and chr2) and now have the mapping start positions (no gaps).
// Store the start position in a String alignPosList: 7, 100, 172, 272
String<unsigned> alignPosList;
resize(alignPosList, 4);
alignPosList[0] = 7;
alignPosList[1] = 100;
alignPosList[2] = 172;
alignPosList[3] = 272;
// Optional
// Bisulfite conversion
// Assume chr1 is beeing bisulfate treated: Copy chr1 to a new genome bsChr1 and exchange every 'C' with a 'T'
DnaString bsChr1;
assign(bsChr1, chr1);
for (unsigned i = 0; i < length(bsChr1); ++i)
if (bsChr1[i] == 'C')
bsChr1[i] = 'T';
// Print alignments using Segment: Do the same as above, but instead of using a for loop to build the fragment,
// use the Segment class to build an infix of bsChr1.
// Note: Because printAlign uses templates, we don't have to change the function even though the type of
// genomeFragment is different.
std::cout << " \n Print alignment using Segment: " << std::endl;
for (unsigned i = 0; i < length(readList); ++i)
{
// Begin and end position of a given alignment between the read and the genome
unsigned beginPosition = alignPosList[i];
unsigned endPosition = beginPosition + length(readList[i]);
// Build infix
Infix<DnaString>::Type genomeFragment = infix(chr1, beginPosition, endPosition);
// Call of our function to print the simple alignment
printAlign(genomeFragment, readList[i]);
}
return 0;
}
|
