Type

Review

Objective

Compute two global alignments between the DNA sequences"AAATGACGGATTG"."AGTCGGATCTACTG" using the Gotoh algorithm[Got82], implementing the Affine Gap model, with the following scoring parameters:match=4,mismatch=-2,gapOpen=-4 andgapExtend=-2.Store the alignments in two Align objects and print them together with the scores.

Hints

The Gotoh algorithm uses the Affine Gap function. In SeqAn you can switch between Linear, Affine and Dynamic gap functions by customizing your scoring scheme with one of the three tagsLinearGaps(),AffineGaps() orDynamicGaps() and relative penalty valuesgapOpen andgapExtend. When a single gap value is provided the Linear Gap model is selected as default while the Affine Gap model is chosen as standard when two different gap costs are set. If the Dynamic Gap model[UPA+14] is required the relative tag must be supplied.Have a look on theScoring Schemes section if you are not sure about the correct ordering.

Solution

First we have to define the body of our program.This includes the definition of the library headers that we want to use.In this case it is theiostream from the STL and the<seqan/align.h>header file defining all algorithms and data structures we want to use.After we added the namespace and opened themain body we define our types we want to use in this function.We use anString with theDna alphabet, since we know that we work with DNA sequences.The second type is ourAlign object storing the alignment later on.

#include<iostream>#include<seqan/align.h>usingnamespaceseqan2;intmain(){typedefString<Dna>TSequence;// sequence typetypedefAlign<TSequence,ArrayGaps>TAlign;// align type

In the next step we initialize our objects.This includes the both input sequencesseq1 andseq2 andalign.We resize the underlying set ofalign that manages the separateGaps data structures.Finally, we assign the input sequences as sources to the corresponding row ofalign.

TSequenceseq1="AAATGACGGATTG";TSequenceseq2="AGTCGGATCTACTG";TAlignalign;resize(rows(align),2);assignSource(row(align,0),seq1);assignSource(row(align,1),seq2);

Now we compute the alignment using a scoring scheme with affine gap costs.The first parameter corresponds to thematch value, the second to themismatch value, the third to thegapextend value and the last one to thegapopen value.We store the computed score of the best alignment in the equally named variablescore.In the end we print the score and the alignment using print methods provided by theiostream module of the STL.

intscore=globalAlignment(align,Score<int,Simple>(4,-2,-2,-4),AffineGaps());std::cout<<"Score: "<<score<<std::endl;std::cout<<align<<std::endl;return0;}

Congratulation!You have computed an alignment using affine gap costs.Here the result of the program:

Score: 16      0     .    :    .        AAATGACGGAT----TG        |   | |||||    ||        A---GTCGGATCTACTG

Type

Review

Objective

Compute a semi-global alignment between the sequencesAAATGACGGATTG andTGGGA using the costs 1, -1, -1 for match, mismatch and gap, respectively.Use an AlignmentGraph to store the alignment.Print the score and the resulting alignment to the standard output.

Hint

A semi-global alignment is a special form of an overlap alignment often used when aligning short sequences against a long sequence.Here we only allow free end-gaps at the beginning and the end of the shorter sequence.

Solution

First we have to define the body of our program.This includes the definition of the library headers that we want to use.In this case we include theiostream header from the STL and the<seqan/align.h> header, which defines all algorithms and data structures we want to use.After we added the namespace and opened themain function body we define our types we want to use in this function.We use anString with theDna alphabet, since we know that we work with DNA sequences.We use an additionalStringSet to store the input sequences.In this scenario we use anAlignmentGraph to store the alignment.Remember, that the AlignmentGraph uses anDependentStringSet to map the vertices to the correct input sequences.

#include<iostream>#include<seqan/align.h>usingnamespaceseqan2;intmain(){typedefString<Dna>TSequence;// sequence typetypedefStringSet<TSequence>TStringSet;// container for stringstypedefStringSet<TSequence,Dependent<>>TDepStringSet;// dependent string settypedefGraph<Alignment<TDepStringSet>>TAlignGraph;// alignment graph

