frontends.cc 15.2 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
#include "computations/frontends4.h"
Damien Leroux's avatar
Damien Leroux committed
19
#include <unordered_set>
20
#include "output_impl.h"
Damien Leroux's avatar
Damien Leroux committed
21

Damien Leroux's avatar
Damien Leroux committed
22
#if 0
Damien Leroux's avatar
Damien Leroux committed
23
24
25
26
27
28
29
30
MatrixXd
ftest_along_chromosome(chromosome_value chr, const locus_key& lk,
                       const model& M0, const std::vector<double> pos)
{
    collection<model_block_type>
        popl = compute_parental_origins_multipop(all_populations(), chr, lk, pos);
    
    computation_along_chromosome cac;
Damien Leroux's avatar
Damien Leroux committed
31
    /*compute_along_chromosome(cac, FTest, Test, M0, M0, popl);*/
Damien Leroux's avatar
Damien Leroux committed
32
33
34
35
    return cac.ftest_pvalue;
}


36
37
38



Damien Leroux's avatar
Damien Leroux committed
39
40
41
42
std::pair<bool, double>
detect_strongest_qtl(chromosome_value chr, const locus_key& lk, double threshold,
                     const model& M0, const std::vector<double> pos)
{
Damien Leroux's avatar
Damien Leroux committed
43
    active_settings->set_title(MESSAGE("QTL detection along chromosome " << chr->name << " given " << lk));
Damien Leroux's avatar
Damien Leroux committed
44
    /*MatrixXd ftac = ftest_along_chromosome(chr, lk, M0, pos);*/
Damien Leroux's avatar
Damien Leroux committed
45
46
47
    int max = max_col(ftac);
    return {ftac(0, max) >= threshold, pos[max]};
}
Damien Leroux's avatar
Damien Leroux committed
48
#endif
Damien Leroux's avatar
Damien Leroux committed
49

Damien Leroux's avatar
Damien Leroux committed
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
87
88
89
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
#if 0
locus_selection
forward(model_manager& mm,
        const std::map<std::string, std::vector<std::string>>& markers)
{
    mm.init_loci_by_marker_list(markers);
    locus_selection cofactors;
    while (true) {
        auto result = mm.find_max_over_all_chromosomes();
        if (result.over_threshold) {
            /*cofactors.push_back({result.chrom, result.locus});*/
            cofactors.add_locus(result.chrom->name, result.locus);
            locus_key lk(new locus_key_struc({result.locus}));
            model_block_key k({{result.chrom, lk}});
            mm.Mcurrent.add_block(k, result.block);
        }
    }
    return cofactors;
}

model_manager&
manual_skeleton(model_manager& mm,
                const std::map<std::string, std::vector<std::string>>& markers)
{
    mm.init_loci_by_marker_list(markers);
    return mm;
}


model_manager&
init_skeleton(model_manager& mm)
{
    std::map<std::string, std::vector<double>> skeleton;
    if (active_settings->skeleton_mode == "auto") {
        for (const chromosome& chr: active_settings->map) {
            double accept = -1.;
            for (double l: chr.condensed.marker_locus) {
                if (l > accept) {
                    skeleton[chr.name].push_back(l);
                    accept = l + active_settings->skeleton_interval;
                }
            }
        }
    } else if (active_settings->skeleton_mode == "manual") {
        for (const auto& name: active_settings->skeleton_markers) {
            bool found = false;
            for (const chromosome& chr: active_settings->map) {
                auto li = chr.condensed.marker_locus.begin();
                auto lj = chr.condensed.marker_locus.end();
                auto ni = chr.condensed.marker_name.begin();
                for (; li != lj; ++li, ++ni) {
                    if (*ni == name) {
                        skeleton[chr.name].push_back(*li);
                        found = true;
                        break;
                    }
                }
                if (found) {
                    break;
                }
            }
        }
    }
    mm.init_loci_by_hand(skeleton);

    return mm;
}
118
#endif
119
120
121
122
123
124
125



model_manager&
forward(model_manager& mm)
{
    /*MSG_DEBUG("FORWARD");*/
126
127
    while (true) {
        const auto& result = mm.search_new_best();
128
129
        /*MSG_DEBUG(result);*/
        if (result.over_threshold) {
130
            result.select(mm);
131
            MSG_INFO("Select " << result.chrom->name << ':' << result.locus << " as cofactor");
132
133
        } else {
            break;
134
        }
135
136
    }
    MSG_INFO("Cofactors: " << mm.get_selection());
137
138
    return mm;
}
Damien Leroux's avatar
Damien Leroux committed
139

