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
78
79
80
81
82
83
84
85
86
|
// The comment lines containing ![fragment-line] are there for the
// documentation system. You can ignore them when reading this file.
//![full]
//![one]
#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()
{
// 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;
// 1. 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
//![one]
String<unsigned> alignPosList;
resize(alignPosList, 4);
alignPosList[0] = 7;
alignPosList[1] = 100;
alignPosList[2] = 172;
alignPosList[3] = 272;
// 2. 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';
//![two]
// 3. Print alignments of the reads with chr1 (or bsChr1) sequence using the function printAlign
// and the positions in alignPosList.
// To do that, you have to create a copy of the fragment in chr1 (bsChr1) that is aligned to the read.
std::cout << " \n Print alignment: " << std::endl;
for (unsigned i = 0; i < length(readList); ++i)
{
// Begin position beginPosition of a given alignment between the read and the genome
//![two]
unsigned beginPosition = alignPosList[i];
//![three]
// Genome fragment
DnaString genomeFragment;
//![three]
// We have to create a copy of the corresponding fragment of the genome, where the read aligns to
for (unsigned j = 0; j < length(readList[i]); ++j)
appendValue(genomeFragment, chr1[beginPosition + j]);
//![four]
// Call of our function to print the simple alignment
printAlign(genomeFragment, readList[i]);
}
return 0;
}
//![four]
//![full]
|
