frontends2.h 40.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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
44
45
46
47
48
49
50
51
52
53
54
55
#ifndef _SPEL_FRONTENDS_H_
#define _SPEL_FRONTENDS_H_

#include "cache2.h"
#include "basic_data.h"
#include "model.h"

std::pair<bool, double>
detect_strongest_qtl(chromosome_value chr, const locus_key& lk,
                     const model& M0, const std::vector<double> pos);

MatrixXd
ftest_along_chromosome(chromosome_value chr, const locus_key& lk,
                       const model& M0, const std::vector<double> pos);


/* Definitions:
 * - cofactor: isolated POP (single chromosome, single locus)
 * - QTLs: joint POP (single chromosome, single or multiple loc(i)us)
 *
 *
 * Configurations:
 * - with/without Dominance                          D
 * - with/without Constraints (can/can't estimate)   WC
 * - Joint/Single POP computation mode               JS
 *
 *
 * Steps:                           D   JS  WC
 * - establish skeleton                      
 *   - manual (marker list)                  
 *   - by step                               
 * - discover cofactors                 S    
 *   - manual                           S    
 *   - forward                          S    
 *   - backward                         S    
 * - detect QTLs                    ?   J   Y
 *   - CIM-                         ?   J   Y
 *   - iQTLm                        ?   J   Y
 *   - iQTLm++                      ?   J   Y
 * - OPTIONALLY analyze epistasis   ?   J   Y
 * - estimate parameters            
 *
 *
 * Operations:
 * - select chromosome
 * - cofactors to QTLs for the current chromosome
 * - QTLs to cofactors for the current chromosome
 * - test along the chromosome
 * - test along all chromosomes
 * - add cofactor (if current chromosome in product probability mode)
 * - add QTL (if current chromosome in joint probability mode)
 * - remove cofactor/QTL
 */


Damien Leroux's avatar
WIP.    
Damien Leroux committed
56
struct signal_display {
57
#ifdef SIGNAL_DISPLAY_ONELINER
Damien Leroux's avatar
WIP.    
Damien Leroux committed
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
    static const char* tick(double x)
    {
        static const char* ticks[9] = { " ", "\u2581", "\u2582", "\u2583", "\u2584", "\u2585", "\u2586", "\u2587", "\u2588" };
        return ticks[x < 0. ? 0
                            : x >= 1. ? 8
                                      : ((int) floor(x * 9))];
    }

    VectorXd values;
    int imax_;
    bool above_;

    signal_display(const VectorXd& v, int imax, bool above)
        : values(v.innerSize()), imax_(imax), above_(above)
    {
        values = v;
74
75
76
77
#if 0
        int sig_cols = msg_handler_t::termcols() - 3;
        MSG_DEBUG("values.innerSize = " << values.innerSize());
        MSG_QUEUE_FLUSH();
Damien Leroux's avatar
WIP.    
Damien Leroux committed
78
        while (values.innerSize() >= sig_cols) {
79
80
81
82
83
            if (values.innerSize() & 1) {
                int sz = values.innerSize();
                values.conservativeResize(sz + 1);
                values(sz) = values(sz - 1);
            }
Damien Leroux's avatar
WIP.    
Damien Leroux committed
84
85
            int i = values.innerSize() >> 1;
            values = values.transpose() * kroneckerProduct(MatrixXd::Identity(i, i), MatrixXd::Constant(1, 2, .5));
86
87
            MSG_DEBUG("values.innerSize = " << values.innerSize());
            MSG_QUEUE_FLUSH();
Damien Leroux's avatar
WIP.    
Damien Leroux committed
88
        }
89
#endif
Damien Leroux's avatar
WIP.    
Damien Leroux committed
90
91
        double vmin = values.minCoeff();
        double vmax = values.maxCoeff();
92
93
94
95
96
        if (vmin == vmax) {
            values = (values.array() - vmin).matrix();
        } else {
            values = ((values.array() - vmin) / (vmax - vmin)).matrix();
        }
Damien Leroux's avatar
WIP.    
Damien Leroux committed
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
    }

    friend std::ostream& operator << (std::ostream& os, const signal_display& sd)
    {
        os << _WHITE << '[';
        for (int i = 0; i < sd.values.innerSize(); ++i) {
            if (i == sd.imax_) {
                os << (sd.above_ ? _GREEN : _RED);
            }
            os << tick(sd.values(i));
            if (i == sd.imax_) {
                os << _WHITE;
            }
        }
        return os << ']' << _NORMAL;
    }
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
#else
    braille_grid grid;

    signal_display(const chromosome& chr, const std::vector<double>& X, const VectorXd& y, int imax, double threshold)
        : grid(build(chr, X, y, imax, threshold))
    {}

    braille_grid
        build(const chromosome& chr, const std::vector<double>& X, const VectorXd& y, int imax, double threshold)
        {
            std::vector<double> Y(y.data(), y.data() + y.size());
            int padding_left = 0;
            int W = (int) (msg_handler_t::termcols() * .8);
            braille_grid chr_map = chr.pretty_print(W, {}, {}, padding_left, false);

            braille_plot plot(W - padding_left, 5, 0, X.back(), 0, std::max(threshold, y(imax)));
            plot.plot(X, Y);
            plot.hline(threshold, 1, 1, 0, 255, 0);
            bool above = y(imax) > threshold;
            plot.vline(X[imax], 1, 0, above ? 0 : 255, above ? 255 : 0, 0);
            return plot.compose_vert(true, chr_map, false);
        }