Damien Leroux's avatar
Damien Leroux committed
140
141
142
143

model_manager&
backward(model_manager& mm)
{
144
145
146
147
148
    mm.select_all_loci();
    while (true) {
        auto worst_kb = mm.weakest_link();
        if (worst_kb.second >= mm.threshold) {
            break;
Damien Leroux's avatar
Damien Leroux committed
149
        }
150
151
152
153
        const model_block_key& mbk = worst_kb.first.first;
        assert(mbk->type == mbk_POP);
        for (double d: mbk->loci) {
            mm.deselect(mbk->chr, d);
Damien Leroux's avatar
Damien Leroux committed
154
        }
155
156
        /*MSG_DEBUG("Removed " << result.block_key << " which scored only " << result.test_value << " and now model has " << mm.Mcurrent.m_blocks.size() << " blocks");*/
    }
Damien Leroux's avatar
Damien Leroux committed
157
158
159
160
161
162
163
164
    return mm;
}



model_manager&
qtl_detect_cim(model_manager& mm)
{
165

Damien Leroux's avatar
Damien Leroux committed
166
167
168
169
170
171
172
    return mm;
}


model_manager&
qtl_detect_cim_minus(model_manager& mm)
{
173
174
175
176
177
    auto chrom_best = mm.search_new_best_per_chromosome(true);
    mm.clear_selection();
    for (const auto& cb: chrom_best) {
        if (cb.second.over_threshold) {
            cb.second.select(mm);
Damien Leroux's avatar
Damien Leroux committed
178
        }
179
    }
Damien Leroux's avatar
Damien Leroux committed
180
181
182
183
    return mm;
}


184
185
struct loop_exception {};

186
187
188

bool
iqtlm_backward_one(model_manager& mm)
189
{
190
191
192
193
194
195
196
197
198
199
    std::vector<std::pair<bool, test_result>> results;
    size_t sz = 0;

    for (const auto& chr_sel: mm.get_selection()) { sz += chr_sel.second->depth(); }
    results.reserve(sz);

    bool not_modified = true;
    
    for (const auto& chr_sel: mm.get_selection()) {
        for (double d: chr_sel.second) {
200
201
202
            std::pair<bool, test_result> result = mm.challenge_qtl(chr_sel.first, d);
            MSG_DEBUG("Challenging QTL at " << chr_sel.first << ':' << d << " gives " << result.first << ' ' << result.second);
            results.emplace_back(result);
203
            not_modified &= results.back().first;
204
205
        }
    }
206

207
208
209
210
211
212
213
214
215
    {
        std::stringstream dbg;
        dbg << "backward";
        for (const auto& bt: results) {
            dbg << " {" << (bt.first ? "true" : "false") << ' ' << bt.second << '}';
        }
        MSG_DEBUG(dbg.str());
    }

216
217
218
219
220
221
222
223
224
    if (not_modified) {
        return false;
    }

    mm.clear_selection();

    for (const auto& good_result: results) {
        if (good_result.second.over_threshold) {
            good_result.second.select(mm);
225
226
        }
    }
227
228
229
230

    return true;
}

231
232

bool
233
in_history(std::set<std::map<chromosome_value, locus_key>>& history, const std::map<chromosome_value, locus_key>& selection)
234
{
235
    return !history.insert(selection).second;
236
237
238
}


239
240
bool
iqtlm_backward(model_manager& mm, std::set<std::map<chromosome_value, locus_key>>& history)
241
{
242
    DUMP_FILE_LINE();
243
    while (iqtlm_backward_one(mm)) {
244
    DUMP_FILE_LINE();
245
        if (in_history(history, mm.get_selection())) {
246
    DUMP_FILE_LINE();
247
248
249
            return false;
        }
    }
250
    DUMP_FILE_LINE();
251
    return true;
252
253
254
}


255
256
bool
iqtlm_forward(model_manager& mm, std::set<std::map<chromosome_value, locus_key>>& history)
257
{
258
    DUMP_FILE_LINE();
259
    while (true) {
260
    DUMP_FILE_LINE();
261
262
263
        auto results = mm.search_new_best_per_chromosome();
        bool modified = false;
        for (const auto& result: results) {
264
    DUMP_FILE_LINE();
265
266
            modified |= result.second.over_threshold;
            if (result.second.over_threshold) {
267
    DUMP_FILE_LINE();
268
                result.second.select(mm);
269
                MSG_DEBUG("added " << result);
270
271
            }
        }
272
    DUMP_FILE_LINE();
273
        if (!modified) {
274
    DUMP_FILE_LINE();
275
276
277
            return true;
        }
        if (in_history(history, mm.get_selection())) {
278
    DUMP_FILE_LINE();
279
280
            return false;
        }
281
282
283
284
    }
}


