| 1 | // =============================================================== // |
|---|
| 2 | // // |
|---|
| 3 | // File : ad_spec.cxx // |
|---|
| 4 | // Purpose : Create and modify species and SAI. // |
|---|
| 5 | // // |
|---|
| 6 | // Institute of Microbiology (Technical University Munich) // |
|---|
| 7 | // http://www.arb-home.de/ // |
|---|
| 8 | // // |
|---|
| 9 | // =============================================================== // |
|---|
| 10 | |
|---|
| 11 | #include "NT_local.h" |
|---|
| 12 | #include "map_viewer.h" |
|---|
| 13 | |
|---|
| 14 | #include <dbui.h> |
|---|
| 15 | |
|---|
| 16 | #include <www.hxx> |
|---|
| 17 | |
|---|
| 18 | #include <aw_awars.hxx> |
|---|
| 19 | #include <aw_msg.hxx> |
|---|
| 20 | #include <aw_root.hxx> |
|---|
| 21 | |
|---|
| 22 | #include <arb_progress.h> |
|---|
| 23 | #include <arb_defs.h> |
|---|
| 24 | |
|---|
| 25 | #include <cctype> |
|---|
| 26 | |
|---|
| 27 | static const char * const SAI_COUNTED_CHARS = "COUNTED_CHARS"; |
|---|
| 28 | |
|---|
| 29 | void NT_count_different_chars(AW_window*, GBDATA *gb_main) { |
|---|
| 30 | // @@@ extract algorithm and call extracted from testcode |
|---|
| 31 | ARB_ERROR error = GB_begin_transaction(gb_main); // open transaction |
|---|
| 32 | |
|---|
| 33 | if (!error) { |
|---|
| 34 | char *alignment_name = GBT_get_default_alignment(gb_main); |
|---|
| 35 | if (!alignment_name) { |
|---|
| 36 | error = GB_await_error(); |
|---|
| 37 | } |
|---|
| 38 | else { |
|---|
| 39 | int alignment_len = GBT_get_alignment_len(gb_main, alignment_name); |
|---|
| 40 | if (alignment_len<=0) { |
|---|
| 41 | error = GB_await_error(); |
|---|
| 42 | } |
|---|
| 43 | |
|---|
| 44 | if (!error) { |
|---|
| 45 | int is_amino = GBT_is_alignment_protein(gb_main, alignment_name); |
|---|
| 46 | |
|---|
| 47 | const int MAXLETTER = 256; |
|---|
| 48 | const int FIRSTLETTER = 0; // cppcheck-suppress variableScope |
|---|
| 49 | |
|---|
| 50 | typedef bool letterOccurs[MAXLETTER]; |
|---|
| 51 | |
|---|
| 52 | letterOccurs *occurs = ARB_alloc<letterOccurs>(alignment_len); |
|---|
| 53 | for (int i = 0; i<MAXLETTER; ++i) { // LOOP_VECTORIZED[!<6.0,!>8.0] |
|---|
| 54 | for (int p = 0; p<alignment_len; ++p) { // LOOP_VECTORIZED[!<8.0] |
|---|
| 55 | occurs[p][i] = false; |
|---|
| 56 | } |
|---|
| 57 | } |
|---|
| 58 | |
|---|
| 59 | // loop over all marked species |
|---|
| 60 | { |
|---|
| 61 | long all_marked = GBT_count_marked_species(gb_main); |
|---|
| 62 | arb_progress progress("Counting different characters", all_marked); |
|---|
| 63 | |
|---|
| 64 | for (GBDATA *gb_species = GBT_first_marked_species(gb_main); |
|---|
| 65 | gb_species && !error; |
|---|
| 66 | gb_species = GBT_next_marked_species(gb_species)) |
|---|
| 67 | { |
|---|
| 68 | GBDATA *gb_ali = GB_entry(gb_species, alignment_name); // search the sequence database entry ( ali_xxx/data ) |
|---|
| 69 | if (gb_ali) { |
|---|
| 70 | GBDATA *gb_data = GB_entry(gb_ali, "data"); |
|---|
| 71 | if (gb_data) { |
|---|
| 72 | const char * const seq = GB_read_char_pntr(gb_data); |
|---|
| 73 | if (seq) { |
|---|
| 74 | for (int i=0; i< alignment_len; ++i) { |
|---|
| 75 | unsigned char c = seq[i]; |
|---|
| 76 | if (!c) break; |
|---|
| 77 | |
|---|
| 78 | occurs[i][c-FIRSTLETTER] = true; |
|---|
| 79 | } |
|---|
| 80 | } |
|---|
| 81 | } |
|---|
| 82 | } |
|---|
| 83 | progress.inc_and_check_user_abort(error); |
|---|
| 84 | } |
|---|
| 85 | } |
|---|
| 86 | |
|---|
| 87 | if (!error) { |
|---|
| 88 | |
|---|
| 89 | char filter[256]; |
|---|
| 90 | if (is_amino) for (int c = 0; c<256; ++c) filter[c] = isupper(c) && !strchr("BJOUZ", c); |
|---|
| 91 | else for (int c = 0; c<256; ++c) filter[c] = bool(strchr("ACGTU", c)); |
|---|
| 92 | |
|---|
| 93 | char result[alignment_len+1]; |
|---|
| 94 | for (int i=0; i<alignment_len; i++) { |
|---|
| 95 | int sum = 0; |
|---|
| 96 | for (int c = 'A'; c < 'Z'; ++c) { |
|---|
| 97 | if (filter[c]) { |
|---|
| 98 | sum += (occurs[i][c] || occurs[i][tolower(c)]); |
|---|
| 99 | } |
|---|
| 100 | } |
|---|
| 101 | result[i] = sum<10 ? '0'+sum : 'A'-10+sum; |
|---|
| 102 | } |
|---|
| 103 | result[alignment_len] = 0; |
|---|
| 104 | |
|---|
| 105 | { |
|---|
| 106 | GBDATA *gb_sai = GBT_find_or_create_SAI(gb_main, SAI_COUNTED_CHARS); |
|---|
| 107 | if (!gb_sai) error = GB_await_error(); |
|---|
| 108 | else { |
|---|
| 109 | GBDATA *gb_data = GBT_add_data(gb_sai, alignment_name, "data", GB_STRING); |
|---|
| 110 | if (!gb_data) error = GB_await_error(); |
|---|
| 111 | else error = GB_write_string(gb_data, result); |
|---|
| 112 | } |
|---|
| 113 | } |
|---|
| 114 | } |
|---|
| 115 | |
|---|
| 116 | free(occurs); |
|---|
| 117 | } |
|---|
| 118 | |
|---|
| 119 | free(alignment_name); |
|---|
| 120 | } |
|---|
| 121 | } |
|---|
| 122 | |
|---|
| 123 | GB_end_transaction_show_error(gb_main, error, aw_message); |
|---|
| 124 | } |
|---|
| 125 | |
|---|
| 126 | void NT_create_sai_from_pfold(AW_window *aww) { |
|---|
| 127 | /*! \brief Creates an SAI from protein secondary structure of a selected species. |
|---|
| 128 | * |
|---|
| 129 | * \param[in] aww AW_window |
|---|
| 130 | * |
|---|
| 131 | * The function takes the currently selected species and searches for the field |
|---|
| 132 | * "sec_struct". A new SAI is created using the data in this field. A simple input |
|---|
| 133 | * window allows the user to change the default name ([species name]_pfold) for |
|---|
| 134 | * the new SAI. |
|---|
| 135 | * |
|---|
| 136 | * \note The import filter "dssp_all.ift" allows for importing the amino acid sequence |
|---|
| 137 | * as well as the protein secondary structure from a dssp file and the structure |
|---|
| 138 | * is stored in the field "sec_struct". That way, secondary structure can be |
|---|
| 139 | * aligned along with the sequence manually and can later be extracted to create |
|---|
| 140 | * an SAI. |
|---|
| 141 | * |
|---|
| 142 | * \attention The import filter "dssp_2nd_struct.ift" extracts only the protein |
|---|
| 143 | * secondary structure which is stored as alignment data. SAIs can simply |
|---|
| 144 | * be created from these species via move_species_to_extended(). |
|---|
| 145 | */ |
|---|
| 146 | |
|---|
| 147 | GB_begin_transaction(GLOBAL.gb_main); |
|---|
| 148 | |
|---|
| 149 | GB_ERROR error = NULp; |
|---|
| 150 | char *sai_name = NULp; |
|---|
| 151 | char *sec_struct = NULp; |
|---|
| 152 | bool canceled = false; |
|---|
| 153 | |
|---|
| 154 | // get the selected species |
|---|
| 155 | char *species_name = aww->get_root()->awar(AWAR_SPECIES_NAME)->read_string(); |
|---|
| 156 | GBDATA *gb_species = NULp; |
|---|
| 157 | if (!strcmp(species_name, "") || !(gb_species = GBT_find_species(GLOBAL.gb_main, species_name))) { |
|---|
| 158 | error = "Please select a species first."; |
|---|
| 159 | } |
|---|
| 160 | else { |
|---|
| 161 | // search for the field "sec_struct" |
|---|
| 162 | GBDATA *gb_species_sec_struct = GB_entry(gb_species, "sec_struct"); |
|---|
| 163 | if (!gb_species_sec_struct) { |
|---|
| 164 | error = "Field \"sec_struct\" not found or empty. Please select another species."; |
|---|
| 165 | } |
|---|
| 166 | else if (!(sec_struct = GB_read_string(gb_species_sec_struct))) { |
|---|
| 167 | error = "Couldn't read field \"sec_struct\". Is it empty?"; |
|---|
| 168 | } |
|---|
| 169 | else { |
|---|
| 170 | // generate default name and name input field for the new SAI |
|---|
| 171 | { |
|---|
| 172 | char *sai_default_name = GBS_global_string_copy("%s%s", species_name, strstr(species_name, "_pfold") ? "" : "_pfold"); |
|---|
| 173 | sai_name = aw_input("Name of SAI to create:", sai_default_name); |
|---|
| 174 | free(sai_default_name); |
|---|
| 175 | } |
|---|
| 176 | |
|---|
| 177 | if (!sai_name) { |
|---|
| 178 | canceled = true; |
|---|
| 179 | } |
|---|
| 180 | else if (strspn(sai_name, " ") == strlen(sai_name)) { |
|---|
| 181 | error = "Name of SAI is empty. Please enter a valid name."; |
|---|
| 182 | } |
|---|
| 183 | else { |
|---|
| 184 | GBDATA *gb_sai_data = GBT_get_SAI_data(GLOBAL.gb_main); |
|---|
| 185 | GBDATA *gb_sai = GBT_find_SAI_rel_SAI_data(gb_sai_data, sai_name); |
|---|
| 186 | |
|---|
| 187 | if (gb_sai) { |
|---|
| 188 | error = "SAI with the same name already exists. Please enter another name."; |
|---|
| 189 | } |
|---|
| 190 | else { |
|---|
| 191 | char *ali_name = GBT_get_default_alignment(GLOBAL.gb_main); |
|---|
| 192 | if (!ali_name) { |
|---|
| 193 | error = GB_await_error(); |
|---|
| 194 | } |
|---|
| 195 | else { |
|---|
| 196 | // create SAI container and copy fields from the species to the SAI |
|---|
| 197 | gb_sai = GB_create_container(gb_sai_data, "extended"); |
|---|
| 198 | GBDATA *gb_species_field = GB_child(gb_species); |
|---|
| 199 | |
|---|
| 200 | while (gb_species_field && !error) { |
|---|
| 201 | char *key = GB_read_key(gb_species_field); |
|---|
| 202 | GBDATA *gb_sai_field = GB_search(gb_sai, GB_read_key(gb_species_field), GB_read_type(gb_species_field)); |
|---|
| 203 | |
|---|
| 204 | if (strcmp(key, "name") == 0) { // write the new name |
|---|
| 205 | error = GB_write_string(gb_sai_field, sai_name); |
|---|
| 206 | } |
|---|
| 207 | else if (strcmp(key, "sec_struct") == 0) { // write contents from the field "sec_struct" to the alignment data |
|---|
| 208 | GBDATA *gb_sai_ali = GB_search(gb_sai, ali_name, GB_CREATE_CONTAINER); |
|---|
| 209 | if (!gb_sai_ali) error = GB_await_error(); |
|---|
| 210 | else error = GBT_write_string(gb_sai_ali, "data", sec_struct); |
|---|
| 211 | } |
|---|
| 212 | else if (strcmp(key, "acc") != 0 && strcmp(key, ali_name) != 0) { // don't copy "acc" and the old alignment data |
|---|
| 213 | error = GB_copy_dropProtectMarksAndTempstate(gb_sai_field, gb_species_field); |
|---|
| 214 | } |
|---|
| 215 | gb_species_field = GB_nextChild(gb_species_field); |
|---|
| 216 | free(key); |
|---|
| 217 | } |
|---|
| 218 | |
|---|
| 219 | // generate accession number and delete field "sec_struct" from the SAI |
|---|
| 220 | if (!error) { |
|---|
| 221 | // TODO: is it necessary that a new acc is generated here? |
|---|
| 222 | GBDATA *gb_sai_acc = GB_search(gb_sai, "acc", GB_FIND); |
|---|
| 223 | if (gb_sai_acc) { |
|---|
| 224 | GB_delete(gb_sai_acc); |
|---|
| 225 | GBT_gen_accession_number(gb_sai, ali_name); |
|---|
| 226 | } |
|---|
| 227 | GBDATA *gb_sai_sec_struct = GB_search(gb_sai, "sec_struct", GB_FIND); |
|---|
| 228 | if (gb_sai_sec_struct) GB_delete(gb_sai_sec_struct); |
|---|
| 229 | aww->get_root()->awar(AWAR_SAI_NAME)->write_string(sai_name); |
|---|
| 230 | } |
|---|
| 231 | free(ali_name); |
|---|
| 232 | } |
|---|
| 233 | } |
|---|
| 234 | } |
|---|
| 235 | } |
|---|
| 236 | } |
|---|
| 237 | |
|---|
| 238 | if (canceled) error = "Aborted by user"; |
|---|
| 239 | |
|---|
| 240 | GB_end_transaction_show_error(GLOBAL.gb_main, error, aw_message); |
|---|
| 241 | |
|---|
| 242 | if (!error) { |
|---|
| 243 | AW_window *sai_info = NT_create_extendeds_window(aww->get_root()); |
|---|
| 244 | // TODO: why doesn't info box show anything on first startup? proper refresh needed? |
|---|
| 245 | sai_info->activate(); |
|---|
| 246 | } |
|---|
| 247 | |
|---|
| 248 | free(species_name); |
|---|
| 249 | free(sai_name); |
|---|
| 250 | free(sec_struct); |
|---|
| 251 | } |
|---|
| 252 | |
|---|
| 253 | void launch_MapViewer_cb(GBDATA *gbd, AD_MAP_VIEWER_TYPE type) { |
|---|
| 254 | // Note: sync with ../PARSIMONY/PARS_main.cxx@PARS_map_viewer |
|---|
| 255 | |
|---|
| 256 | GB_ERROR error = GB_push_transaction(GLOBAL.gb_main); |
|---|
| 257 | |
|---|
| 258 | if (!error) { |
|---|
| 259 | const char *species_name = ""; |
|---|
| 260 | GBDATA *gb_species_data = GBT_get_species_data(GLOBAL.gb_main); |
|---|
| 261 | |
|---|
| 262 | if (gbd && GB_get_father(gbd) == gb_species_data) { |
|---|
| 263 | species_name = GBT_get_name_or_description(gbd); |
|---|
| 264 | } |
|---|
| 265 | |
|---|
| 266 | AW_root::SINGLETON->awar(AWAR_SPECIES_NAME)->write_string(species_name); |
|---|
| 267 | } |
|---|
| 268 | |
|---|
| 269 | if (!error && gbd && type == ADMVT_WWW) { |
|---|
| 270 | error = awt_openDefaultURL_with_item(GLOBAL.aw_root, GLOBAL.gb_main, gbd); |
|---|
| 271 | } |
|---|
| 272 | |
|---|
| 273 | error = GB_end_transaction(GLOBAL.gb_main, error); |
|---|
| 274 | if (error) aw_message(error); |
|---|
| 275 | } |
|---|
| 276 | |
|---|
| 277 | |
|---|
| 278 | // -------------------------------------------------------------------------------- |
|---|
| 279 | |
|---|
| 280 | #ifdef UNIT_TESTS |
|---|
| 281 | #include <test_unit.h> |
|---|
| 282 | #include <arb_unit_test.h> |
|---|
| 283 | |
|---|
| 284 | static uint32_t counted_chars_checksum(GBDATA *gb_main) { |
|---|
| 285 | GB_transaction ta(gb_main); |
|---|
| 286 | |
|---|
| 287 | GBDATA *gb_sai; |
|---|
| 288 | GBDATA *gb_ali; |
|---|
| 289 | GBDATA *gb_counted_chars; |
|---|
| 290 | |
|---|
| 291 | char *ali_name = GBT_get_default_alignment(gb_main); // potential error will be detected by TEST_EXPECT_RESULT__NOERROREXPORTED below |
|---|
| 292 | |
|---|
| 293 | TEST_EXPECT_RESULT__NOERROREXPORTED(gb_sai = GBT_expect_SAI(gb_main, SAI_COUNTED_CHARS)); |
|---|
| 294 | TEST_EXPECT_RESULT__NOERROREXPORTED(gb_ali = GB_entry(gb_sai, ali_name)); |
|---|
| 295 | TEST_EXPECT_RESULT__NOERROREXPORTED(gb_counted_chars = GB_entry(gb_ali, "data")); |
|---|
| 296 | |
|---|
| 297 | const char *data = GB_read_char_pntr(gb_counted_chars); |
|---|
| 298 | |
|---|
| 299 | free(ali_name); |
|---|
| 300 | |
|---|
| 301 | return GBS_checksum(data, 0, NULp); |
|---|
| 302 | } |
|---|
| 303 | |
|---|
| 304 | void TEST_count_chars() { |
|---|
| 305 | // calculate SAI for test DBs |
|---|
| 306 | |
|---|
| 307 | arb_suppress_progress silence; |
|---|
| 308 | GB_shell shell; |
|---|
| 309 | |
|---|
| 310 | for (int prot = 0; prot<2; ++prot) { |
|---|
| 311 | GBDATA *gb_main; |
|---|
| 312 | TEST_EXPECT_RESULT__NOERROREXPORTED(gb_main = GB_open(prot ? "TEST_prot.arb" : "TEST_nuc.arb", "rw")); |
|---|
| 313 | |
|---|
| 314 | GBT_mark_all(gb_main, 1); |
|---|
| 315 | NT_count_different_chars(NULp, gb_main); |
|---|
| 316 | |
|---|
| 317 | uint32_t expected = prot ? 0x9cad14cc : 0xefb05e4e; |
|---|
| 318 | TEST_EXPECT_EQUAL(counted_chars_checksum(gb_main), expected); |
|---|
| 319 | |
|---|
| 320 | GB_close(gb_main); |
|---|
| 321 | } |
|---|
| 322 | } |
|---|
| 323 | void TEST_SLOW_count_chars() { |
|---|
| 324 | // calculate a real big test alignment |
|---|
| 325 | // |
|---|
| 326 | // the difference to TEST_count_chars() is just in size of alignment. |
|---|
| 327 | // NT_count_different_chars() crashes for big alignments when running in gdb |
|---|
| 328 | arb_suppress_progress silence; |
|---|
| 329 | GB_shell shell; |
|---|
| 330 | { |
|---|
| 331 | arb_unit_test::test_alignment_data data_source[] = { |
|---|
| 332 | { 1, "s1", "ACGT" }, |
|---|
| 333 | { 1, "s2", "ACGTN" }, |
|---|
| 334 | { 1, "s3", "NANNAN" }, |
|---|
| 335 | { 1, "s4", "GATTACA" }, |
|---|
| 336 | }; |
|---|
| 337 | |
|---|
| 338 | const int alilen = 50000; |
|---|
| 339 | const int count = ARRAY_ELEMS(data_source); |
|---|
| 340 | |
|---|
| 341 | char *longSeq[count]; |
|---|
| 342 | for (int c = 0; c<count; ++c) { |
|---|
| 343 | char *dest = longSeq[c] = ARB_alloc<char>(alilen+1); |
|---|
| 344 | |
|---|
| 345 | const char *source = data_source[c].data; |
|---|
| 346 | int len = strlen(source); |
|---|
| 347 | |
|---|
| 348 | for (int p = 0; p<alilen; ++p) { |
|---|
| 349 | dest[p] = source[p%len]; |
|---|
| 350 | } |
|---|
| 351 | dest[alilen] = 0; |
|---|
| 352 | |
|---|
| 353 | data_source[c].data = dest; |
|---|
| 354 | } |
|---|
| 355 | |
|---|
| 356 | ARB_ERROR error; |
|---|
| 357 | GBDATA *gb_main = TEST_CREATE_DB(error, "ali_test", data_source, false); |
|---|
| 358 | |
|---|
| 359 | TEST_EXPECT_NO_ERROR(error.deliver()); |
|---|
| 360 | |
|---|
| 361 | NT_count_different_chars(NULp, gb_main); |
|---|
| 362 | |
|---|
| 363 | // TEST_EXPECT_EQUAL(counted_chars_checksum(gb_main), 0x1d34a14f); |
|---|
| 364 | // TEST_EXPECT_EQUAL(counted_chars_checksum(gb_main), 0x609d788b); |
|---|
| 365 | TEST_EXPECT_EQUAL(counted_chars_checksum(gb_main), 0xccdfa527); |
|---|
| 366 | |
|---|
| 367 | for (int c = 0; c<count; ++c) { |
|---|
| 368 | free(longSeq[c]); |
|---|
| 369 | } |
|---|
| 370 | |
|---|
| 371 | GB_close(gb_main); |
|---|
| 372 | } |
|---|
| 373 | } |
|---|
| 374 | |
|---|
| 375 | #endif // UNIT_TESTS |
|---|
| 376 | |
|---|
