Skip to content
GitLab
Commit 82f7cf46 authored by Damien Leroux's avatar Damien Leroux
Browse files

Now adding dependant blocks to estimate epistatic effects.

parent d2270c93
Hide whitespace changes
Inline Side-by-side
3rd-party/ThreadPool/ThreadPool.h
View file @ 82f7cf46
......@@ -31,6 +31,12 @@ public:
title = s;
}
const std::string& get_title()
{
std::lock_guard<std::mutex> lock(title_mutex);
return title;
}
private:
// need to keep track of threads so we can join them
std::vector< std::thread > workers;
......
include/braille_plot.h
View file @ 82f7cf46
......@@ -6,6 +6,7 @@
struct palette {
struct rgb { int r, g, b; };
typedef int palette_type[256][3];
static const palette_type& _256()
......@@ -119,6 +120,49 @@ struct braille_grid {
braille_grid&
set_background(int r, int g, int b) { bgR = r; bgG = g; bgB = b; return *this; }
braille_grid
compose_vert(bool leftright, const braille_grid& other, bool other_leftright)
{
int w = 2 * std::max(m_data.cols(), other.m_data.cols());
int h = 4 * (m_data.rows() + other.m_data.rows());
int this_x = ((leftright ? w - m_width : 0) & ((int) ~1)) >> 1;
int this_y = 0;
int other_x = ((other_leftright ? w - other.m_width : 0) & ((int) ~1)) >> 1;
int other_y = ((m_height + 3) & ((int) ~3)) >> 2;
int r, c;
braille_grid ret(w, h);
ret.set_background((bgR + other.bgR) >> 1, (bgG + other.bgG) >> 1, (bgB + other.bgB) >> 1);
for (int y = 0; y < m_data.rows(); ++y) {
for (int x = 0; x < m_data.cols(); ++x) {
r = this_y + y;
c = this_x + x;
ret.m_colors[r][c] = m_colors[y][x];
}
}
for (int y = 0; y < other.m_data.rows(); ++y) {
for (int x = 0; x < other.m_data.cols(); ++x) {
r = other_y + y;
c = other_x + x;
ret.m_colors[r][c] = other.m_colors[y][x];
}
}
if (leftright) {
ret.m_pixel_count.topRightCorner(m_pixel_count.rows(), m_pixel_count.cols()) = m_pixel_count;
ret.m_data.topRightCorner(m_data.rows(), m_data.cols()) = m_data;
} else {
ret.m_pixel_count.topLeftCorner(m_pixel_count.rows(), m_pixel_count.cols()) = m_pixel_count;
ret.m_data.topLeftCorner(m_data.rows(), m_data.cols()) = m_data;
}
if (other_leftright) {
ret.m_pixel_count.bottomRightCorner(other.m_pixel_count.rows(), other.m_pixel_count.cols()) = other.m_pixel_count;
ret.m_data.bottomRightCorner(m_data.rows(), m_data.cols()) = m_data;
} else {
ret.m_pixel_count.bottomLeftCorner(other.m_pixel_count.rows(), other.m_pixel_count.cols()) = other.m_pixel_count;
ret.m_data.bottomLeftCorner(other.m_data.rows(), other.m_data.cols()) = other.m_data;
}
return ret;
}
const std::vector<int>
get_background() { return {bgR, bgG, bgB}; }
......@@ -330,6 +374,27 @@ struct braille_grid {
.line(x1, y0, x0, y0, dash_length, space_length, r, g, b);
}
braille_grid&
filled_triangle(int x0, int y0, int x1, int y1, int x2, int y2, int r=255, int g=255, int b=255)
{
std::vector<std::pair<int, int>> V = {{x0, y0}, {x1, y1}, {x2, y2}};
std::sort(V.begin(), V.end(), [] (const std::pair<int, int>& v1, const std::pair<int, int>& v2) { return v1.second < v2.second; });
std::tie(x0, y0) = V[0];
std::tie(x1, y1) = V[1];
std::tie(x2, y2) = V[2];
/* draw two triangles with a horizontal base using horizontal lines */
auto X = [] (int x0, int y0, int x1, int y1, int y) { return x0 + (x1 - x0) * (y - y0) / (y1 - y0); };
int y = y0;
for (; y <= y1; ++y) {
line(X(x0, y0, x1, y1, y), y, X(x0, y0, x2, y2, y), y, 1, 0, r, g, b);
}
for (; y <= y2; ++y) {
line(X(x1, y1, x2, y2, y), y, X(x0, y0, x2, y2, y), y, 1, 0, r, g, b);
}
return *this;
}
braille_grid&
filled_box(int x0, int y0, int x1, int y1, int r=255, int g=255, int b=255)
{
......
include/computations/frontends2.h
View file @ 82f7cf46
......@@ -54,6 +54,7 @@ ftest_along_chromosome(chromosome_value chr, const locus_key& lk,
struct signal_display {
#ifdef SIGNAL_DISPLAY_ONELINER
static const char* tick(double x)
{
static const char* ticks[9] = { " ", "\u2581", "\u2582", "\u2583", "\u2584", "\u2585", "\u2586", "\u2587", "\u2588" };
......@@ -109,6 +110,36 @@ struct signal_display {
}
return os << ']' << _NORMAL;
}
#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
};
......@@ -116,6 +147,8 @@ struct signal_display {
enum probability_mode { Joint, Single };
struct model_manager {
/* name of the studied trait */
std::string trait_name;
/* populations used in this model */
collection<population_value> all_pops;
/* POP matrices per locus per population (in the order of active_settings->populations) */
......@@ -154,13 +187,17 @@ struct model_manager {
/* slight hack. */
std::vector<double> testpos;
/* additional data for output and postprocessing */
std::map<std::string, std::map<double, std::pair<double, double>>> qtl_confidence_interval;
/* Construction
*/
model_manager(const collection<population_value>& colpops,
model_manager(const std::string& trait, const collection<population_value>& colpops,
const value<MatrixXd>& y, double th,
ComputationType ct = ComputationType::FTest,
SolverType st = SolverType::QR)
: all_pops(colpops)
: trait_name(trait)
, all_pops(colpops)
, pop_blocks()
, n_par()
, Mcurrent(y, colpops, st)
......@@ -178,9 +215,12 @@ struct model_manager {
, output_test(false)
, output_model(false)
, testpos()
, qtl_confidence_interval()
{
Mcurrent.add_block({}, cross_indicator(all_pops));
Mcurrent.compute();
MSG_DEBUG("INITIAL MODEL SIZE: " << MATRIX_SIZE(Mcurrent.X()));
MSG_QUEUE_FLUSH();
for (const auto& kpop: active_settings->populations) {
context_key ck(new context_key_struc(&kpop.second,
&active_settings->map[0],
......@@ -310,11 +350,25 @@ struct model_manager {
return chrom_under_study;
}
void
set_joint_mode()
{
MSG_DEBUG("[Model Manager] JOINT MODE");
pmode = Joint;
}
void
set_product_mode()
{
MSG_DEBUG("[Model Manager] PRODUCT MODE");
pmode = Single;
}
/* cofactors <-> QTLs for selected chromosome
*/
locus_key cofactors_to_QTLs()
{
pmode = Joint;
set_joint_mode();
model_block_collection mbc = Mcurrent.extract_blocks_with_chromosome(chrom_under_study);
model_block_key qtl_key = model_block_key::join(mbc);
if (qtl_key.empty()) {
......@@ -340,7 +394,7 @@ struct model_manager {
void QTLs_to_cofactors()
{
pmode = Single;
set_product_mode();
model_block_collection mbc = Mcurrent.extract_blocks_with_chromosome(chrom_under_study);
/* There shall be only ONE key in mbc. */
if (mbc.size() == 0) {
......@@ -473,8 +527,8 @@ struct model_manager {
}
testpos.clear();
_recompute(testpos, minpos, maxpos);
/*MSG_DEBUG("Rank(Mcurrent " << Mcurrent.keys() << ") = " << Mcurrent.rank());*/
/*MSG_DEBUG("Rank(Mbase " << Mbase->keys() << ") = " << Mbase->rank());*/
MSG_DEBUG("Rank(Mcurrent " << Mcurrent.keys() << ") = " << Mcurrent.rank());
MSG_DEBUG("Rank(Mbase " << Mbase->keys() << ") = " << Mbase->rank());
compute_along_chromosome(cac, ct, cr,
Mcurrent, Mbase,
/*Mcurrent, Mcurrent,*/
......@@ -575,7 +629,7 @@ struct model_manager {
double highest_test = 0;
double highest_locus = -1;
auto pmode_backup = pmode;
pmode = Single;
set_product_mode();
model m0 = Mcurrent;
m0.add(chrom_under_study, locus);
......@@ -627,8 +681,13 @@ struct model_manager {
/*MSG_DEBUG((*last_computation));*/
/*MSG_QUEUE_FLUSH();*/
#ifdef SIGNAL_DISPLAY_ONELINER
signal_display sd(last_computation->transpose(), i_max, max > threshold);
MSG_DEBUG("[COMPUTATION] " << loci->front() << sd << loci->back() << " max=" << max << " at " << (*loci)[i_max]);
#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
return {chrom_under_study,
testpos[i_max], max, i_max, max > threshold,
......@@ -709,6 +768,26 @@ struct model_manager {
return {};
}
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;
}
}
void set_selection(const model_block_collection& mbk) {
/*MSG_DEBUG(__LINE__ << ": set_selection pmode=" << pmode << " selection=" << Mcurrent.keys());*/
if (pmode == Joint) {
......@@ -769,6 +848,57 @@ struct model_manager {
Mcurrent.add_block(mbk, vmb);
}
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;
}
void add_all_loci()
{
for (const auto& kv: all_loci) {
......@@ -838,6 +968,70 @@ struct model_manager {
model::printable_keys keys() const { return Mcurrent.keys(); }
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];
/*MSG_INFO("Confidence interval for " << chromosome << ':' << locus << " {" << left << ':' << right << '}');*/
/*MSG_DEBUG_DEDENT;*/
return {left, right};
}
};
std::vector<QTL>
QTLs()
{
std::set<model_block_key> uniq_loci;
for (const auto& mb: Mcurrent.m_blocks) {
for (const auto& k: mb.first.split_loci()) {
uniq_loci.insert(k);
}
}
std::vector<QTL> ret;
for (const auto& ul: uniq_loci) {
/*chrom_under_study = ul.selection.begin()->first;*/
double loc = *std::begin(ul.selection.begin()->second);
select_chromosome(ul.selection.begin()->first);
model_manager sub(*this);
sub.output_model = sub.output_rank = sub.output_rss = sub.output_test = false;
sub.remove(loc);
/*testpos = *all_loci[chrom_under_study];*/
value<model> Mbase = sub.Mcurrent;
active_settings->set_title(MESSAGE("Computing confidence interval for QTL " << chrom_under_study->name << ':' << loc));
auto lod = sub.custom_test_along_chromosome(ComputationType::FTestLOD, ComputationResults::Test, Mbase, 0, sub.max_testpos()).ftest_lod;
/*MSG_DEBUG("ALL_LOCI " << testpos);*/
/*MSG_DEBUG("LOD " << lod);*/
ret.emplace_back(sub.chrom_under_study->name, loc, testpos, lod);
}
return ret;
}
private:
void _recompute(std::vector<double>& testpos, double minpos, double maxpos)
{
......@@ -851,6 +1045,9 @@ private:
void _recompute(const locus_key& lk, std::vector<double>& testpos, double minpos, double maxpos)
{
if (pmode == Joint) {
threshold = get_threshold();
}
/*testpos = compute_effective_loci(lk, *loci);*/
std::vector<double>::const_iterator mini, maxi;
std::vector<double>::const_iterator i = loci->begin(), j = loci->end();
......
include/computations/model.h
View file @ 82f7cf46
......@@ -171,7 +171,7 @@ init_connected_block_builder(
index.insert({l, sz});
}
}
#if 1
#if 0
MSG_DEBUG("index size " << index.size());
std::stringstream s;
s << '{';
......@@ -260,6 +260,8 @@ qtl_threshold(const std::string& trait_name, chromosome_value chr, const locus_k
double
qtl_threshold_all_chromosomes(const std::string& trait_name, double quantile, int n_permut);
double
qtl_threshold_all_chromosomes_for_model(const std::string& trait_name, double quantile, int n_permut, const model& Mperm);
std::vector<double>
compute_effective_loci(const locus_key& lk, const std::vector<double>& loci);
......
include/computations/to_md5.h
View file @ 82f7cf46
......@@ -41,6 +41,13 @@ md5_digest& operator << (md5_digest& md5, const context_key& ck)
return md5 << ck->pop->name << ck->pop->gen->name << ck->loci;
}
/*static inline*/
/*md5_digest& operator << (md5_digest& md5, const model& m)*/
/*{*/
/*return md5 << m.keys();*/
/*}*/
#if 0
static inline
md5_digest& operator << (md5_digest& md5, const impl::generation_rs::segment_computer_t& sc)
......
include/data/chromosome.h
View file @ 82f7cf46
<
......@@ -24,6 +24,45 @@ struct marker_name_and_position {
}
};
struct compute_xya {
template <typename StartFun, typename EndFun, typename Payload>
void
operator () (int padding_left, double x_factor,
const std::map<double, Payload>& poi,
StartFun&& get_start, EndFun&& get_end,
std::vector<int>& X, std::vector<int>& Y,
std::vector<int>& anchor, size_t& ymax)
{
std::vector<int> xend;
for (const auto& pl: poi) {
anchor.push_back(round(padding_left + pl.first * x_factor));
/*int x = anchor.back() - 2 * (1 + pl.second.size());*/
int x = get_start(anchor.back(), pl.second);
int y;
auto yi = std::find_if(xend.begin(), xend.end(), [&] (int xend) { return xend < x; });
if (yi == xend.end()) {
/*xend.push_back(anchor.back() + 2);*/
xend.push_back(get_end(anchor.back(), pl.second));
y = xend.size() + 1;
yi = xend.end();
} else {
/**yi = anchor.back() + 2;*/
*yi = get_end(anchor.back(), pl.second);
y = 2 + (yi - xend.begin());
}
for (auto xe = xend.begin(); xe < yi; ++xe) { *xe = get_end(anchor.back(), pl.second); }
X.push_back(x);
Y.push_back(y);
}
/*MSG_DEBUG("X: " << X);*/
/*MSG_DEBUG("Y: " << Y);*/
ymax = xend.size();
/*MSG_DEBUG("ymax: " << ymax);*/
}
};
struct chromosome {
std::string name;
marker_name_and_position raw;
......@@ -128,7 +167,14 @@ struct chromosome {
haplotype_iterator end() const { return {*this, true}; }
braille_grid
pretty_print(size_t width, const std::map<double, std::string>& POI) const
pretty_print(size_t width, const std::map<double, std::string>& POI, const std::map<double, std::pair<double, double>>& ROI) const
{
int padding_left;
return pretty_print(width, POI, ROI, padding_left);
}
braille_grid
pretty_print(size_t width, const std::map<double, std::string>& POI, const std::map<double, std::pair<double, double>>& ROI, int& padding_left, bool title=true) const
{
std::map<double, std::string> marker_labels;
double length = raw.marker_locus.back();
......@@ -141,55 +187,78 @@ struct chromosome {
}
marker_labels[pos] = label.str();
}
MSG_DEBUG("marker_labels begin");
for (const auto& ml: marker_labels) { MSG_DEBUG(ml.first << ": " << ml.second); }
MSG_DEBUG("marker_labels end");
size_t llsz = marker_labels.begin()->second.size();
if (POI.size()) {
llsz = std::max(llsz, POI.begin()->second.size());
}
int padding_left = 2 * (llsz + 1);
double x_factor = (2 * (width - 1) - padding_left) / length;
auto compute_xya
= [&] (const std::map<double, std::string>& poi, std::vector<int>& X, std::vector<int>& Y, std::vector<int>& anchor, size_t& ymax)
/*MSG_DEBUG("marker_labels begin");*/
/*for (const auto& ml: marker_labels) { MSG_DEBUG(ml.first << ": " << ml.second); }*/
/*MSG_DEBUG("marker_labels end");*/
/* compute padding_left */
padding_left = 0;
double x0 = 0, x1 = length;
auto compute_padding_left
= [&] (const std::map<double, std::string>& poi)
{
std::vector<int> xend;
for (const auto& pl: poi) {
anchor.push_back(round(padding_left + pl.first * x_factor));
int x = anchor.back() - 2 * (1 + pl.second.size());
int y;
auto yi = std::find_if(xend.begin(), xend.end(), [&] (int xend) { return xend < x; });
if (yi == xend.end()) {
xend.push_back(anchor.back() + 2);
y = xend.size() + 1;
yi = xend.end();
} else {
*yi = anchor.back() + 2;
y = 2 + (yi - xend.begin());
size_t L = pl.second.size() + 1;
double xL = pl.first;
int min_pad = (int) ceil((x1 - x0) * L / (x1 - xL) + (x0 - xL) * width / (x1 - xL));
if (padding_left < min_pad) {
padding_left = min_pad;
}
for (auto xe = xend.begin(); xe < yi; ++xe) { *xe = anchor.back() + 2; }
X.push_back(x);
Y.push_back(y);
}
ymax = xend.size();
};
compute_padding_left(marker_labels);
compute_padding_left(POI);
padding_left *= 2;
/*size_t llsz = marker_labels.begin()->second.size();*/
/*if (POI.size()) {*/
/*llsz = std::max(llsz, POI.begin()->second.size());*/
/*}*/
/*int padding_left = 2 * (llsz + 1);*/
double x_factor = ((2 * width - 1) - padding_left) / length;
std::vector<int> label_x, label_y, label_anchor;
std::vector<int> poi_x, poi_y, poi_anchor;
size_t label_ymax, poi_ymax;
compute_xya(marker_labels, label_x, label_y, label_anchor, label_ymax);
compute_xya(POI, poi_x, poi_y, poi_anchor, poi_ymax);
braille_grid grid(2 * width, 12 + 4 * (label_ymax + poi_ymax));
std::vector<int> roi_x, roi_y, roi_anchor;
size_t label_ymax, poi_ymax, roi_ymax;
compute_xya()(padding_left, x_factor,
marker_labels,
[] (int a, const std::string& s) { return a - 2 * (1 + s.size()); },
[] (int a, const std::string&) { return a + 2; },