    friend
        std::ostream&
        operator << (std::ostream& os, const signal_display& sd)
        {
            return os << sd.grid;
        }
#endif
Damien Leroux's avatar
WIP.    
Damien Leroux committed
143
144
145
};


146
147
148
149

enum probability_mode { Joint, Single };

struct model_manager {
150
151
    /* name of the studied trait */
    std::string trait_name;
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
    /* populations used in this model */
    collection<population_value> all_pops;
    /* POP matrices per locus per population (in the order of active_settings->populations) */
    std::vector<collection<parental_origin_per_locus_type>> pop_blocks;
    /* number of parents (for locus_key reduction of POP matrices) per population (same order) */
    std::vector<size_t> n_par;
    /* thy reference model */
    model Mcurrent;
    /* thy cofactor/qtl acceptance threshold */
    double threshold;
    /* the test to perform (FTest, Chi2...) */
    ComputationType test_type;
    /* the lists of loci behind the expression "compute X along the chromosome" per chromosome */
    std::map<const chromosome*, value<std::vector<double>>> all_loci;
    /* the chrrently selected chromosome */
    const chromosome* chrom_under_study;
    /* the list of loci for the chromosome currently under study */
    value<std::vector<double>> loci;
    /* the model blocks along the chromosome currently under study */
    collection<model_block_type> locus_blocks;
    /* the keys for each model block in locus_blocks */
    /*std::vector<model_block_key> locus_keys;*/
    model_block_key locus_base_key;
    /* structure to cache the results of the last computation along the chromosome */
    computation_along_chromosome cac;
    /* a direct access to the result of the last computation (&cac.TEST_TYPE) */
    MatrixXd* last_computation;
    /* the current probability mode */
    probability_mode pmode;

Damien Leroux's avatar
Damien Leroux committed
182
183
184
185
186
    bool output_rss;
    bool output_rank;
    bool output_test;
    bool output_model;

187
188
189
    /* slight hack. */
    std::vector<double> testpos;

190
191
192
    /* additional data for output and postprocessing */
    std::map<std::string, std::map<double, std::pair<double, double>>> qtl_confidence_interval;

193
194
/* Construction
 */
195
    model_manager(const std::string& trait, const collection<population_value>& colpops,
196
197
198
                  const value<MatrixXd>& y, double th,
                  ComputationType ct = ComputationType::FTest,
                  SolverType st = SolverType::QR)
199
200
        : trait_name(trait)
        , all_pops(colpops)
201
202
        , pop_blocks()
        , n_par()
Damien Leroux's avatar
Damien Leroux committed
203
        , Mcurrent(y, colpops, st)
204
205
206
207
208
209
210
211
212
        , threshold(th)
        , test_type(ct)
        , all_loci()
        , chrom_under_study(NULL)
        , loci()
        , locus_blocks()
        , cac()
        , last_computation(NULL)
        , pmode(Single)
Damien Leroux's avatar
Damien Leroux committed
213
214
215
216
        , output_rss(false)
        , output_rank(false)
        , output_test(false)
        , output_model(false)
217
        , testpos()
218
        , qtl_confidence_interval()
219
220
221
    {
        Mcurrent.add_block({}, cross_indicator(all_pops));
        Mcurrent.compute();
222
223
        MSG_DEBUG("INITIAL MODEL SIZE: " << MATRIX_SIZE(Mcurrent.X()));
        MSG_QUEUE_FLUSH();
224
225
226
227
228
229
230
231
232
233
234
235
236
        for (const auto& kpop: active_settings->populations) {
            context_key ck(new context_key_struc(&kpop.second,
                                                 &active_settings->map[0],
                                                 std::vector<double>()));
            n_par.push_back(
                    make_value<Mem>(compute_state_to_parental_origin,
                                    as_value(ck),
                                    as_value(locus_key()))->innerSize());
        }
        pop_blocks.resize(active_settings->populations.size());
    }


237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
    model_manager(const model_manager& mm)
        : all_pops(mm.all_pops)
        , pop_blocks(mm.pop_blocks)
        , n_par(mm.n_par)
        , Mcurrent(mm.Mcurrent)
        , threshold(mm.threshold)
        , test_type(mm.test_type)
        , all_loci(mm.all_loci)
        , chrom_under_study(mm.chrom_under_study)
        , loci(mm.loci)
        , locus_blocks(mm.locus_blocks)
        , cac()
        , last_computation(NULL)
        , pmode(Single)
        , output_rss(mm.output_rss)
        , output_rank(mm.output_rank)
        , output_test(mm.output_test)
        , output_model(mm.output_model)
        , testpos(mm.testpos)
    {}


259
260
261
262
263
264
265
266
/* Initializing test positions (TODO: optional exclusion window)
 */

    void init_loci_by_step(double step)
    {
        for (const chromosome& c: active_settings->map) {
            all_loci[&c]
                = compute_steps(
Damien Leroux's avatar
Damien Leroux committed
267
                    c.condensed.marker_locus,
268
                    step);
269
            /*MSG_DEBUG("test loci for chromosome " << c << std::endl << all_loci[&c]);*/
270
        }
Damien Leroux's avatar
Damien Leroux committed
271
        pop_blocks.clear();
272
273
274
275
276
277
278
279
280
281
282
    }