In the next step we initialize our input StringSetstrings and pass it as argument to the constructor of the AlignmentGraphalignG.

TSequenceseq1="AAATGACGGATTG";TSequenceseq2="TGGGA";TStringSetsequences;appendValue(sequences,seq1);appendValue(sequences,seq2);TAlignGraphalignG(sequences);

Now we compute the alignment using the Levenshtein distance and a AlignConfig object to set the correct free end-gaps.In this example we put the shorter sequence on the vertical axis of our alignment matrix.Hence, we have to use free end-gaps in the first and last row, which corresponds to the first and the last parameter in the AlignConfig object.If you add the shorter sequence at first tostrings, then you simply have to flip thebool values of the AlignConfig object.

intscore=globalAlignment(alignG,Score<int,Simple>(1,-1,-1),AlignConfig<true,false,false,true>());std::cout<<"Score: "<<score<<std::endl;std::cout<<alignG<<std::endl;return0;}

Here the result of the program.

Score: 3Alignment matrix:      0     .    :        AAATGACGGATTG           ||  |||        ---TG--GGA---

Type

Application

Objective

Write a program that computes a global alignment between theRna sequencesAAGUGACUUAUUG andAGUCGGAUCUACUG using the Myers-Hirschberg variant. You should use the Align data structure to store the alignment.Print the score and the alignment. Additionally, output for each row of the Align object the view positions of the gaps.

Hint

You can use an iterator to iterate over a row.Use the metafunctionRow to get the type of the row used by the Align object.Use the functionisGap to check whether the current value of the iterator is a gap or not.The gaps are already in the view space.

Solution

As usual, first the necessary includes and typedefs.Our sequence type isString<Rna>.TAlign andTRow are defined as in the previous example program.The typeIterator<TRow>::Type will be used to iterate over the rows of the alignment.

#include<iostream>#include<seqan/align.h>usingnamespaceseqan2;intmain(){typedefString<Rna>TSequence;typedefAlign<TSequence,ArrayGaps>TAlign;typedefRow<TAlign>::TypeTRow;typedefIterator<TRow>::TypeTRowIterator;

In the next step we initialize our Align objectalign with the corresponding source files.

TSequenceseq1="AAGUGACUUAUUG";TSequenceseq2="AGUCGGAUCUACUG";TAlignalign;resize(rows(align),2);assignSource(row(align,0),seq1);assignSource(row(align,1),seq2);

Now we compute the alignment using Myers-Hirschberg algorithm by specifying the correct tag at the end of the function.

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

Finally, we iterate over both gap structures and print the view positions of the gaps within the sequences.

unsignedaliLength=_max(length(row(align,0)),length(row(align,1)));for(unsignedi=0;i<length(rows(align));++i){TRowIteratorit=iter(row(align,i),0);TRowIteratoritEnd=iter(row(align,i),aliLength);unsignedpos=0;std::cout<<"Row "<<i<<" contains gaps at positions: ";std::cout<<std::endl;while(it!=itEnd){if(isGap(it))std::cout<<pos<<std::endl;++it;++pos;}}return0;}

The output of the program is as follows.

Score: -6      0     .    :    .        AAGU--GA-CUUAUUG        | ||  || || | ||        A-GUCGGAUCU-ACUGRow 0 contains gaps at positions:458Row 1 contains gaps at positions:111

Type

Review

Objective

Write a program which computes the 3 best local alignments of the twoAminoAcid sequences “PNCFDAKQRTASRPL” and “CFDKQKNNRTATRDTA” using the following scoring parameters:match=3,mismatch=-2,gapopen=-5,gapextension=-1.

Hint

Use an extra variable to enumerate the k best alignments.

Solution

The usual includes.

#include<iostream>#include<seqan/align.h>usingnamespaceseqan2;intmain(){

The initialization of theAlign object.

Align<String<AminoAcid>>ali;resize(rows(ali),2);assignSource(row(ali,0),"PNCFDAKQRTASRPL");assignSource(row(ali,1),"CFDKQKNNRTATRDTA");

