main.cc 6.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/* Spell-QTL  Software suite for the QTL analysis of modern datasets.
 * Copyright (C) 2016,2017  Damien Leroux <damien.leroux@inra.fr>, Sylvain Jasson <sylvain.jasson@inra.fr>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

Damien Leroux's avatar
Damien Leroux committed
18
#include "pedigree.h"
19
/*#include "bayes/factor_var4.h"*/
20
#include "bayes/output.h"
Damien Leroux's avatar
Damien Leroux committed
21
#include "cli.h"
22
#include "dispatch.h"
23
#include <fenv.h>
24
#include <input/marker_obs_formats.h>
Damien Leroux's avatar
Damien Leroux committed
25
26
27
28
29
30

void print_usage();


int main(int argc, const char** argv)
{
31
32
    feenableexcept(FE_ALL_EXCEPT & ~FE_INEXACT);

33
    msg_handler_t::set_color(true);
34
    msg_handler_t::debug_enabled() = false;
35
36
37
38
39
40
41
    bn_settings_t* settings = NULL;
    try {
         settings = bn_settings_t::from_args(argc, argv);
    } catch (file::error& fe) {
        MSG_ERROR("An error happened while reading input files. " << fe.what(), "");
        return -1;
    }
Damien Leroux's avatar
Damien Leroux committed
42
43
44
    if (!settings) {
        print_usage();
    } else {
45
46
        std::vector<std::string> cli(argv, argv + argc);
        MSG_INFO("Commandline: " << cli);
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
        settings->pedigree_filename = settings->work_directory + "/" + settings->pop_name + ".cache/" + settings->pop_name + ".spell-pedigree.data";
        if (!check_file(settings->pedigree_filename, false, false)) {
            MSG_ERROR("Couldn'f find the spell-pedigree output file.", "Please run spell-pedigree first with the same name and work_directory parameters.");
            return -1;
        }
        settings->pedigree.load(settings->pedigree_filename, false);
        MSG_DEBUG_INDENT_EXPR("[Pedigree] ");
        MSG_DEBUG("loaded from " << settings->pedigree_filename);
        MSG_DEBUG("" << settings->pedigree);
        MSG_DEBUG_DEDENT;

        for (const auto& kv: settings->observed_mark) {
            settings->marker_observation_specs[kv.second.format_name]
                    = marker_obs_formats::get_format(settings->pedigree, kv.second.format_name);
        }

63
64
65
66
67
68
        std::map<std::string, ObservationDomain> obs_gen;
        for (const auto& kv: settings->observed_mark) {
            obs_gen[kv.first] = settings->marker_observation_specs[kv.second.format_name].domain;
            MSG_DEBUG("OBSERVED GEN " << kv.first << ": " << obs_gen[kv.first]);
        }

69
        /*
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
        size_t n_alleles = 1;
        for (const auto& kv: settings->marker_observation_specs) {
            if (kv.second.domain == ODAllele) {
                if (n_alleles == 1) {
                    n_alleles = kv.second.domain_size;
                } else {
                    MSG_ERROR("Multiple and inconsistent formats (domain: allele)", "Consolidate your data and use at most ONE format per domain (alleles/ancestors) for your marker observations");
                }
            }
            settings->compiled_obs_specs[kv.first] = kv.second.compile();
        }
        for (const auto& kvfmt: settings->compiled_obs_specs) {
            MSG_DEBUG("COMPILED FORMAT " << kvfmt.first);
            for (const auto& kvscore: kvfmt.second) {
                MSG_DEBUG("  score '" << kvscore.first << "': " << kvscore.second.transpose());
            }
        }

        pedigree_bayesian_network bn = make_bn(pedigree, obs_gen, n_alleles, settings->noise, settings->tolerance);
        settings->bn = &bn;
90
        */
91
92
93
94
95
96
97
98

        std::set<std::string> tmp;
        for (const auto& obs: settings->observed_mark) {
            for (const auto& kv: obs.second.observations.data) {
                tmp.insert(kv.first);
            }
        }
        settings->marker_names.assign(tmp.begin(), tmp.end());
Damien Leroux's avatar
Damien Leroux committed
99
        /*MSG_DEBUG("MARKER NAMES " << settings->marker_names);*/
100
101
102
103
104
105

        for (size_t mn = 0; mn < settings->marker_names.size(); ++mn) {
            settings->marker_index[settings->marker_names[mn]] = mn;
        }

        if (settings->is_master()) {
106
            if (do_the_job(settings, "compute-alleles-per-marker")) {
107
108
109
110
111
                size_t n_mark = settings->count_markers();
                settings->alleles_per_marker.clear();
                settings->alleles_per_marker.resize(n_mark);
                std::set<size_t> unique_n_alleles;
                for (size_t mark = 0; mark < n_mark; ++mark) {
112
                    ifile result(settings->job_filename("compute-alleles-per-marker", mark));
113
114
115
116
117
118
119
120
                    /*auto& apm = settings->alleles_per_marker[mark];*/
                    std::set<char> temp;
                    rw_base()(result, temp);
                    unique_n_alleles.insert(temp.size());
                    /*MSG_DEBUG("Read an allele count of " << settings->alleles_per_marker[mark].size());*/
                    MSG_DEBUG("Read an allele count of " << temp.size());
                }
                settings->unique_n_alleles.assign(unique_n_alleles.begin(), unique_n_alleles.end());
121
                {
122
                    ofile unafs(settings->job_filename("unique_n_alleles", 0));
123
124
                    rw_base()(unafs, settings->unique_n_alleles);
                }
125
                MSG_DEBUG("Computed n_alleles: " << settings->unique_n_alleles);
Damien Leroux's avatar
Damien Leroux committed
126
127
                if (do_the_job(settings, "compute-factor-graphs")
                        && do_the_job(settings, "compute-LV")
128
129
                        && (settings->output_mode & bn_settings_t::OutputPopData)
                        && do_the_job(settings, "collect-LV")) {
130
                    ifile output(settings->job_filename("output", 0));
Damien Leroux's avatar
Damien Leroux committed
131
132
                    std::string filename;
                    rw_base() (output, filename);
133
                    settings->closing_messages.push_back(MESSAGE(GREEN << "Marker data for population `" << settings->pop_name << "' written in file " << WHITE << filename << NORMAL << std::endl));
Damien Leroux's avatar
Damien Leroux committed
134
                }
135
            }
136
137
138
            /*do_the_job(settings, "dummy-one");*/
            /*do_the_job(settings, "dummy");*/
            /*do_the_job(settings, "dummy-one");*/
139
            settings->cleanup_job_files();
140
            CREATE_MESSAGE(msg_channel::Out, MESSAGE(std::endl << std::endl));
141
142
143
144
            for (auto& msg: settings->closing_messages) {
                CREATE_MESSAGE(msg_channel::Out, std::move(msg));
            }
            MSG_QUEUE_FLUSH();
145
        } else {
146
            (void) do_the_job(settings, settings->job_name);
147
148
        }

Damien Leroux's avatar
Damien Leroux committed
149
150
151
152
153
154
        delete settings;
    }
    return 0;
}