    void init_loci_by_hand(const std::map<std::string, std::vector<double>>& locmap)
    {
        for (auto mark_list: locmap) {
            const chromosome* c = active_settings->find_chromosome(mark_list.first);
            value<std::vector<double>>& vml = all_loci[c];
            vml = std::vector<double>();
            *vml = mark_list.second;
            std::sort(vml->begin(), vml->end());
        }
Damien Leroux's avatar
Damien Leroux committed
283
        pop_blocks.clear();
284
285
286
287
288
289
290
291
292
293
294
295
296
297
    }

    void init_loci_by_marker_list(const std::map<std::string, std::vector<std::string>>& markers)
    {
        for (auto mark_list: markers) {
            const chromosome* c = active_settings->find_chromosome(mark_list.first);
            value<std::vector<double>>& vml = all_loci[c];
            vml = std::vector<double>();
            vml->reserve(mark_list.second.size());
            for (const auto& mn: mark_list.second) {
                vml->push_back(c->raw.locus_by_name(mn));
            }
            std::sort(vml->begin(), vml->end());
        }
Damien Leroux's avatar
Damien Leroux committed
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
        pop_blocks.clear();
    }

    void init_loci_by_step_with_exclusion(double step, double exclusion_window = 0)
    {
        /* exclude all loci already in the model block keys */
        std::map<const chromosome*, std::vector<std::pair<double, double>>> excl;
        for (const auto& kb: Mcurrent.m_blocks) {
            for (const auto& ch_loc: kb.first.selection) {
                auto& exvec = excl[ch_loc.first];
                for (double l: ch_loc.second->to_vec()) {
                    exvec.emplace_back(l - exclusion_window, l + exclusion_window);
                }
            }
        }
        /* sort all exclusion windows */
        for (auto& cv: excl) {
            std::sort(cv.second.begin(), cv.second.end(),
                      [] (const std::pair<double, double>& a, const std::pair<double, double>& b)
                      {
                          return a.first < b.first;
                      });
        }
        /* initialize all_loci */
        for (const chromosome& c: active_settings->map) {
            all_loci[&c]
                = compute_steps(
                    c.condensed.marker_locus,
                    step,
                    excl[&c]);
        }
        pop_blocks.clear();
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
    }

/* Chromosome selection
 */
    void select_chromosome(const std::string& name)
    {
        select_chromosome(active_settings->find_chromosome(name));
    }

    void select_chromosome(const chromosome* c)
    {
        if (chrom_under_study != c) {
            chrom_under_study = c;
            pop_blocks.clear();
        }
Damien Leroux's avatar
Damien Leroux committed
345
        loci = all_loci[chrom_under_study];
346
347
348
349
350
351
352
    }

    const chromosome* selected_chromosome() const
    {
        return chrom_under_study;
    }

353
354
355
356
357
358
359
360
361
362
363
364
365
366
    void
        set_joint_mode()
        {
            MSG_DEBUG("[Model Manager] JOINT MODE");
            pmode = Joint;
        }

    void
        set_product_mode()
        {
            MSG_DEBUG("[Model Manager] PRODUCT MODE");
            pmode = Single;
        }

367
368
/* cofactors <-> QTLs for selected chromosome
 */
Damien Leroux's avatar
Damien Leroux committed
369
    locus_key cofactors_to_QTLs()
370
    {
371
        set_joint_mode();
372
373
        model_block_collection mbc = Mcurrent.extract_blocks_with_chromosome(chrom_under_study);
        model_block_key qtl_key = model_block_key::join(mbc);
Damien Leroux's avatar
Damien Leroux committed
374
375
376
        if (qtl_key.empty()) {
            return {};
        }
377
378
        const locus_key& lk = qtl_key[chrom_under_study];
        std::vector<double> loc{lk->locus};
Damien Leroux's avatar
Damien Leroux committed
379
380
381
        /*MSG_DEBUG(__LINE__ <<":chromosome " << chrom_under_study->name << " loci " << (*loci)); MSG_QUEUE_FLUSH();*/
        std::vector<double> allpos = lk->to_vec();
        std::vector<double> testpos = {lk->locus};
382
383
384
385
386
387
        Mcurrent.add_block(
            qtl_key,
            compute_parental_origins_multipop(
                all_pops,
                as_value(chrom_under_study),
                as_value(lk->parent),
Damien Leroux's avatar
Damien Leroux committed
388
389
                as_value(allpos),
                testpos)[0]);
390
        Mcurrent.compute();
Damien Leroux's avatar
Damien Leroux committed
391
392
        locus_base_key = qtl_key;
        return qtl_key.selection.begin()->second;
393
394
395
396
    }