Computing the three best alignments with the desired scoring parameters:

Score<int>sc(3,-2,-1,-5);unsignedcount=0;LocalAlignmentEnumerator<Score<int>,Unbanded>enumerator(sc);while(nextLocalAlignment(ali,enumerator)&&count<3){std::cout<<"Score = "<<getScore(enumerator)<<std::endl;std::cout<<ali;++count;}return0;}

The resulting output is as follows.

Score = 21      0     .    :        CFDAKQ---RTASR        ||| ||   ||| |        CFD-KQKNNRTATRScore = 8      0     .        AKQR-TA        |  | ||        AT-RDTAScore = 5      0        D-A        | |        DTA

Type

Transfer

Objective

Write an approximate pattern matching algorithm using alignment algorithms.Report the positions of all hits where the pattern matches the text with at most2 errors.Output the number of total edits used to match the pattern and print the corresponding cigar string of the alignment without leading and trailing gaps in the pattern.Text: “MISSISSIPPIANDMISSOURI” Pattern: “SISSI”

Hint

• The first step would be to verify at which positions in the text the pattern matches with at most 2 errors.

• Use theinfix function to return a subsequence of a string.

• A CIGAR string is a different representation of an alignment.It consists of a number followed by an operation.The number indicates how many operations are executed.Operations can beM for match or mismatch,I for insertion andD for deletion.Here is an example:

  ref:AC--GTCATTTr01:ACGTCTCA---Cigarofr01:2M2I4M3D

Solution (Step 1)

#include<iostream>#include<seqan/align.h>usingnamespaceseqan2;intmain(){typedefString<char>TSequence;TSequencetext="MISSISSIPPIANDMISSOURI";TSequencepattern="SISSI";String<int>locations;for(unsignedi=0;i<length(text)-length(pattern);++i){// Compute the MyersBitVector in current window of text.TSequencetmp=infix(text,i,i+length(pattern));// Report hits with at most 2 errors.if(globalAlignmentScore(tmp,pattern,MyersBitVector())>=-2){appendValue(locations,i);}}return0;}

Solution (Step 2)

#include<iostream>#include<seqan/align.h>usingnamespaceseqan2;intmain(){typedefString<char>TSequence;typedefGaps<TSequence,ArrayGaps>TGaps;typedefIterator<String<int>>::TypeTIterator;TSequencetext="MISSISSIPPIANDMISSOURI";TSequencepattern="SISSI";String<int>locations;for(unsignedi=0;i<length(text)-length(pattern);++i){// Compute the MyersBitVector in current window of text.TSequencetmp=infix(text,i,i+length(pattern));// Report hits with at most 2 errors.if(globalAlignmentScore(tmp,pattern,MyersBitVector())>=-2){appendValue(locations,i);}}TGapsgapsText;TGapsgapsPattern;assignSource(gapsPattern,pattern);std::cout<<"Text: "<<text<<"\tPattern: "<<pattern<<std::endl;for(TIteratorit=begin(locations);it!=end(locations);++it){// Clear previously computed gaps.clearGaps(gapsText);clearGaps(gapsPattern);// Only recompute the area within the current window over the text.TSequencetextInfix=infix(text,*it,*it+length(pattern));assignSource(gapsText,textInfix);// Use semi-global alignment since we do not want to track leading/trailing gaps in the pattern.// Restirct search space using a band allowing at most 2 errors in vertical/horizontal direction.intscore=globalAlignment(gapsText,gapsPattern,Score<int>(0,-1,-1),AlignConfig<true,false,false,true>(),-2,2);std::cout<<"Hit at position "<<*it<<"\ttotal edits: "<<abs(score)<<std::endl;}return0;}

Solution (Step 3)

