pvp.cxx

Go to the documentation of this file.

// ============================================================= //
//                                                               //
//   File      : pvp.cxx                                         //
//   Purpose   : calculate positional variability (parsimony)    //
//                                                               //
//   Institute of Microbiology (Technical University Munich)     //
//   http://www.arb-home.de/                                     //
//                                                               //
// ============================================================= //

#include "pvp.h"
#include "AP_pos_var.h"

#include <TreeNode.h>
#include <arb_progress.h>

GB_ERROR PVP_calculate(GBDATA *gb_main, const char *ali_name, const char *tree_name, const char *target_SAI_name) {
    GB_ERROR       error   = NULp;
    GB_transaction ta(gb_main);
    long           ali_len = GBT_get_alignment_len(gb_main, ali_name);

    if (ali_len <= 0) {
        error = GB_await_error();
    }
    else {
        const double phase1_fraction = 1.0/650;
        arb_progress progress(WEIGHTED, "Calculating positional variability", phase1_fraction);
        progress.subtitle("Loading Tree");

        // get tree
        TreeNode *tree;
        {
            tree = GBT_read_tree(gb_main, tree_name, new SimpleRoot);
            if (!tree) {
                error = GB_await_error();
            }
            else {
                GBT_link_tree(tree, gb_main, true, NULp, NULp);
            }
        }

        progress.inc_and_check_user_abort(error);

        if (!error) {
            progress.subtitle("Counting mutations");

            GB_alignment_type at = GBT_get_alignment_type(gb_main, ali_name);
            arb_assert(at != GB_AT_UNKNOWN);

            bool       isNUC = at==GB_AT_DNA || at==GB_AT_RNA;
            AP_pos_var pvp(gb_main, ali_name, ali_len, isNUC, tree_name);

            error             = pvp.delete_aliEntry_from_SAI(target_SAI_name);
            if (!error) error = pvp.retrieve(tree);
            if (!error) error = pvp.save_aliEntry_to_SAI(target_SAI_name);

            progress.inc_and_check_user_abort(error);
        }

        destroy(tree);
    }

    ta.close(error);

    return error;
}

// --------------------------------------------------------------------------------

#ifdef UNIT_TESTS
#ifndef TEST_UNIT_H
#include <test_unit.h>
#endif

static arb_test::match_expectation saidata_equal(GBDATA *gb_main, const char *sainame, const char *aliname, const char *expected_data) {
    using namespace arb_test;

    GB_transaction ta(gb_main);

    GBDATA *gb_sai     = GBT_find_SAI(gb_main, sainame);
    GBDATA *gb_saiali  = GB_entry(gb_sai, aliname);
    GBDATA *gb_saidata = GB_entry(gb_saiali, "data");

    const char *saidata = GB_read_char_pntr(gb_saidata);

    return that(saidata).is_equal_to(expected_data);
}
static arb_test::match_expectation saitype_equal(GBDATA *gb_main, const char *sainame, const char *aliname, const char *expected_type) {
    using namespace arb_test;

    GB_transaction ta(gb_main);

    GBDATA *gb_sai     = GBT_find_SAI(gb_main, sainame);
    GBDATA *gb_saiali  = GB_entry(gb_sai, aliname);
    GBDATA *gb_saitype = GB_entry(gb_saiali, "_TYPE");

    const char *saitype = GB_read_char_pntr(gb_saitype);

    return that(saitype).is_equal_to(expected_type);
}

#define TEST_EXPECT_SAIDATA_EQUAL(expected) TEST_EXPECTATION(saidata_equal(gb_main,SAI_name,aliname,expected))
#define TEST_EXPECT_SAITYPE_EQUAL(expected) TEST_EXPECTATION(saitype_equal(gb_main,SAI_name,aliname,expected))

