Question: consensus sequence from MSA using seqAn
1
gravatar for Zev.Kronenberg
3.9 years ago by
United States
Zev.Kronenberg11k wrote:

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. 

 

seqan c++ • 1.6k views
ADD COMMENTlink modified 3.9 years ago by Manuel340 • written 3.9 years ago by Zev.Kronenberg11k
3
gravatar for Manuel
3.9 years ago by
Manuel340
Germany
Manuel340 wrote:

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:

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

ADD COMMENTlink written 3.9 years ago by Manuel340
1

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

ADD REPLYlink written 3.9 years ago by Zev.Kronenberg11k
1
gravatar for Zev.Kronenberg
3.9 years ago by
United States
Zev.Kronenberg11k wrote:

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;
    }

 

 

 

ADD COMMENTlink written 3.9 years ago by Zev.Kronenberg11k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 725 users visited in the last hour