ld_matrices.cc 5.48 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
#include "error.h"
#include "data/chromosome.h"
20
/*#include "generation_rs.h"*/
21
#include "settings.h"
22
#include "output_impl.h"
23

24
void read_ld(settings_t* settings, const std::string& qtl_gen, const std::string& filename, std::istream& is)
25
26
27
28
29
30
{
    if (!settings->map.size()) {
        MSG_ERROR("Reading " << filename << ": No genetic map defined! Can't read LD data", "Define a genetic map before specifying the LD data file");
        return;
    }

31
32
33
34
    /*if (!settings->design) {*/
        /*MSG_ERROR("Reading " << filename << ": No breeding design defined! Can't read LD data", "Define a breeding design before specifying the LD data file");*/
        /*return;*/
    /*}*/
35
36
37
38
39
40
41
42
43
44
45
46
47
48

    std::string tmp;
    LD_matrices* cur = NULL;
    MatrixXd* curmat;
    int n_parents = -1;
    while (!is.eof()) {
        is >> tmp;
        if (tmp == "PARENTS") {
            is >> n_parents;
            if (n_parents <= 0) {
                MSG_ERROR("Reading " << filename << ": Invalid count of parents", "Verify the LD file wasn't altered. If you created the file by hand, please use clusthaplo instead.");
                return;
            }
            std::map<std::string, bool> ancestors;
49
50
51
52
            /* FIXME */
            /*for (const auto& an: settings->design->ancestor_names()) {*/
                /*ancestors[an] = false;*/
            /*}*/
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
            for (int i = 0; i < n_parents; ++i) {
                std::string tmp;
                is >> tmp;
                if (tmp == "CHROMOSOME" || tmp == "LOCUS") {
                    MSG_ERROR("Reading " << filename << ": Not enough parent names", "Verify the LD file wasn't altered. If you created the file by hand, please use clusthaplo instead.");
                    return;
                }
                auto it = ancestors.find(tmp);
                if (it == ancestors.end()) {
                    MSG_ERROR("Reading " << filename << ": Ancestor " << tmp << " isn't part of the breeding design.", "Verify you are using the correct LD file.");
                    return;
                }
                if (it->second) {
                    MSG_ERROR("Reading " << filename << ": Ancestor " << tmp << " appears more than once.", "Verify the LD file wasn't altered. If you created the file by hand, please use clusthaplo instead.");
                    return;
                }
                it->second = true;
70
                settings->ld_parents[qtl_gen].emplace_back(tmp);
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
            }
            for (const auto& kv: ancestors) {
                if (!kv.second) {
                    MSG_ERROR("Reading " << filename << ": Ancestor " << kv.first << " doesn't appear in the LD file.", "Verify you are using the correct LD file.");
                    return;
                }
            }
        } else if (tmp == "CHROMOSOME") {
            if (n_parents == -1) {
                MSG_ERROR("Reading " << filename << ": Expected PARENTS before CHROMOSOME", "Verify the LD file wasn't altered. If you created the file by hand, please use clusthaplo instead.");
                return;
            }
            is >> tmp;
            const chromosome* chr = settings->find_chromosome(tmp);
            if (!chr) {
                MSG_ERROR("Reading " << filename << ": Chromosome " << chr << " doesn't exist in the genetic map!", "Verify you are using the correct LD file.");
                return;
            }
89
            cur = &settings->ld_data[qtl_gen][chr];
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
        } else if (tmp == "LOCUS") {
            if (cur == NULL) {
                MSG_ERROR("Reading " << filename << ": Expected CHROMOSOME before LOCUS", "Verify the LD file wasn't altered. If you created the file by hand, please use clusthaplo instead.");
                return;
            }
            double l;
            int r, c;
            char v;
            is >> l >> r >> c;
            cur->loci.push_back(l);
            cur->ld.emplace_back(r, c);
            curmat = &cur->ld.back();
            for (int i = 0; i < r; ++i) {
                for (int j = 0; j < c; ++j) {
                    is >> v;
                    if (v != '0' && v != '1') {
                        MSG_ERROR("Reading " << filename << ": Illegal character '" << v << "' in matrix data", "Verify the LD file wasn't altered. If you created the file by hand, please use clusthaplo instead.");
                        return;
                    }
                    (*curmat)(i, j) = (double) (v == '1');
                }
            }
        } else {
            MSG_ERROR("Reading " << filename << ": Expected LOCUS or CHROMOSOME, read " << tmp << " instead", "Verify you are using the correct LD file and check it wasn't altered. If you created the file by hand, please use clusthaplo instead.");
            return;
        }
    }
}