    void QTLs_to_cofactors()
    {
397
        set_product_mode();
398
399
        model_block_collection mbc = Mcurrent.extract_blocks_with_chromosome(chrom_under_study);
        /* There shall be only ONE key in mbc. */
Damien Leroux's avatar
Damien Leroux committed
400
401
402
        if (mbc.size() == 0) {
            return;
        }
403
404
405
        std::vector<model_block_key> qtl_keys = mbc.begin()->first.split_loci();
        for (const auto& mbk: qtl_keys) {
            const locus_key& lk = mbk[chrom_under_study];
Damien Leroux's avatar
Damien Leroux committed
406
407
408
409
            std::vector<double> loc = lk->to_vec();
            loc.push_back(lk->locus);
            std::sort(loc.begin(), loc.end());
            /*MSG_DEBUG(__LINE__ <<":chromosome " << chrom_under_study->name << " loci " << loc); MSG_QUEUE_FLUSH();*/
410
411
412
413
414
415
            Mcurrent.add_block(
                mbk,
                compute_parental_origins_multipop(
                    all_pops,
                    as_value(chrom_under_study),
                    as_value(lk->parent),
Damien Leroux's avatar
Damien Leroux committed
416
417
                    as_value(loc),
                    loc)[0]);
418
419
            Mcurrent.compute();
        }
Damien Leroux's avatar
Damien Leroux committed
420
        locus_base_key.selection.clear();
421
422
423
424
425
426
427
428
429
430
431
432
433
434
    }

/* Test along the currently selected chromosome
 * TODO
 */
    struct test_result {
        const chromosome* chrom;
        double locus;
        double test_value;
        int index;
        bool over_threshold;
        model_block_key block_key;
        value<model_block_type> block;

Damien Leroux's avatar
Damien Leroux committed
435
436
437
438
439
440
441
442
        test_result()
            : chrom(NULL), locus(0), test_value(0), index(0), over_threshold(false), block_key(), block()
        {}

        test_result(const chromosome* c, double l, double tv, int i, bool ot, const model_block_key& mbk, const value<model_block_type>& mb)
            : chrom(c), locus(l), test_value(tv), index(i), over_threshold(ot), block_key(mbk), block(mb)
        {}

443
444
445
446
447
448
        test_result(const test_result& tr)
            : chrom(tr.chrom), locus(tr.locus), test_value(tr.test_value),
              index(tr.index), over_threshold(tr.over_threshold),
              block_key(tr.block_key), block(tr.block)
        {}

Damien Leroux's avatar
Damien Leroux committed
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
        test_result&
            operator = (const test_result& tr)
            {
                chrom = tr.chrom;
                locus = tr.locus;
                test_value = tr.test_value;
                index = tr.index;
                over_threshold = tr.over_threshold;
                block_key = tr.block_key;
                block = tr.block;
                return *this;
            }

        void reset()
        {
                chrom = NULL;
                locus = 0;
                test_value = 0;
                index = 0;
                over_threshold = false;
                block_key.selection.clear();
                block = value<model_block_type>();
        }

473
474
475
476
477
478
479
480
481
482
483
484
485
486
        friend
            std::ostream& operator << (std::ostream& os, const test_result& tr)
            {
                os << "<result chrom=" << (tr.chrom ? tr.chrom->name : "nil")
                    << " locus=" << tr.locus
                    << " test=" << tr.test_value
                    << " at=" << tr.index
                    << " over?=" << tr.over_threshold
                    << " block_key=" << tr.block_key
                    << '>';
                return os;
            }
    };

487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
    /* todo: study_maha(locus)
     * creates base model by removing locus from proper block and computes maha distance cleanly
     */

    MatrixXd study_maha(const chromosome* chr, double locus, double minpos, double maxpos)
    {
        model_manager tmp(*this);
        tmp.select_chromosome(chr);
        tmp.remove(locus);
        value<model> Mbase = tmp.Mcurrent;
        return custom_test_along_chromosome(ComputationType::Mahalanobis, ComputationResults::Test, Mbase, minpos, maxpos).mahalanobis;
    }

/*enum ComputationType { NoTest=0, FTest=1, FTestLOD=2, R2=4, Chi2=8, Mahalanobis=16 };*/
    std::string
        computation_type_to_string(ComputationType ct)
        {
            std::stringstream ret;
            if (ct & FTest) { ret << "_FTest"; }
            if (ct & FTestLOD) { ret << "_FTestLOD"; }
            if (ct & R2) { ret << "_R2"; }
            if (ct & Chi2) { ret << "_Chi2"; }
            if (ct & Mahalanobis) { ret << "_Mahalanobis"; }
            return ret.str();
        }

