| 1 | // =============================================================== // |
|---|
| 2 | // // |
|---|
| 3 | // File : distanalyse.cxx // |
|---|
| 4 | // Purpose : // |
|---|
| 5 | // // |
|---|
| 6 | // Institute of Microbiology (Technical University Munich) // |
|---|
| 7 | // http://www.arb-home.de/ // |
|---|
| 8 | // // |
|---|
| 9 | // =============================================================== // |
|---|
| 10 | |
|---|
| 11 | #include "di_matr.hxx" |
|---|
| 12 | #include <AP_seq_dna.hxx> |
|---|
| 13 | #include <AP_filter.hxx> |
|---|
| 14 | #include <aw_awar.hxx> |
|---|
| 15 | #include <aw_msg.hxx> |
|---|
| 16 | #include <aw_root.hxx> |
|---|
| 17 | |
|---|
| 18 | |
|---|
| 19 | void DI_MATRIX::analyse() { |
|---|
| 20 | long row; |
|---|
| 21 | |
|---|
| 22 | long act_gci, mean_gci=0; |
|---|
| 23 | float act_gc, min_gc=9999.9, max_gc=0.0; |
|---|
| 24 | long act_len, mean_len=0, min_len=9999999, max_len=0; |
|---|
| 25 | |
|---|
| 26 | if (is_AA) { |
|---|
| 27 | if (nentries> 100) { |
|---|
| 28 | aw_message("A lot of sequences!\n ==> fast Kimura selected! (instead of PAM)"); |
|---|
| 29 | aw_root->awar(AWAR_DIST_CORR_TRANS)->write_int(DI_TRANSFORMATION_KIMURA); |
|---|
| 30 | } |
|---|
| 31 | else { |
|---|
| 32 | aw_message("Only limited number of sequences!\n" |
|---|
| 33 | " ==> slow PAM selected! (instead of Kimura)"); |
|---|
| 34 | aw_root->awar(AWAR_DIST_CORR_TRANS)->write_int(DI_TRANSFORMATION_PAM); |
|---|
| 35 | } |
|---|
| 36 | } |
|---|
| 37 | else { |
|---|
| 38 | // calculate meanvalue of sequencelength: |
|---|
| 39 | for (row=0; row<nentries; row++) { |
|---|
| 40 | act_gci = 0; |
|---|
| 41 | act_len = 0; |
|---|
| 42 | |
|---|
| 43 | const char *sequ = entries[row]->sequence_parsimony->get_sequence(); |
|---|
| 44 | |
|---|
| 45 | size_t flen = aliview->get_length(); |
|---|
| 46 | for (size_t pos=0; pos<flen; pos++) { |
|---|
| 47 | char ch = sequ[pos]; |
|---|
| 48 | if (ch == AP_C || ch == AP_G) act_gci++; |
|---|
| 49 | if (ch == AP_A || ch == AP_C || ch == AP_G || ch == AP_T) act_len++; |
|---|
| 50 | } |
|---|
| 51 | mean_gci += act_gci; |
|---|
| 52 | act_gc = ((float) act_gci) / act_len; |
|---|
| 53 | if (act_gc < min_gc) min_gc = act_gc; |
|---|
| 54 | if (act_gc > max_gc) max_gc = act_gc; |
|---|
| 55 | mean_len += act_len; |
|---|
| 56 | if (act_len < min_len) min_len = act_len; |
|---|
| 57 | if (act_len > max_len) max_len = act_len; |
|---|
| 58 | } |
|---|
| 59 | |
|---|
| 60 | if (min_len * 1.3 < max_len) { |
|---|
| 61 | aw_message("Warning: The length of sequences differs significantly!\n" |
|---|
| 62 | " Be careful: Neighbour Joining is sensitive to\n" |
|---|
| 63 | " this kind of \"error\""); |
|---|
| 64 | } |
|---|
| 65 | mean_len /= nentries; |
|---|
| 66 | |
|---|
| 67 | if (mean_len < 100) { |
|---|
| 68 | aw_message("Too short sequences!\n ==> No correction selected!"); |
|---|
| 69 | aw_root->awar(AWAR_DIST_CORR_TRANS)->write_int(DI_TRANSFORMATION_NONE); |
|---|
| 70 | } |
|---|
| 71 | else if (mean_len < 300) { |
|---|
| 72 | aw_message("Meanlength shorter than 300\n ==> Jukes Cantor selected!"); |
|---|
| 73 | aw_root->awar(AWAR_DIST_CORR_TRANS)->write_int(DI_TRANSFORMATION_JUKES_CANTOR); |
|---|
| 74 | } |
|---|
| 75 | else if ((mean_len < 1000) || ((max_gc / min_gc) < 1.2)) { |
|---|
| 76 | const char *reason; |
|---|
| 77 | if (mean_len < 1000) reason = "Sequences are too short for Olsen!"; |
|---|
| 78 | else reason = GBS_global_string("Maximal GC (%f) : Minimal GC (%f) < 1.2", max_gc, min_gc); |
|---|
| 79 | |
|---|
| 80 | reason = GBS_global_string("%s ==> Felsenstein selected!", reason); |
|---|
| 81 | aw_message(reason); |
|---|
| 82 | aw_root->awar(AWAR_DIST_CORR_TRANS)->write_int(DI_TRANSFORMATION_FELSENSTEIN); |
|---|
| 83 | } |
|---|
| 84 | else { |
|---|
| 85 | aw_message("Olsen selected!"); |
|---|
| 86 | aw_root->awar(AWAR_DIST_CORR_TRANS)->write_int(DI_TRANSFORMATION_OLSEN); |
|---|
| 87 | } |
|---|
| 88 | } |
|---|
| 89 | } |
|---|
