source: tags/svn.1.5.4/RNACMA/AlignedSequenceLoader.cxx

Last change on this file was 8322, checked in by westram, 14 years ago
  • clang-fix: replace auto-arrays of non-POD types by dynamic memory
File size: 5.0 KB
Line 
1/*
2 * AlignedSequenceLoader.cxx
3 *
4 *  Created on: Feb 15, 2010
5 *      Author: Breno Faria
6 */
7
8
9#include "AlignedSequenceLoader.h"
10
11#include <arbdb.h>
12#include <arbdbt.h>
13#include "dbconn.h"
14
15/**
16 * Loads the marked sequences aligned from Arb's DB.
17 * This loader only considers letters given by the following regular
18 * expression: [ACGTUacgtu]
19 */
20AlignedSequenceLoader::AlignedSequenceLoader() {
21
22    GBDATA *gb_main = runningDatabase();
23
24    GB_ERROR error = GB_push_transaction(gb_main);
25    if (error)
26        exit(1);
27
28    seqs = new VecVecType(0);
29
30    char *al_name = GBT_get_default_alignment(gb_main);
31    MSA_len = GBT_get_alignment_len(gb_main, al_name);
32
33    size_t occurrences[MSA_len];
34    for (size_t i = 0; i < MSA_len; i++)
35        occurrences[i] = 0;
36
37    cout << "loading marked species: ";
38    flush( cout);
39    for (GBDATA *gb_species = GBT_first_marked_species(gb_main); gb_species; gb_species
40             = GBT_next_marked_species(gb_species)) {
41
42        GBDATA *gb_ali = GB_entry(gb_species, al_name);
43
44        if (gb_ali) { // existing alignment for this species
45            GBDATA *gb_data = GB_entry(gb_ali, "data");
46
47            if (gb_data) {
48
49                cout << GBT_read_name(gb_species) << " ";
50                flush(cout);
51
52                string sequence = GB_read_string(gb_data);
53
54                string *seq_as_vec = new string[MSA_len];
55                int     k          = 0;
56                for (string::iterator i = sequence.begin(); i != sequence.end(); ++i) {
57                    switch (*i) {
58                        case 'A':
59                        case 'a':
60                            seq_as_vec[k] = "A";
61                            occurrences[k] += 1;
62                            break;
63                        case 'C':
64                        case 'c':
65                            seq_as_vec[k] = "C";
66                            occurrences[k] += 1;
67                            break;
68                        case 'G':
69                        case 'g':
70                            seq_as_vec[k] = "G";
71                            occurrences[k] += 1;
72                            break;
73                        case 'T':
74                        case 't':
75                        case 'U':
76                        case 'u':
77                            seq_as_vec[k] = "T";
78                            occurrences[k] += 1;
79                            break;
80                        default:
81                            seq_as_vec[k] = "-";
82                            break;
83                    }
84                    k++;
85                }
86
87                assert((size_t)k == MSA_len);
88                vector<string> seq_vector(&seq_as_vec[0], &seq_as_vec[k]);
89                delete [] seq_as_vec;
90
91                seqs->push_back(seq_vector);
92            }
93        }
94
95    }
96
97    GB_pop_transaction(gb_main);
98
99    cout << "done. Total number of species: " << seqs->size() << endl;
100    flush(cout);
101
102    cleanSeqs(occurrences, MSA_len);
103}
104
105/**
106 * Returns the aligned seqs.
107 *
108 * @return the aligned seqs.
109 */
110VecVecType* AlignedSequenceLoader::getSequences() {
111    return seqs;
112}
113
114/**
115 * Destructor.
116 */
117AlignedSequenceLoader::~AlignedSequenceLoader() {
118    delete seqs;
119}
120
121/**
122 * Getter for the MSA_len.
123 *
124 * @return the MSA length.
125 */
126size_t AlignedSequenceLoader::getMsaLen() {
127    return MSA_len;
128}
129
130/**
131 * Getter for the position map. The position map maps positions in the clean
132 * sequence to the original positions in the alignment.
133 *
134 * @return the position map.
135 */
136vector<size_t> * AlignedSequenceLoader::getPositionMap() {
137    return position_map;
138}
139
140/**
141 * This method cleans-up the empty positions of the MSA.
142 *
143 *
144 * @param occurrences: a list gathering the number of occurrences of bases at
145 *                     each position of the MSA.
146 * @param len: the length of occurrences.
147 */
148void AlignedSequenceLoader::cleanSeqs(size_t *occurrences, long len) {
149
150    cout << "cleaning-up sequences of empty positions... " << endl;
151    flush( cout);
152
153    size_t num_of_bases = 0;
154    for (int i = 0; i < len; i++) {
155        if (occurrences[i] != 0) {
156            num_of_bases++;
157        }
158    }
159
160    cout << "number of non-empty positions in MSA: " << num_of_bases
161         << ". Filtered out " << len - num_of_bases << " positions." << endl;
162
163    VecVecType *clean_seqs = new VecVecType(0);
164
165    cout << "computing position map..." << endl;
166    position_map = new vector<size_t> (num_of_bases, 0);
167    int j = 0;
168    for (int i = 0; i < len; i++) {
169        if (occurrences[i] != 0) {
170            position_map->at(j) = i;
171            j++;
172        }
173    }
174
175    for (VecVecType::iterator seq = seqs->begin(); seq != seqs->end(); ++seq) {
176        //for every sequence
177        vector<string> sequence(num_of_bases, "");
178        int jj = 0;
179        for (int i = 0; i < len; ++i) {
180            if (occurrences[i] != 0) {
181                sequence.at(jj) = seq->at(i);
182                jj++;
183            }
184        }
185        assert(sequence.size() == num_of_bases);
186        clean_seqs->push_back(sequence);
187    }
188
189    seqs = clean_seqs;
190
191    cout << "clean-up done." << endl;
192
193}
194
Note: See TracBrowser for help on using the repository browser.