    double max_testpos() const { return max_testpos(chrom_under_study); }
    double max_testpos(const chromosome* c) const { return all_loci.find(c)->second->back(); }

Damien Leroux's avatar
Damien Leroux committed
516
    const computation_along_chromosome&
517
        custom_test_along_chromosome(ComputationType ct, ComputationResults cr, value<model> Mbase, double minpos, double maxpos)
Damien Leroux's avatar
Damien Leroux committed
518
519
520
        {
            Mcurrent.compute();
            if (output_model) {
521
522
                Mcurrent.output_X_to_file();
                Mcurrent.output_XtX_inv_to_file();
Damien Leroux's avatar
Damien Leroux committed
523
524
            }
            if (!loci->size()) {
Damien Leroux's avatar
Damien Leroux committed
525
                static computation_along_chromosome dummy_cac;  // = {{}, {}, {}, {}, {}, {}, {}, {}, {}};
Damien Leroux's avatar
Damien Leroux committed
526
527
                return dummy_cac;
            }
528
            testpos.clear();
529
            _recompute(testpos, minpos, maxpos);
530
531
            MSG_DEBUG("Rank(Mcurrent " << Mcurrent.keys() << ") = " << Mcurrent.rank());
            MSG_DEBUG("Rank(Mbase " << Mbase->keys() << ") = " << Mbase->rank());
Damien Leroux's avatar
Damien Leroux committed
532
533
            compute_along_chromosome(cac, ct, cr,
                                     Mcurrent, Mbase,
534
                                     /*Mcurrent, Mcurrent,*/
Damien Leroux's avatar
Damien Leroux committed
535
536
537
538
539
540
                                     locus_base_key,
                                     chrom_under_study,
                                     testpos,
                                     locus_blocks);

            if (output_test | output_rss | output_rank) {
541
542
                /*MSG_DEBUG(MATRIX_SIZE(cac.ftest_pvalue));*/
                /*MSG_DEBUG(MATRIX_SIZE(cac.rss));*/
Damien Leroux's avatar
Damien Leroux committed
543
544
                filename_stream filename;
                filename << Mcurrent.keys();
545
                if (output_test) { filename << computation_type_to_string(ct); }
Damien Leroux's avatar
Damien Leroux committed
546
547
                if (output_rss) { filename << '_' << "RSS"; }
                if (output_rank) { filename << '_' << "Rank"; }
548
                filename << (pmode == Joint ? "_Joint" : "");
Damien Leroux's avatar
Damien Leroux committed
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
                filename << ".txt";
                std::ofstream f(filename);
                if (output_test) { f << '\t' << "Test"; }
                if (output_rss) { for (int i = 0; i < cac.rss.innerSize(); ++i) { f << '\t' << "RSS"; } }
                if (output_rank) { f << '\t' << "Rank"; }
                f << std::endl;
                for (size_t i = 0; i < testpos.size(); ++i) {
                    f << testpos[i];
                    if (output_test) {
                        switch(test_type) {
                            case ComputationType::FTest:
                                f << '\t' << cac.ftest_pvalue(0, i);
                                break;
                            case ComputationType::FTestLOD:
                                f << '\t' << cac.ftest_lod(0, i);
                                break;
                            case ComputationType::Chi2:
                                f << '\t' << cac.chi2(0, i);
                                break;
                            default:
                                last_computation = NULL;
                        };
                    }
                    if (output_rss) {
                        for (int j = 0; j < cac.rss.innerSize(); ++j) {
                            f << '\t' << cac.rss(j, i);
                        }
                    }
                    if (output_rank) {
                        f << '\t' << cac.rank(i);
                    }
                    f << std::endl;
                }
            }
583

Damien Leroux's avatar
Damien Leroux committed
584
            return cac;
585
        }
Damien Leroux's avatar
Damien Leroux committed
586
587

    test_result
588
        test_along_chromosome(double min_pos, double max_pos)
Damien Leroux's avatar
Damien Leroux committed
589
590
591
592
        {
            if (!loci->size()) {
                last_computation = NULL;
                return test_result();
593
            }
Damien Leroux's avatar
Damien Leroux committed
594
595
596
            ComputationResults cr = Test;
            if (output_rss) { cr = cr | RSS; }
            if (output_rank) { cr = cr | Rank; }
597
598
599
600
601
602
603
604
605
606
607

            value<model> Mbase;
            if (pmode == Joint) {
                Mbase = Mcurrent.reduce(chrom_under_study);
                Mbase->compute();
            } else {
                Mbase = Mcurrent.clone();
            }

            custom_test_along_chromosome(test_type, cr, Mbase, min_pos, max_pos);

Damien Leroux's avatar
Damien Leroux committed
608
609
610
611
612
613
614
615
616
617
618
619
620
            switch(test_type) {
                case ComputationType::FTest:
                    last_computation = &cac.ftest_pvalue;
                    break;
                case ComputationType::FTestLOD:
                    last_computation = &cac.ftest_lod;
                    break;
                case ComputationType::Chi2:
                    last_computation = &cac.chi2;
                    break;
                default:
                    last_computation = NULL;
            };
Damien Leroux's avatar
WIP.    
Damien Leroux committed
621

Damien Leroux's avatar
Damien Leroux committed
622
            return find_max();
623
        }
Damien Leroux's avatar
Damien Leroux committed
624

625
626
627
628
629
630
631
#if 0
    double is_ghost_qtl(double locus)
    {
        locus_key selection = get_selection(chrom_under_study);
        double highest_test = 0;
        double highest_locus = -1;
        auto pmode_backup = pmode;
632
        set_product_mode();
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
        model m0 = Mcurrent;
        m0.add(chrom_under_study, locus);

        for (double l: selection) {
            model mA = m0;
            mA.add(chrom_under_study, l);
            model mB = mA;
            mB.ghost_constraint(chrom_under_study, locus, l);
            MSG_DEBUG("GHOST CHECK");
            MSG_DEBUG(mA);
            MSG_DEBUG(mB);
            /* FIXME: which test? */
            double test = ftest(mA, mB);
            if (test > highest_test) {
                highest_test = test;
            }
        }
        pmode = pmode_backup;
        return highest_test;
    }
#endif

Damien Leroux's avatar
Damien Leroux committed
655
656
657
658
659
660
    test_result
        find_max()
        {
            if (last_computation == NULL) {
                return {};
            }
661
662
663
            if (testpos.size() != (size_t)(last_computation->outerSize())) {
                MSG_ERROR("LOCI INCONSISTENT WITH COMPUTATION RESULT (" << loci->size() << " vs " << last_computation->outerSize() << ")", "");
            }
Damien Leroux's avatar
Damien Leroux committed
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
            int i_max = -1;
            double max = -1;
            for (int i = 0; i < last_computation->outerSize(); ++i) {
                if ((*last_computation)(0, i) >= max) {
                    max = (*last_computation)(0, i);
                    i_max = i;
                }
            }
            if (i_max == -1) {
                return {};
            }
            /*model_block_key k = locus_base_key;*/
            model_block_key mbk = locus_base_key;
            mbk += std::make_pair(chrom_under_study, (*loci)[i_max]);
            /*MSG_DEBUG("locus_base_key " << locus_base_key << " mbk " << mbk);*/
679
680
681
682
            /*MSG_DEBUG("last_computation@" << last_computation);*/
            /*MSG_QUEUE_FLUSH();*/
            /*MSG_DEBUG((*last_computation));*/
            /*MSG_QUEUE_FLUSH();*/
Damien Leroux's avatar
WIP.    
Damien Leroux committed
683

684
#ifdef SIGNAL_DISPLAY_ONELINER
685
            signal_display sd(last_computation->transpose(), i_max, max > threshold);
Damien Leroux's avatar
WIP.    
Damien Leroux committed
686
            MSG_DEBUG("[COMPUTATION] " << loci->front() << sd << loci->back() << " max=" << max << " at " << (*loci)[i_max]);
687
688
689
690
#else
            signal_display sd(*chrom_under_study, testpos, last_computation->transpose(), i_max, threshold);
            MSG_DEBUG("[COMPUTATION] " << loci->front() << " ... " << loci->back() << " max=" << max << " at " << (*loci)[i_max] << std::endl << sd);
#endif
Damien Leroux's avatar
WIP.    
Damien Leroux committed
691

Damien Leroux's avatar
Damien Leroux committed
692
            return {chrom_under_study,
693
                testpos[i_max], max, i_max, max > threshold,
694
695
                mbk,
                locus_blocks[i_max]};
Damien Leroux's avatar
Damien Leroux committed
696
        }
697
698
699
700

