1. we discussed in detail three methods for pairwise alignment that utilize the

Question

1. we discussed in detail three methods for pairwise alignment that utilize the dynamic programming algorithm: Needleman-Wunsch, Waterman-Smith-Beyer, and Smith-Waterman algorithms. For each work and example by hand to implement the pairwise alignment for two short sequences.

2. Write a program two of the above algorithms for pairwise alignment. The program should take two input sequences and output all optimal alignments. The program should also have the ability to easily adjust the scoring function.

3. Try different scoring functions in your program to see how this changes the results.

4. Choose two example nucleotide sequences to align from http://www.ncbi.nlm.nih.gov/nucleotide/ and compare the alignment with the alignment generated by one of their alignment tools (or another tool) with the alignment generated by both of the algorithms you coded. Describe and explain any differences.

Explanation / Answer

#include<iostream>

#include<seqan/align.h>

using namespace seqan;

int main()

{

//Let’s start with initializing the Align object to contain the two sequences.

   Align<String<char> > ali;

    resize(rows(ali), 2);

    assignSource(row(ali, 0), "aphilologicaltheorem");

    assignSource(row(ali, 1), "bizarreamphibology");

/* Now the best alignment given the scoring parameters is computed using the Dynamic Gap model by the function localAlignment. The returned score value is printed directly, and the alignment itself in the next line. The functions clippedBeginPosition and clippedEndPosition can be used to retrieve the begin and end position of the matching subsequences within the original sequences.*/

   std::cout << "Score = " << localAlignment(ali, Score<int>(3, -3, -2, -2), DynamicGaps()) << std::endl;

    std::cout << ali;

    std::cout << "Aligns Seq1[" << clippedBeginPosition(row(ali, 0)) << ":" << (clippedEndPosition(row(ali, 0)) - 1) << "]";

    std::cout << " and Seq2[" << clippedBeginPosition(row(ali, 1)) << ":" << (clippedEndPosition(row(ali, 1)) - 1) << "]" << std::endl << std::endl;

/* Next, several local alignments of the two given DNA sequences are going to be computed. First, the Align object is created.*/

   Align<String<Dna> > ali2;

    resize(rows(ali2), 2);

    assignSource(row(ali2, 0), "ataagcgtctcg");

    assignSource(row(ali2, 1), "tcatagagttgc");

/*A LocalAlignmentEnumerator object needs to be initialized on the Align object. In addition to the Align object and the scoring scheme, we now also pass the finder and a minimal score value, 4 in this case, to the localAlignment function. The WatermanEggert tag specifies the desired Waterman-Eggert algorithm. While the score of the local alignment satisfies the minimal score cutoff, the alignments are printed with their scores and the subsequence begin and end positions.*/

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

    LocalAlignmentEnumerator<Score<int>, Unbanded> enumerator(scoring, 5);

    while (nextLocalAlignment(ali2, enumerator))

    {

        std::cout << "Score = " << getScore(enumerator) << std::endl;

        std::cout << ali2;

        std::cout << "Aligns Seq1[" << clippedBeginPosition(row(ali2, 0)) << ":" << (clippedEndPosition(row(ali2, 0)) - 1) << "]";

        std::cout << " and Seq2[" << clippedBeginPosition(row(ali2, 1)) << ":" << (clippedEndPosition(row(ali2, 1)) - 1) << "]" << std::endl << std::endl;

    }

    return 0;

}

Related Questions

Navigate

• Browse (All)
• Browse 1
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.