source: branches/help/RNACMA/AlignedSequenceLoader.cxx

/*
 * AlignedSequenceLoader.cxx
 *
 *  Created on: Feb 15, 2010
 *      Author: Breno Faria
 *
 *  Institute of Microbiology (Technical University Munich)
 *  http://www.arb-home.de/
 */


#include "AlignedSequenceLoader.h"

#include <arbdbt.h>

/**
 * Loads the marked sequences aligned from Arb's DB.
 * This loader only considers letters given by the following regular
 * expression: [ACGTUacgtu]
 */
AlignedSequenceLoader::AlignedSequenceLoader(GBDATA *gb_main) :
    error(NULp),
    seqs(NULp)
{
    error = GB_push_transaction(gb_main);
    if (!error) {
        char *al_name = GBT_get_default_alignment(gb_main);
        int   al_len  = GBT_get_alignment_len(gb_main, al_name);

        if (al_len == -1) {
            error = GB_await_error();
        }

        if (!error) {
            seqs    = new VecVecType(0);
            MSA_len = al_len;

            size_t occurrences[MSA_len];
            for (size_t i = 0; i < MSA_len; i++) occurrences[i] = 0;

            cout << "loading marked species: ";
            flush(cout);
            for (GBDATA *gb_species = GBT_first_marked_species(gb_main);
                 gb_species && !error;
                 gb_species = GBT_next_marked_species(gb_species))
            {
                GBDATA *gb_data = GBT_find_sequence(gb_species, al_name);
                if (gb_data) { // existing alignment data for this species
                    cout << GBT_get_name_or_description(gb_species) << " ";
                    flush(cout);

                    string sequence = GB_read_char_pntr(gb_data);

                    string *seq_as_vec = new string[MSA_len];
                    int     k          = 0;
                    for (string::iterator i = sequence.begin(); i != sequence.end(); ++i) {
                        switch (*i) {
                            case 'A':
                            case 'a':
                                seq_as_vec[k] = "A";
                                occurrences[k] += 1;
                                break;
                            case 'C':
                            case 'c':
                                seq_as_vec[k] = "C";
                                occurrences[k] += 1;
                                break;
                            case 'G':
                            case 'g':
                                seq_as_vec[k] = "G";
                                occurrences[k] += 1;
                                break;
                            case 'T':
                            case 't':
                            case 'U':
                            case 'u':
                                seq_as_vec[k] = "T";
                                occurrences[k] += 1;
                                break;
                            default:
                                seq_as_vec[k] = "-";
                                break;
                        }
                        k++;
                    }

                    arb_assert((size_t)k == MSA_len);
                    vector<string> seq_vector(&seq_as_vec[0], &seq_as_vec[k]);
                    delete [] seq_as_vec;

                    seqs->push_back(seq_vector);
                }
                else {
                    error = GBS_global_string("species '%s' has no data in selected alignment '%s'",
                                              GBT_get_name_or_description(gb_species),
                                              al_name);
                }
            }

            cout << "done. Total number of species: " << seqs->size() << endl;
            flush(cout);

            cleanSeqs(occurrences, MSA_len);
        }

        free(al_name);
    }

    error = GB_end_transaction(gb_main, error);
}

/**
 * Returns the aligned seqs.
 *
 * @return the aligned seqs.
 */
VecVecType* AlignedSequenceLoader::getSequences() {
    arb_assert(!has_error());
    arb_assert(seqs);
    return seqs;
}

/**
 * Destructor.
 */
AlignedSequenceLoader::~AlignedSequenceLoader() {
    delete seqs;
}

/**
 * Getter for the MSA_len.
 *
 * @return the MSA length.
 */
size_t AlignedSequenceLoader::getMsaLen() {
    arb_assert(!has_error());
    return MSA_len;
}

/**
 * Getter for the position map. The position map maps positions in the clean
 * sequence to the original positions in the alignment.
 *
 * @return the position map.
 */
const vector<size_t>& AlignedSequenceLoader::getPositionMap() {
    arb_assert(!has_error());
    return position_map;
}

/**
 * This method cleans-up the empty positions of the MSA.
 *
 *
 * @param occurrences: a list gathering the number of occurrences of bases at
 *                     each position of the MSA.
 * @param len: the length of occurrences.
 */
void AlignedSequenceLoader::cleanSeqs(const size_t *occurrences, long len) {

    cout << "cleaning-up sequences of empty positions... " << endl;
    flush( cout);

    size_t num_of_bases = 0;
    for (int i = 0; i < len; i++) {
        if (occurrences[i] != 0) {
            num_of_bases++;
        }
    }

    cout << "number of non-empty positions in MSA: " << num_of_bases
         << ". Filtered out " << len - num_of_bases << " positions." << endl;

    VecVecType *clean_seqs = new VecVecType(0);

    cout << "computing position map..." << endl;
    position_map.resize(num_of_bases, 0);

    int j = 0;
    for (int i = 0; i < len; i++) {
        if (occurrences[i] != 0) {
            position_map.at(j) = i;
            j++;
        }
    }

    for (VecVecType::iterator seq = seqs->begin(); seq != seqs->end(); ++seq) {
        // for every sequence
        vector<string> sequence(num_of_bases, "");
        int jj = 0;
        for (int i = 0; i < len; ++i) {
            if (occurrences[i] != 0) {
                sequence.at(jj) = seq->at(i);
                jj++;
            }
        }
        arb_assert(sequence.size() == num_of_bases);
        clean_seqs->push_back(sequence);
    }

    delete seqs; // before overwriting it
    seqs = clean_seqs;

    cout << "clean-up done." << endl;

}