    test_result
        find_max_over_all_chromosomes()
        {
Damien Leroux's avatar
Damien Leroux committed
701
            test_result ret;
702
703
704
705
706
            for (const chromosome& c: active_settings->map) {
                select_chromosome(&c);
                if (!loci->size()) {
                    continue;
                }
707
                auto tmp = test_along_chromosome(0, max_testpos());
708
709
710
711
712
713
714
715
716
717
718
719
                if (tmp.test_value > ret.test_value) {
                    ret = tmp;
                }
            }
            return ret;
        }

/* Adding/Removing cofactors and QTLs
 */
    void add(const test_result& tr)
    {
        add(tr.block_key, tr.block);
Damien Leroux's avatar
Damien Leroux committed
720
721
722
723
        /*remove_test_locus(tr.chrom, tr.locus);*/
        if (pmode == Joint) {
            locus_base_key = tr.block_key;
        }
724
725
    }

Damien Leroux's avatar
Damien Leroux committed
726
    void remove_test_locus(const chromosome* chr, double loc)
727
728
729
730
731
732
733
734
735
736
737
738
    {
        std::vector<double>& llist = *all_loci[chr];
        auto it = std::find(llist.begin(), llist.end(), loc);
        if (it != llist.end()) {
            llist.erase(it);
            /*MSG_DEBUG("after remove_locus: " << llist);*/
        } else {
            MSG_ERROR("Locus " << loc << " doesn't exist in the list for chromosome " << chr->name, "");
            /*MSG_INFO("Available loci are " << llist);*/
        }
    }

Damien Leroux's avatar
Damien Leroux committed
739
740
741
742
743
744
745
    void add_test_locus(const chromosome* chr, double loc)
    {
        std::vector<double>& llist = *all_loci[chr];
        llist.push_back(loc);
        std::sort(llist.begin(), llist.end());
    }

746
747
748
749
750
751
752
753
754
755
756
757
758
    void add(const model_block_key& mbk, const value<model_block_type>& mb)
    {
        /*MSG_DEBUG("ADD " << mbk << " columns " << mb->column_labels);*/
        if (pmode == Joint) {
            /* replace previous block with the new one */
            Mcurrent.remove_blocks_with_chromosome(chrom_under_study);
        }
        /* add new block */
        Mcurrent.add_block(mbk, mb);
        Mcurrent.compute();
        /*MSG_DEBUG(Mcurrent.X());*/
    }

Damien Leroux's avatar
Damien Leroux committed
759
760
761
762
763
764
765
766
767
768
769
770
#if 1
    locus_key get_selection(const chromosome* chr)
    {
        if (pmode == Joint) {
            auto bwc = Mcurrent.blocks_with_chromosome(chr);
            if (bwc.begin() != bwc.end()) {
                return bwc.begin()->first[chr];
            }
        }
        return {};
    }

771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
    double
        get_threshold()
        {
            if (pmode == Joint) {
                model Mperm(trait_permutations(trait_name, active_settings->n_permutations), Mcurrent);
                Mperm.compute();
                std::string old_title = active_settings->set_title(MESSAGE("Computing threshold for trait " << trait_name << " given selection " << Mperm.keys() << " using " << active_settings->n_permutations << " and quantile " << active_settings->qtl_threshold_quantile));
                double ret = *make_value<Disk|Sync>(qtl_threshold_all_chromosomes_for_model,
                                                    as_value(trait_name),
                                                    as_value(active_settings->qtl_threshold_quantile),
                                                    as_value(active_settings->n_permutations),
                                                    as_value(Mperm));
                MSG_INFO("New threshold computed given " << Mcurrent.keys() << ": " << ret);
                active_settings->set_title(old_title);
                return ret;
            } else {
                return threshold;
            }
        }

Damien Leroux's avatar
Damien Leroux committed
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
    void set_selection(const model_block_collection& mbk) {
        /*MSG_DEBUG(__LINE__ << ": set_selection pmode=" << pmode << " selection=" << Mcurrent.keys());*/
        if (pmode == Joint) {
            /* replace previous block with the new one */
            Mcurrent.remove_blocks_with_chromosome(chrom_under_study);
            locus_base_key = mbk.begin()->first;
        }
        /*MSG_DEBUG(__LINE__ << ": set_selection pmode=" << pmode << " selection=" << Mcurrent.keys());*/
        Mcurrent.add_blocks(mbk);
        /*MSG_DEBUG(__LINE__ << ": set_selection pmode=" << pmode << " selection=" << Mcurrent.keys());*/
        Mcurrent.compute();
    }