285
286
bool
iqtlm_forward_gw(model_manager& mm, std::set<std::map<chromosome_value, locus_key>>& history)
287
{
288
289
290
291
292
293
294
295
296
297
298
299
    while (true) {
        auto result = mm.search_new_best();
        bool modified = result.over_threshold;
        if (result.over_threshold) {
            result.select(mm);
        }
        if (!modified) {
            return true;
        }
        if (in_history(history, mm.get_selection())) {
            return false;
        }
300
    }
301
}
302
303


304
305
306
model_manager&
qtl_detect_iqtlm(model_manager& mm)
{
307
    MSG_INFO("iQTLm starts with: " << mm.get_selection());
308
309
310
    std::set<std::map<chromosome_value, locus_key>> history;
    std::map<chromosome_value, locus_key> last_sel, new_sel;
    new_sel = mm.get_selection();
311
    do {
312
313
314
315
        if (iqtlm_backward(mm, history) && iqtlm_forward(mm, history)) {
            last_sel = new_sel;
            new_sel = mm.get_selection();
        } else {
316
317
            break;
        }
318
    } while (new_sel != last_sel);
319
320
    return mm;
}
321
322


323
324
325
326
327
328
329
330
331
332
333
334
335
336
model_manager&
qtl_detect_iqtlm_gw(model_manager& mm)
{
    std::set<std::map<chromosome_value, locus_key>> history;
    std::map<chromosome_value, locus_key> last_sel, new_sel;
    new_sel = mm.get_selection();
    while (new_sel != last_sel) {
        if (iqtlm_backward(mm, history) && iqtlm_forward_gw(mm, history)) {
            last_sel = new_sel;
            new_sel = mm.get_selection();
        } else {
            break;
        }
    }
337
338
339
    return mm;
}

340
341