#include<iostream>#include<seqan/align.h>usingnamespaceseqan2;intmain(){typedefString<char>TSequence;typedefGaps<TSequence,ArrayGaps>TGaps;typedefIterator<TGaps>::TypeTGapsIterator;typedefIterator<String<int>>::TypeTIterator;TSequencetext="MISSISSIPPIANDMISSOURI";TSequencepattern="SISSI";String<int>locations;for(unsignedi=0;i<length(text)-length(pattern);++i){// Compute the MyersBitVector in current window of text.TSequencetmp=infix(text,i,i+length(pattern));// Report hits with at most 2 errors.if(globalAlignmentScore(tmp,pattern,MyersBitVector())>=-2){appendValue(locations,i);}}TGapsgapsText;TGapsgapsPattern;assignSource(gapsPattern,pattern);std::cout<<"Text: "<<text<<"\tPattern: "<<pattern<<std::endl;for(TIteratorit=begin(locations);it!=end(locations);++it){// Clear previously computed gaps.clearGaps(gapsText);clearGaps(gapsPattern);// Only recompute the area within the current window over the text.TSequencetextInfix=infix(text,*it,*it+length(pattern));assignSource(gapsText,textInfix);// Use semi-global alignment since we do not want to track leading/trailing gaps in the pattern.// Restirct search space using a band allowing at most 2 errors in vertical/horizontal direction.intscore=globalAlignment(gapsText,gapsPattern,Score<int>(0,-1,-1),AlignConfig<true,false,false,true>(),-2,2);TGapsIteratoritGapsPattern=begin(gapsPattern);TGapsIteratoritGapsEnd=end(gapsPattern);// Remove trailing gaps in pattern.intcount=0;while(isGap(--itGapsEnd))++count;setClippedEndPosition(gapsPattern,length(gapsPattern)-count);// Remove leading gaps in pattern.if(isGap(itGapsPattern)){setClippedBeginPosition(gapsPattern,countGaps(itGapsPattern));setClippedBeginPosition(gapsText,countGaps(itGapsPattern));}std::cout<<"Hit at position "<<*it<<"\ttotal edits: "<<abs(score)<<std::endl;}return0;}

Solution (Step 4)

#include<iostream>#include<seqan/align.h>usingnamespaceseqan2;intmain(){typedefString<char>TSequence;typedefGaps<TSequence,ArrayGaps>TGaps;typedefIterator<TGaps>::TypeTGapsIterator;typedefIterator<String<int>>::TypeTIterator;TSequencetext="MISSISSIPPIANDMISSOURI";TSequencepattern="SISSI";String<int>locations;for(unsignedi=0;i<length(text)-length(pattern);++i){// Compute the MyersBitVector in current window of text.TSequencetmp=infix(text,i,i+length(pattern));// Report hits with at most 2 errors.if(globalAlignmentScore(tmp,pattern,MyersBitVector())>=-2){appendValue(locations,i);}}TGapsgapsText;TGapsgapsPattern;assignSource(gapsPattern,pattern);std::cout<<"Text: "<<text<<"\tPattern: "<<pattern<<std::endl;for(TIteratorit=begin(locations);it!=end(locations);++it){// Clear previously computed gaps.clearGaps(gapsText);clearGaps(gapsPattern);// Only recompute the area within the current window over the text.TSequencetextInfix=infix(text,*it,*it+length(pattern));assignSource(gapsText,textInfix);// Use semi-global alignment since we do not want to track leading/trailing gaps in the pattern.// Restirct search space using a band allowing at most 2 errors in vertical/horizontal direction.intscore=globalAlignment(gapsText,gapsPattern,Score<int>(0,-1,-1),AlignConfig<true,false,false,true>(),-2,2);TGapsIteratoritGapsPattern=begin(gapsPattern);TGapsIteratoritGapsEnd=end(gapsPattern);// Remove trailing gaps in pattern.intcount=0;while(isGap(--itGapsEnd))++count;setClippedEndPosition(gapsPattern,length(gapsPattern)-count);// Remove leading gaps in pattern.if(isGap(itGapsPattern)){setClippedBeginPosition(gapsPattern,countGaps(itGapsPattern));setClippedBeginPosition(gapsText,countGaps(itGapsPattern));}// Reinitilaize the iterators.TGapsIteratoritGapsText=begin(gapsText);itGapsPattern=begin(gapsPattern);itGapsEnd=end(gapsPattern);// Use a stringstream to construct the cigar string.std::stringstreamcigar;while(itGapsPattern!=itGapsEnd){// Count insertions.if(isGap(itGapsText)){intnumGaps=countGaps(itGapsText);cigar<<numGaps<<"I";itGapsText+=numGaps;itGapsPattern+=numGaps;continue;}++itGapsText;++itGapsPattern;}std::cout<<"Hit at position "<<*it<<"\ttotal edits: "<<abs(score)<<std::endl;}return0;}