    void set_selection(const locus_key& lk)
    {
        /*MSG_DEBUG("POUET " << pmode);*/
        /*MSG_DEBUG(__LINE__ << ": set_selection pmode=" << pmode << " selection=" << Mcurrent.keys());*/
        Mcurrent.remove_blocks_with_chromosome(chrom_under_study);
        if (pmode == Joint) {
            /* replace previous block with the new one */
            /*MSG_DEBUG(__LINE__ << ": set_selection pmode=" << pmode << " selection=" << Mcurrent.keys());*/
            if (!lk->is_empty()) {
                locus_base_key.selection.clear();
                locus_base_key.selection[chrom_under_study] = lk;
                /*std::vector<double> l = {lk->locus};*/
                std::vector<double> tmp_all = lk->to_vec();
                std::vector<double> tmp_test = {lk->locus};
                auto vmb = compute_parental_origins_multipop(all_pops,
                                                             as_value(chrom_under_study),
                                                             as_value(lk->parent),
                                                             as_value(tmp_all),
                                                             tmp_test)[0];
                /* add new block */
                /*MSG_DEBUG("ADD " << locus_base_key << " lk=" << lk << std::endl << vmb);*/
                Mcurrent.add_block(locus_base_key, vmb);
            }
        }
        /*MSG_DEBUG(__LINE__ << ": set_selection pmode=" << pmode << " selection=" << Mcurrent.keys());*/
        Mcurrent.compute();
    }

    void add(const chromosome* chr, double loc)
    {
        locus_key sel = get_selection(chr);
        model_block_key mbk;
        mbk.selection[chr] = sel + loc;
        std::vector<double> l = {loc};
        /*MSG_DEBUG(__LINE__ <<":chromosome " << chrom_under_study->name << " loci " << l); MSG_QUEUE_FLUSH();*/
        auto vmb = compute_parental_origins_multipop(all_pops,
                                                     as_value(chr),
                                                     as_value(sel),
                                                     as_value(l),
                                                     l)[0];
        /* add new block */
        if (pmode == Joint) {
            Mcurrent.remove_blocks_with_chromosome(chr);
        }
        Mcurrent.add_block(mbk, vmb);
    }

851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
    void
        add(const chromosome* chr, const locus_key& sel)
        {
            model_block_key mbk;
            std::vector<double> loci = sel->to_vec();
            std::vector<double> l = {loci.back()};
            loci.pop_back();
            mbk.selection[chr] = sel;
            locus_key presel = sel - l.front();
            MSG_DEBUG("Compute parental origins chr=" << chr->name << " sel " << sel << " presel " << presel << " locus " << l);
            auto vmb = compute_parental_origins_multipop(all_pops,
                                                         as_value(chr),
                                                         as_value(presel),
                                                         as_value(l),
                                                         l)[0];
            Mcurrent.add_block(mbk, vmb);
        }

    void
        add_sub_blocks(const chromosome* chr, const locus_key& lk)
        {
            if (!lk || lk->is_empty()) {
                return;
            }
            MSG_DEBUG("ADDING SUB-BLOCKS FOR " << chr->name << " " << lk);
            std::vector<double> vsel = lk->to_vec();
            if (vsel.size() > 1) {
                for (double l: vsel) {
                    locus_key k = lk - l;
                    MSG_DEBUG("sub key " << k);
                    add_sub_blocks(chr, k);
                    add(chr, k);
                }
            }
        }