#if 0
Damien Leroux's avatar
Damien Leroux committed
342
343
344
345
346
347
model_manager&
qtl_detect_iqtlm(model_manager& mm)
{
    bool go_on = true;
    model_manager::test_result result;
    mm.init_loci_by_step(active_settings->step);
348
    /*mm.set_joint_mode();*/
349
    /*mm.pmode = Joint;*/
Damien Leroux's avatar
Damien Leroux committed
350
    std::map<const chromosome*, std::unordered_set<locus_key>> previous_selections;
351
    std::map<const chromosome*, std::vector<std::pair<locus_key, MatrixXd>>> ordered_selections;
352
    std::map<const chromosome*, locus_key> previous;
Damien Leroux's avatar
Damien Leroux committed
353

354
355
356
357
    /*mm.output_test = true;*/
    /*mm.output_rank = true;*/
    /*mm.output_rss = true;*/
    /*mm.output_model = true;*/
Damien Leroux's avatar
Damien Leroux committed
358

359
360
    auto check_loops = [&] (const chromosome* chr)
    {
Damien Leroux's avatar
Damien Leroux committed
361
362
363
        auto& history = previous_selections[chr];
        locus_key current = mm.get_selection(chr);

364
365
366
367
368
369
        if (previous[chr] == current) {
            MSG_INFO("Selection is stable. Proceeding.");
            return false;
        }
        previous[chr] = current;

Damien Leroux's avatar
Damien Leroux committed
370
        MSG_INFO("[iQTLm] pass ended with selection " << mm.keys() << " (chromosome selection is " << current << ')');
371
        /*MSG_INFO("[iQTLm] sizeof(history) = " << history.size());*/
Damien Leroux's avatar
Damien Leroux committed
372
373
374
375
376
377
378
379
380
381
382
383
384
385
        auto i = history.begin(), j = history.end();
        if (i != j) {
            std::stringstream s;
            s << (*i);
            for (++i; i != j; ++i) {
                s << ", " << (*i);
            }
            MSG_DEBUG("History: " << s.str());
        }
        if (current->is_empty()) {
            return false;
        } else {
            auto ret = history.insert(current);
            if (!ret.second) {
386
387
                /*return false;*/
                throw loop_exception();
Damien Leroux's avatar
Damien Leroux committed
388
            }
389
            ordered_selections[chr].push_back({current, mm.Mcurrent.rss() / (mm.Mcurrent.X().innerSize() - mm.Mcurrent.rank())});
Damien Leroux's avatar
Damien Leroux committed
390
391
392
393
        }
        return true;
    };

394
395
396
397
398
399
400
401
402
403
404
405
406
407
    auto fix_loop = [&] (const chromosome* chr)
    {
        auto& history = ordered_selections[chr];
        locus_key current = mm.get_selection(chr);
        int i;
        for (i = history.size() - 1;
             i >= 0 && history[i].first != current;
             --i);
        if (i < 0) {
            MSG_ERROR("Couldn't find selection in history.", "");
        }
        int i_best = i;
        const int i_max = (int) history.size();
        for (; i < i_max; ++i) {
408
            if (history[i].second.lpNorm<1>() < history[i_best].second.lpNorm<1>()) {
409
410
411
412
413
                i_best = i;
            }
        }
        mm.set_selection(history[i_best].first);
    };
Damien Leroux's avatar
Damien Leroux committed
414

415
    for (const chromosome& chr: active_settings->map) {
416
        active_settings->set_title(MESSAGE("iQTLm on " << chr.name));
417
418
        try {
            /* setup */
Damien Leroux's avatar
Damien Leroux committed
419

420
421
            mm.select_chromosome(&chr);
            mm.cofactors_to_QTLs();
422
            previous[&chr] = mm.get_selection(&chr);
Damien Leroux's avatar
Damien Leroux committed
423
            do {
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
                /* step 1: backward */

                bool modified;

                MSG_INFO("[iQTLm] Backward");

                do {
                    /*locus_key loci = mm.Mcurrent.blocks_with_chromosome(&chr).begin()->first[&chr];*/
                    locus_key loci = mm.get_selection(&chr);
                    locus_key tmp_sel;

                    modified = false;

                    for (double l: *loci) {
                        auto tmp_removed = mm.remove(l);
                        /*mm.add_test_locus(&chr, l);*/
                        active_settings->set_title(MESSAGE("iQTLm on " << chr.name << " checking locus " << l << " against " << mm.keys()));
                        result = mm.test_along_chromosome(0, mm.max_testpos());
                        MSG_DEBUG(result);
                        if (result.over_threshold) {
                            model_block_key k1, k2;
                            k1 += std::make_pair(&chr, l);
                            k2 += std::make_pair(&chr, result.locus);
                            /*mm.add(result);*/
                            tmp_sel = tmp_sel + result.locus;
449
450
451
452
453
                            if (k1 != k2) {
                                MSG_DEBUG("[iQTLm] Replaced " << k1 << " with " << k2);
                            } else {
                                MSG_DEBUG("[iQTLm] Kept " << k1);
                            }
454
455
456
457
458
459
460
461
                            modified |= (result.locus != l);
                        } else {
                            model_block_key k1, k2;
                            k1 += std::make_pair(&chr, l);
                            MSG_DEBUG("[iQTLm] Removed " << k1);
                            modified = true;
                        }
                        mm.set_selection(tmp_removed);
Damien Leroux's avatar
Damien Leroux committed
462
                    }
463
464
465
466
467
                    if (modified) {
                        mm.set_selection(tmp_sel);
                    }
                } while (modified && check_loops(&chr));

468
                go_on = modified;
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483

                /* step 2: forward */

                if (!check_loops(&chr)) {
                    active_settings->set_title(MESSAGE("iQTLm on " << chr.name << " forward"));
                    MSG_INFO("[iQTLm] Forward");
                    do {
                        result = mm.test_along_chromosome(0, mm.max_testpos());
                        MSG_DEBUG(result);
                        if (result.over_threshold) {
                            model_block_key k1;
                            k1 += std::make_pair(&chr, result.locus);
                            /*MSG_INFO("[iQTLm] Added " << (model_block_key() + std::make_pair(&chr, result.locus)));*/
                            mm.add(result);
                            MSG_DEBUG("[iQTLm] Added " << k1 << ", current selection is " << mm.keys());
484
                            go_on = true;
485
486
                        }
                    } while (result.over_threshold);
487
                    MSG_INFO("[iQTLm] pass ended with selection " << mm.keys() << " (chromosome selection is " << mm.get_selection(&chr) << ')');
Damien Leroux's avatar
Damien Leroux committed
488
489
                }

490
            } while (go_on && check_loops(&chr));
Damien Leroux's avatar
Damien Leroux committed
491

492
493
494
495
496
497
            /* teardown */
            /*mm.QTLs_to_cofactors();*/
        } catch (const loop_exception&) {
            fix_loop(&chr);
            MSG_WARNING("[iQTLm] Cycle detected while studying chromosome " << chr.name << ". Final selection is " << mm.get_selection(&chr));
        }
Damien Leroux's avatar
Damien Leroux committed
498
499
500
501
    }
    return mm;
}

502
#endif
Damien Leroux's avatar
Damien Leroux committed
503