void TEST_pvp() {
    GB_shell    shell;
    const char *SAI_name = "POS_VAR_BY_PARSIMONY";

    { // nuc data:
        GBDATA     *gb_main = GB_open("TEST_nuc.arb", "rw");
        const char *aliname = "ali_16s";
        TEST_EXPECT_NO_ERROR(PVP_calculate(gb_main, aliname, "tree_nuc", SAI_name));

        TEST_EXPECT_SAITYPE_EQUAL("PVP: Positional Variability by Parsimony: tree 'tree_nuc' ntaxa 12");
        TEST_EXPECT_SAIDATA_EQUAL(".----222----7------32-7-7----774-----77-----747-77-7-47-----------73--433-.3134--1-41-4-4011222337--"
                                  "....-632320124--11141-4-..24312434--3----7-77---7----774-4-4-34--44-77-732777424--7-----3247------44"
                                  "33----------2--73063442221633737777326110244-3----73--2--4-773--4-7--3774777137-7774437274-3372--7--"
                                  "--7------47-7---24--73--------2-4--7374---74-734------43-777-47-7---3---3----------7-3---3----------"
                                  "-----7-----4-7----7--4-4------3--3---3----4-4---4-47-------77---4-----447-7--4-7-----377--7-7434373-"
                                  "-77--5412243422747--7-441244321337--777--3434333-77---777---------7------------7-----------7---77--4"
                                  "-7----477----4---7-----------4-774-4---7--742237-7-777-7---4---327447----4----43477423-----447-7--77"
                                  "21234--7---7477-77777-727777---7-3--4---77------7------7---4--4--4-4---7---7------------47727--47777"
                                  "4477-----7--4-7747-------47---7-7-7--------------------------74--7----------43-7--7477-774734-6----7"
                                  "47---6-42-447-774-4--------77-----7-34-------7-------747------77-7---------4------------4-----------"
                                  "7-----------------------4--4---7-.---7--------74437-------7-3424-732233427--7-72333127-6.--7777---6-"
                                  "-424347-------7-------47--------------77----77---------------7------------------77437-77------4-743-"
                                  "77-3-77---7----7-3377-------7-7-74--7-7---------7-7---4-------74---7----7------7-3447----7----------"
                                  "--------2443----7777-33--2-4444-7--44437---37---7746----34447-47---7----7--4--4----------4--4--7----"
                                  "44-----7---------7-4777-----3-77774-7------------7----47---------------------7--7-47---47-774777744-"
                                  "--34-7-44227----3...7712.1----123-7742-4-777774-7774----4-3--7--------------------77--447--744--77-7"
                                  "---------------3-...");

        // tests reported errors:
        TEST_EXPECT_ERROR_CONTAINS(PVP_calculate(gb_main, "ali_61s", "tree_nuc",  SAI_name), "alignment 'ali_61s' not found");
        TEST_EXPECT_ERROR_CONTAINS(PVP_calculate(gb_main, "ali_16s", "tree_miss", SAI_name), "tree not found");

        GB_close(gb_main);
    }

    { // protein data:
        GBDATA     *gb_main = GB_open("TEST_prot.arb", "rw");
        const char *aliname = "ali_tuf_pro";
        TEST_EXPECT_NO_ERROR(PVP_calculate(gb_main, aliname, "tree_prot", SAI_name));

        TEST_EXPECT_SAITYPE_EQUAL("PVP: Positional Variability by Parsimony: tree 'tree_prot' ntaxa 11");
        TEST_EXPECT_SAIDATA_EQUAL("-----4-634363-4-3-4666----6--6-6-6666466666664666466666-43-34366633334446-4666-364-6--66666664-44646"
                                  "666-4---4-66---6--66-4-66-6666444-6666466666--4--4--6466633346-62--3-3436343243266464433446433434662"
                                  "2-4634-663466644346314-23364666646243343-6366-34632-0632-46646-666-43-6-6---46-46-3-3432-32634464334"
                                  "443636666-66461612636----64664642421634-3-33444443632334-1-3-64-44664646634-433366666623624426246444"
                                  "4436330334446-63434344622344666320664664--666643--66-446341446--6------6--466----------");

        GB_close(gb_main);
    }
}

#endif // UNIT_TESTS

// --------------------------------------------------------------------------------