    model_manager
        model_for_estimation()
        {
            model_manager ret(*this);
            for (const auto& kb: Mcurrent.m_blocks) {
                for (const auto& chr_sel: kb.first.selection) {
                    if (chr_sel.second->is_empty()) { continue; }
                    ret.add_sub_blocks(chr_sel.first, chr_sel.second);
                }
            }
            MSG_DEBUG("New keys: " << ret.Mcurrent.keys());
            ret.Mcurrent.compute();
            return ret;
        }

Damien Leroux's avatar
Damien Leroux committed
902
903
904
905
906
907
908
909
910
911
912
913
914
915
    void add_all_loci()
    {
        for (const auto& kv: all_loci) {
            if (!kv.second) {
                continue;
            }
            for (double loc: *kv.second) {
                add(kv.first, loc);
            }
        }
    }
#endif

    model_block_collection remove(double loc)
916
917
918
919
920
921
922
923
924
925
926
927
928
929
    {
        model_block_collection
            mbc = Mcurrent.extract_blocks_with_chromosome_and_locus(
                    chrom_under_study, loc);
        /* FIXME: there should be a special treatment for epistasis blocks on multiple chromosomes */
        for (auto& kv: mbc) {
            model_block_key k = kv.first - std::make_pair(chrom_under_study, loc);
            if (!k.empty()) {
                /* If this block still has any locus in it, perform reduction of the joint probabilities */
                locus_key lk = kv.first[chrom_under_study];
                locus_key lk2 = k[chrom_under_study];
                auto pop_it = active_settings->populations.begin();
                auto npar_it = n_par.begin();
                /* disassemble population blocks */
Damien Leroux's avatar
Damien Leroux committed
930
931
932
                auto pb = disassemble_parental_origins_multipop(chrom_under_study, lk->parent, *kv.second, all_pops);
                std::vector<collection<parental_origin_per_locus_type>> all_popl;
                all_popl.reserve(pb.size());
933
934
                for (auto& vmat: pb) {
                    /* reduce each population block */
935
                    const qtl_pop_type* pop = &pop_it->second;
936
937
                    context_key ck(new context_key_struc(pop, chrom_under_study, std::vector<double>()));
                    MatrixXd red = lk->reduce(*npar_it, loc);
Damien Leroux's avatar
Damien Leroux committed
938
939
940
941
942
943
944
945
946
                    /*MSG_DEBUG(MATRIX_SIZE(vmat->data));*/
                    /*MSG_DEBUG(MATRIX_SIZE(red));*/
                    /*MSG_DEBUG("vmat BEFORE red" << std::endl << vmat);*/
                    MatrixXd data = vmat->data * red;
                    vmat->data = data;
                    vmat->column_labels = get_stpom_data(ck, lk2->parent)->row_labels;
                    /*MSG_DEBUG("vmat AFTER red" << std::endl << vmat);*/
                    all_popl.emplace_back();
                    all_popl.back().push_back(vmat);
947
948
949
                }
                /* reassemble model block */
                auto new_block = assemble_parental_origins_multipop(as_value(chrom_under_study),
Damien Leroux's avatar
Damien Leroux committed
950
951
952
953
                                                                    as_value(lk2->parent),
                                                                    all_popl,
                                                                    all_pops)[0];
                /*MSG_DEBUG("Resulting block" << std::endl << new_block);*/
954
955
956
                Mcurrent.add_block(k, new_block);
            }
        }
Damien Leroux's avatar
Damien Leroux committed
957
958
959
960
961
962
963
964
965
966
        if (pmode == Joint) {
            locus_base_key -= std::make_pair(chrom_under_study, loc);
        }
        Mcurrent.compute();
        return mbc;
    }

    void add(const model_block_collection& mbc)
    {
        Mcurrent.add_blocks(mbc);
967
968
    }

Damien Leroux's avatar
Damien Leroux committed
969
970
    model::printable_keys keys() const { return Mcurrent.keys(); }

971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
    struct QTL {
        std::string chromosome;
        double locus;
        std::vector<double> LOD_loci;
        std::vector<double> LOD;

        QTL(const std::string& n, double l, const std::vector<double>& x, const MatrixXd& y)
            : chromosome(n), locus(l), LOD_loci(x), LOD(y.data(), y.data() + y.size())
        {
            /*MSG_DEBUG("QTL at " << chromosome << ':' << locus);*/
            /*MSG_DEBUG(y);*/
            /*MSG_DEBUG(MATRIX_SIZE(y));*/
            /*MSG_DEBUG("" << LOD);*/
        }

        std::pair<double, double>
            confidence_interval() const
            {
                /*MSG_DEBUG_INDENT_EXPR("[Confidence interval] ");*/
                /*MSG_DEBUG("LOD: " << LOD);*/
                double maxLOD = *std::max_element(LOD.begin(), LOD.end());
                double lod_cap = maxLOD - 1.5;
                /*MSG_DEBUG("max=" << maxLOD << " threshold=" << lod_cap);*/
                int i;
                for (i = 0; i < (int) LOD_loci.size() && LOD_loci[i] < locus && LOD[i] < lod_cap; ++i);
                /*MSG_DEBUG("LEFT i=" << i);*/
                double left = LOD_loci[i];
                for (i = LOD_loci.size() - 1; i >= 0 && LOD_loci[i] > locus && LOD[i] < lod_cap; --i);
                /*MSG_DEBUG("RIGHT i=" << i);*/
                double right = LOD_loci[i];
For faster browsing, not all history is shown. View entire blame