consensus sequence from MSA using seqAn

Question

[Zev.Kronenberg]

I'm trying to create a consensus sequencing using seqAn. I've followed the nice guide to create a multiple sequence alignment:

http://seqan.readthedocs.org/en/latest/Tutorial/MultipleSequenceAlignment.html

I am looking for an explicit example of how turn the alignment into a consensus. I am also having trouble using the "alignmentEvaluation()" function. I made the assumptions that gapCount and other parameters are ints, but compilation doesn't agree.

Here is my working code so far:

  {
    using namespace seqan;

    typedef String< Dna > TSequence;
    StringSet<TSequence> seq;

    for(vector<string>::iterator seqs = s.begin();
        seqs != s.end(); seqs++
        ){
      appendValue(seq, *seqs);
    }
    Graph<Alignment<StringSet<TSequence, Dependent<> > > > aliG(seq);
    globalMsaAlignment(aliG, Blosum62(-1, -11));
    cerr << aliG << endl;

  }

dave.weese, any pointers or code snippets would be appreciated.

Answer 1

[Manuel]

Hi Zev, sorry for taking so long to answer.

Today, I had the time to flesh out the MSA tutorial and API documentation of the globalMsaAlignment() function more. Note that the links below go to the documentation of the develop branch but the functionality for the MSA is the same in the master branch and the 1.4.2 release:

• updated documentation of globalMsaAlignment function
• updated multiple sequence alignment tutorial

I updated the tutorial to perform a MSA directly into an Align object and added a second section that shows how to use ProfileChar for first a profile and then a consensus computation. Please indicate any unclear points so we can further improve the documentation.

Please do not hesitate to report more problems with the documentation, no matter whether it's bugs or unclear information. You will get the most responsive answers on our mailing list.

HTH,
Manuel

[Zev.Kronenberg]

This is ideal software support. Thank you kindly. SeqAn is awesome!

Answer 2

[Zev.Kronenberg]

This is an ugly solution, but it worked:

Creating the alignment:

   typedef String< Dna > TSequence;
    typedef Graph<Alignment<StringSet<TSequence, Dependent<> > > > TGraph;

    StringSet<TSequence> seq;

    for(vector<string>::iterator seqs = s.begin();
    seqs != s.end(); seqs++
        ){
      appendValue(seq, *seqs);
    }
    TGraph aliG(seq);
    globalMsaAlignment(aliG, Blosum62(-1, -11));

Looping over the alignment to create the consensus:

  String<char> align;
    convertAlignment(aliG, align);

    unsigned int nseq   = length(seq);
    unsigned int colLen = length(align) / nseq;

    stringstream con;

    //    cerr << nseq << "\t" << colLen << endl;

    for(unsigned int z = 0 ; z < colLen; z++){
      map<char, int> columnBases;
      for(unsigned int s = 0; s < nseq; s++){
        //      cerr << align[z + (s*colLen)] ;
        if(align[z + (s*colLen)] != gapValue<char>()){
          columnBases[align[z + (s*colLen)]]++;
        }
      }
      if(columnBases.size() == 1){

   con << columnBases.begin()->first;

      }
      else{
        con << "N";
      }
      // cerr << endl;
    }