Solution (Step 5)

#include<iostream>#include<seqan/align.h>usingnamespaceseqan2;intmain(){typedefString<char>TSequence;typedefGaps<TSequence,ArrayGaps>TGaps;typedefIterator<TGaps>::TypeTGapsIterator;typedefIterator<String<int>>::TypeTIterator;TSequencetext="MISSISSIPPIANDMISSOURI";TSequencepattern="SISSI";String<int>locations;for(unsignedi=0;i<length(text)-length(pattern);++i){// Compute the MyersBitVector in current window of text.TSequencetmp=infix(text,i,i+length(pattern));// Report hits with at most 2 errors.if(globalAlignmentScore(tmp,pattern,MyersBitVector())>=-2){appendValue(locations,i);}}TGapsgapsText;TGapsgapsPattern;assignSource(gapsPattern,pattern);std::cout<<"Text: "<<text<<"\tPattern: "<<pattern<<std::endl;for(TIteratorit=begin(locations);it!=end(locations);++it){// Clear previously computed gaps.clearGaps(gapsText);clearGaps(gapsPattern);// Only recompute the area within the current window over the text.TSequencetextInfix=infix(text,*it,*it+length(pattern));assignSource(gapsText,textInfix);// Use semi-global alignment since we do not want to track leading/trailing gaps in the pattern.// Restirct search space using a band allowing at most 2 errors in vertical/horizontal direction.intscore=globalAlignment(gapsText,gapsPattern,Score<int>(0,-1,-1),AlignConfig<true,false,false,true>(),-2,2);TGapsIteratoritGapsPattern=begin(gapsPattern);TGapsIteratoritGapsEnd=end(gapsPattern);// Remove trailing gaps in pattern.intcount=0;while(isGap(--itGapsEnd))++count;setClippedEndPosition(gapsPattern,length(gapsPattern)-count);// Remove leading gaps in pattern.if(isGap(itGapsPattern)){setClippedBeginPosition(gapsPattern,countGaps(itGapsPattern));setClippedBeginPosition(gapsText,countGaps(itGapsPattern));}// Reinitilaize the iterators.TGapsIteratoritGapsText=begin(gapsText);itGapsPattern=begin(gapsPattern);itGapsEnd=end(gapsPattern);// Use a stringstream to construct the cigar string.std::stringstreamcigar;while(itGapsPattern!=itGapsEnd){// Count insertions.if(isGap(itGapsText)){intnumGaps=countGaps(itGapsText);cigar<<numGaps<<"I";itGapsText+=numGaps;itGapsPattern+=numGaps;continue;}// Count deletions.if(isGap(itGapsPattern)){intnumGaps=countGaps(itGapsPattern);cigar<<numGaps<<"D";itGapsText+=numGaps;itGapsPattern+=numGaps;continue;}++itGapsText;++itGapsPattern;}// Output the hit position in the text, the total number of edits and the corresponding cigar string.std::cout<<"Hit at position  "<<*it<<"\ttotal edits: "<<abs(score)<<std::endl;}return0;}

Solution (Step 6)

