convert_to_mapmaker.cc 3.18 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/>.
 */

18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43




void convert_to_mapmaker(const pop_data_type& popdata, const std::vector<chromosome>& gm, const std::string& gen_name, const marker_observation_spec& spec, char unknown, const std::string& data_type, const std::string& filename)
{
    std::vector<char> obs_symbols;
    const auto& labels = popdata.generations[gen_name]->unique_labels();
    MatrixXd obs_to_symbol = MatrixXd::Zero(spec.scores.size(), labels.size());

    obs_symbols.resize(spec.scores.size());

    for (size_t c = 0; c < spec.scores.size(); ++c) {
        std::vector<std::pair<char, char>> geno_vec;
        std::tie(obs_symbols[c], geno_vec) = scores[c];
        double coef = 1. / geno_vec.size();
        for (size_t g = 0; g < geno_vec.size(); ++g) {
            allele_pair l(geno_vec[g].first, geno_vec[g].second);
            for (size_t r = 0; r < obs_to_symbol.rows(); ++r) {
                if (labels[r] == l) {
                    obs_to_symbol(r, c) = coef;
                }
            }
        }
    }

44
    ofile output(filename);
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

    size_t n_ind = popdata.LV.data.begin()->second.find(gen_name)->second.size();
    size_t n_mark = 0;
    for (const auto& chrom: gm) {
        n_mark += chrom.raw.marker_names.size();
    }

    output << "data type " << data_type << std::endl;
    output << n_ind << ' ' << n_mark << " 0 0" << std::endl;
    
    for (const auto& chrom: gm) {
        size_t m = 0;
        for (const auto& mark: chrom.raw.marker_names) {
            output << '*' << mark << ' ';
            for (size_t i = 0; i < n_ind; ++i) {
                const auto& lv = popdata.LV(chrom.name, gen_name, i);
                const VectorXd obs = lv.col(i);
                VectorXd obs_like = obs_to_symbol * obs;
                int max_i = 0;
                double max = 0;
                bool fail = true;
                for (int i = 0; i < obs_like.size(); ++i) {
                    if (obs_like(i) > max) {
                        fail = false;
                        max_i = i;
                        max = obs_like(i);
                    } else if (obs_like(i) == max) {
                        fail = true;
                    }
                    if (fail) {
                        output << unknown;
                    } else {
                        output << obs_symbols[max_i];
                    }
                }
            }
            output << std::endl;
            ++m;
        }
    }
}