#include<iostream>#include<seqan/align.h>usingnamespaceseqan2;intmain(){typedefString<char>TSequence;typedefGaps<TSequence,ArrayGaps>TGaps;typedefIterator<TGaps>::TypeTGapsIterator;typedefIterator<String<int>>::TypeTIterator;TSequencetext="MISSISSIPPIANDMISSOURI";TSequencepattern="SISSI";String<int>locations;for(unsignedi=0;i<length(text)-length(pattern);++i){// Compute the MyersBitVector in current window of text.TSequencetmp=infix(text,i,i+length(pattern));// Report hits with at most 2 errors.if(globalAlignmentScore(tmp,pattern,MyersBitVector())>=-2){appendValue(locations,i);}}TGapsgapsText;TGapsgapsPattern;assignSource(gapsPattern,pattern);std::cout<<"Text: "<<text<<"\tPattern: "<<pattern<<std::endl;for(TIteratorit=begin(locations);it!=end(locations);++it){// Clear previously computed gaps.clearGaps(gapsText);clearGaps(gapsPattern);// Only recompute the area within the current window over the text.TSequencetextInfix=infix(text,*it,*it+length(pattern));assignSource(gapsText,textInfix);// Use semi-global alignment since we do not want to track leading/trailing gaps in the pattern.// Restirct search space using a band allowing at most 2 errors in vertical/horizontal direction.intscore=globalAlignment(gapsText,gapsPattern,Score<int>(0,-1,-1),AlignConfig<true,false,false,true>(),-2,2);TGapsIteratoritGapsPattern=begin(gapsPattern);TGapsIteratoritGapsEnd=end(gapsPattern);// Remove trailing gaps in pattern.intcount=0;while(isGap(--itGapsEnd))++count;setClippedEndPosition(gapsPattern,length(gapsPattern)-count);// Remove leading gaps in pattern.if(isGap(itGapsPattern)){setClippedBeginPosition(gapsPattern,countGaps(itGapsPattern));setClippedBeginPosition(gapsText,countGaps(itGapsPattern));}// Reinitilaize the iterators.TGapsIteratoritGapsText=begin(gapsText);itGapsPattern=begin(gapsPattern);itGapsEnd=end(gapsPattern);// Use a stringstream to construct the cigar string.std::stringstreamcigar;intnumMatchAndMismatches=0;while(itGapsPattern!=itGapsEnd){// Count insertions.if(isGap(itGapsText)){intnumGaps=countGaps(itGapsText);cigar<<numGaps<<"I";itGapsText+=numGaps;itGapsPattern+=numGaps;continue;}// Count deletions.if(isGap(itGapsPattern)){intnumGaps=countGaps(itGapsPattern);cigar<<numGaps<<"D";itGapsText+=numGaps;itGapsPattern+=numGaps;continue;}// Count matches and  mismatches.while(itGapsPattern!=itGapsEnd){if(isGap(itGapsPattern)||isGap(itGapsText))break;++numMatchAndMismatches;++itGapsText;++itGapsPattern;}if(numMatchAndMismatches!=0)cigar<<numMatchAndMismatches<<"M";numMatchAndMismatches=0;}// Output the hit position in the text, the total number of edits and the corresponding cigar string.std::cout<<"Hit at position  "<<*it<<"\ttotal edits: "<<abs(score)<<"\tcigar: "<<cigar.str()<<std::endl;}return0;}

Complete Solution (and more explanations)

Write themain body of the program with type definition and initialization of the used data structures.

#include<iostream>#include<seqan/align.h>usingnamespaceseqan2;intmain(){typedefString<char>TSequence;typedefGaps<TSequence,ArrayGaps>TGaps;typedefIterator<TGaps>::TypeTGapsIterator;typedefIterator<String<int>>::TypeTIterator;TSequencetext="MISSISSIPPIANDMISSOURI";TSequencepattern="SISSI";

In the first part of the algorithm we implement an alignment based verification process to identify positions in thesubject sequence at which we can find our pattern with at most2 errors.We slide the5*5 alignment matrix position by position over thesubject sequence and use theMeyersBitVector to verify the hits.If the score is greater or equal than-2, then we have found a hit.We store the begin position of the hit inlocations.

String<int>locations;for(unsignedi=0;i<length(text)-length(pattern);++i){// Compute the MyersBitVector in current window of text.TSequencetmp=infix(text,i,i+length(pattern));// Report hits with at most 2 errors.if(globalAlignmentScore(tmp,pattern,MyersBitVector())>=-2){appendValue(locations,i);}}

In the second part of the algorithm we iterate over all reported locations.This time we compute a semi-global alignment since we won’t penalize gaps at the beginning and at the end of our pattern.We also compute a band allowing at most2 errors in either direction.Don’t forget to clear the gaps in each iteration, otherwise we might encounter wrong alignments.

TGapsgapsText;TGapsgapsPattern;assignSource(gapsPattern,pattern);std::cout<<"Text: "<<text<<"\tPattern: "<<pattern<<std::endl;for(TIteratorit=begin(locations);it!=end(locations);++it){// Clear previously computed gaps.clearGaps(gapsText);clearGaps(gapsPattern);// Only recompute the area within the current window over the text.TSequencetextInfix=infix(text,*it,*it+length(pattern));assignSource(gapsText,textInfix);// Use semi-global alignment since we do not want to track leading/trailing gaps in the pattern.// Restirct search space using a band allowing at most 2 errors in vertical/horizontal direction.intscore=globalAlignment(gapsText,gapsPattern,Score<int>(0,-1,-1),AlignConfig<true,false,false,true>(),-2,2);

In the next part we determine the cigar string for the matched pattern.We have to remove leading and trailing gaps in thegapsPattern object using the functionssetClippedBeginPosition andsetClippedEndPosition.We also need to set the clipped begin position for thegapsText object such that both Gaps begin at the same position.

TGapsIteratoritGapsPattern=begin(gapsPattern);TGapsIteratoritGapsEnd=end(gapsPattern);// Remove trailing gaps in pattern.intcount=0;while(isGap(--itGapsEnd))++count;setClippedEndPosition(gapsPattern,length(gapsPattern)-count);// Remove leading gaps in pattern.if(isGap(itGapsPattern)){setClippedBeginPosition(gapsPattern,countGaps(itGapsPattern));setClippedBeginPosition(gapsText,countGaps(itGapsPattern));}// Reinitilaize the iterators.TGapsIteratoritGapsText=begin(gapsText);itGapsPattern=begin(gapsPattern);itGapsEnd=end(gapsPattern);// Use a stringstream to construct the cigar string.std::stringstreamcigar;intnumMatchAndMismatches=0;while(itGapsPattern!=itGapsEnd){

First, we identify insertions using the functionsisGap andcountGaps.

// Count insertions.if(isGap(itGapsText)){intnumGaps=countGaps(itGapsText);cigar<<numGaps<<"I";itGapsText+=numGaps;itGapsPattern+=numGaps;continue;}

We do the same to identify deletions.

// Count deletions.if(isGap(itGapsPattern)){intnumGaps=countGaps(itGapsPattern);cigar<<numGaps<<"D";itGapsText+=numGaps;itGapsPattern+=numGaps;continue;}

If there is neither an insertion nor a deletion, then there must be a match or a mismatch. As long as we encounter matches and mismatches we move forward in the Gaps structures.Finally we print out the position of the hit, its total number of edits and the corresponding cigar string.

// Count matches and  mismatches.while(itGapsPattern!=itGapsEnd){if(isGap(itGapsPattern)||isGap(itGapsText))break;++numMatchAndMismatches;++itGapsText;++itGapsPattern;}if(numMatchAndMismatches!=0)cigar<<numMatchAndMismatches<<"M";numMatchAndMismatches=0;}// Output the hit position in the text, the total number of edits and the corresponding cigar string.std::cout<<"Hit at position  "<<*it<<"\ttotal edits: "<<abs(score)<<"\tcigar: "<<cigar.str()<<std::endl;}return0;}

Here is the output of this program.

Text: MISSISSIPPIANDMISSOURIPattern: SISSIHit at position  0total edits: 1cigar: 5MHit at position  1total edits: 1cigar: 1I4MHit at position  2total edits: 1cigar: 4M1IHit at position  3total edits: 0cigar: 5MHit at position  4total edits: 1cigar: 1I4MHit at position  6total edits: 2cigar: 5MHit at position  14total edits: 2cigar: 4M1I