Skip to content

Commit

Permalink
reHC*, continue refactoring of old code in order to allow multi-allel…
Browse files Browse the repository at this point in the history 
…ic genotypes and haplotypes.
  • Loading branch information
@yp
yp committed Jun 20, 2012
1 parent ec5138c commit c0994f9
Show file tree
Hide file tree
Showing 7 changed files with 146 additions and 76 deletions.
13 changes: 5 additions & 8 deletions include/haplotypes_genotypes.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,9 +69,6 @@ class single_multiallelic_genotype_t {
public:

static const single_multiallelic_genotype_t MISS;
static const single_multiallelic_genotype_t HOMO1;
static const single_multiallelic_genotype_t HOMO2;
static const single_multiallelic_genotype_t HETER12;

single_multiallelic_genotype_t()
:_allele1(0), _allele2(0)
Expand Down Expand Up @@ -159,11 +156,6 @@ class single_multiallelic_haplotype_t {
public:

static const single_multiallelic_haplotype_t MISS;
static const single_multiallelic_haplotype_t ALLELE1;
static const single_multiallelic_haplotype_t ALLELE2;
static const single_multiallelic_haplotype_t ALLELE3;
static const single_multiallelic_haplotype_t ALLELE4;
static const single_multiallelic_haplotype_t ALLELE5;

single_multiallelic_haplotype_t()
:_allele(0)
Expand All @@ -181,6 +173,11 @@ class single_multiallelic_haplotype_t {
return _allele;
};

static single_multiallelic_haplotype_t ALLELE(const allele_t allele) {
MY_ASSERT_DBG(allele != single_multiallelic_haplotype_t::MISS.allele());
return single_multiallelic_haplotype_t(allele);
};

};

bool
Expand Down
10 changes: 7 additions & 3 deletions include/ped2cnf.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,17 +132,21 @@ class ped2cnf_conv_t:
lit_t e= cnf.get_e(i, l);
lit_t p= cnf.get_p(i, l);
lit_t m= cnf.get_m(i, l);
if (gen == g::HOMO1) {
// Alleles are sorted ! (i.e. gen.allele1() <= gen.allele2())
if (gen.allele1() == 1 && gen.allele2() == 1) {
// gen == HOMO1
// -e p m, e -p, e -m
cnf.add_clause<3>((lit_t[]){-e, p, m});
cnf.add_clause<2>((lit_t[]){ e, -p});
cnf.add_clause<2>((lit_t[]){ e, -m});
} else if (gen == g::HOMO2) {
} else if (gen.allele1() == 2 && gen.allele2() == 2) {
// gen == HOMO2
// -e -p -m, e p, e m
cnf.add_clause<3>((lit_t[]){-e, -p, -m});
cnf.add_clause<2>((lit_t[]){ e, p});
cnf.add_clause<2>((lit_t[]){ e, m});
} else if (gen == g::HETER12) {
} else if (gen.allele1() == 1 && gen.allele2() == 2) {
// gen == HETERO12
#ifdef AVOID_XOR_CLAUSES
cnf.add_clause<3>((lit_t[]){ e, p, m});
cnf.add_clause<3>((lit_t[]){ e, -p, -m});
Expand Down
61 changes: 53 additions & 8 deletions include/pedcnf.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,11 +56,11 @@


class ped_var_kind
:public enum_like_t<ped_var_kind, 8, 8>
:public enum_like_t<ped_var_kind, 10, 10>
{
private:

typedef enum_like_t<ped_var_kind, 8, 8> base;
typedef enum_like_t<ped_var_kind, 10, 10> base;

ped_var_kind(const int val)
:base(val)
Expand All @@ -79,6 +79,8 @@ class ped_var_kind
static const ped_var_kind RP; // Recombination events (from father)
static const ped_var_kind RM; // Recombination events (from mother)
static const ped_var_kind E; // Errors
static const ped_var_kind PM; // Paternal multi allele
static const ped_var_kind MM; // Maternal multi allele
static const ped_var_kind DUMMY; // Dummy

static const int int_values[];
Expand All @@ -96,11 +98,11 @@ class pedcnf_t
// Types
private:

typedef boost::tuple<size_t, size_t> index_var_t;
typedef boost::tuple<size_t, size_t, size_t> index_var_t;

public:

typedef boost::tuple<ped_var_kind, size_t, size_t> pedvar_t;
typedef boost::tuple<ped_var_kind, size_t, size_t, size_t> pedvar_t;
typedef std::map<index_var_t, lit_t> varmap_t;
typedef std::vector<pedvar_t> varvec_t;
typedef std::vector<bool> valvec_t;
Expand All @@ -126,6 +128,8 @@ class pedcnf_t
varmap_t _rp; // Recombination events (from father)
varmap_t _rm; // Recombination events (from mother)
varmap_t _e; // Errors
varmap_t _pm; // Paternal allele
varmap_t _mm; // Maternal allele
size_t _next_dummy;

varvec_t _vars;
Expand Down Expand Up @@ -157,18 +161,39 @@ class pedcnf_t
// Methods
private:

lit_t get_var3(varmap_t& map,
const ped_var_kind& var_kind,
const size_t i1, const size_t i2, const size_t i3);

lit_t get_var3(const varmap_t& map,
const size_t i1, const size_t i2, const size_t i3) const;

bool has_var3(const varmap_t& map,
const size_t i1, const size_t i2, const size_t i3) const;

bool get_val3(const varmap_t& map,
const size_t i1, const size_t i2, const size_t i3) const;

lit_t get_var(varmap_t& map,
const ped_var_kind& var_kind,
const size_t i1, const size_t i2);
const size_t i1, const size_t i2) {
return get_var3(map, var_kind, i1, i2, 0);
};

lit_t get_var(const varmap_t& map,
const size_t i1, const size_t i2) const;
const size_t i1, const size_t i2) const {
return get_var3(map, i1, i2, 0);
};

bool has_var(const varmap_t& map,
const size_t i1, const size_t i2) const;
const size_t i1, const size_t i2) const {
return has_var3(map, i1, i2, 0);
};

bool get_val(const varmap_t& map,
const size_t i1, const size_t i2) const;
const size_t i1, const size_t i2) const {
return get_val3(map, i1, i2, 0);
};


public:
Expand All @@ -195,6 +220,10 @@ class pedcnf_t

lit_t get_rm(const size_t i, const size_t l);

lit_t get_pm(const size_t i, const size_t l, const size_t j);

lit_t get_mm(const size_t i, const size_t l, const size_t j);

lit_t get_e(const size_t i, const size_t l);

bool has_sp(const size_t i, const size_t l) const;
Expand All @@ -209,6 +238,10 @@ class pedcnf_t

bool has_rm(const size_t i, const size_t l) const;

bool has_pm(const size_t i, const size_t l, const size_t j) const;

bool has_mm(const size_t i, const size_t l, const size_t j) const;

bool has_e(const size_t i, const size_t l) const;

lit_t generate_dummy();
Expand All @@ -225,6 +258,10 @@ class pedcnf_t

lit_t get_rm(const size_t i, const size_t l) const;

lit_t get_pm(const size_t i, const size_t l, const size_t j) const;

lit_t get_mm(const size_t i, const size_t l, const size_t j) const;

lit_t get_e(const size_t i, const size_t l) const;

bool sp(const size_t i, const size_t l) const {
Expand All @@ -251,6 +288,14 @@ class pedcnf_t
return get_val(_rm, i, l);
};

bool pm(const size_t i, const size_t l, const size_t j) const {
return get_val3(_pm, i, l, j);
};

bool mm(const size_t i, const size_t l, const size_t j) const {
return get_val3(_mm, i, l, j);
};

bool e(const size_t i, const size_t l) const {
return get_val(_e, i, l);
};
Expand Down
16 changes: 8 additions & 8 deletions include/pedcnf2hc.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,14 +78,14 @@ void compute_reHC_from_SAT(basic_pedigree_t<T_GENOTYPE, T_HAPLOTYPE, T_PHENOTYPE
DEBUG("Not-genotyped individual " << ind.progr_id() <<
" at locus " << locus << " is imputed as homozygous.");
if ( !pil ) {
ind.real_g(locus)= family_t::g::HOMO1;
ind.real_g(locus).set_alleles(1, 1);
} else {
ind.real_g(locus)= family_t::g::HOMO2;
ind.real_g(locus).set_alleles(2, 2);
}
} else {
DEBUG("Not-genotyped individual " << ind.progr_id() <<
" at locus " << locus << " is imputed as heterozygous.");
ind.real_g(locus)= family_t::g::HETER12;
ind.real_g(locus).set_alleles(1, 2);
}
++no_of_imputation;
} else {
Expand All @@ -95,17 +95,17 @@ void compute_reHC_from_SAT(basic_pedigree_t<T_GENOTYPE, T_HAPLOTYPE, T_PHENOTYPE
// Impossible: we performed genotype imputation in the previous step
MY_FAIL;
} else {
ind.hp(locus)= (!pil) ? (family_t::h::ALLELE1) : (family_t::h::ALLELE2);
ind.hm(locus)= (!mil) ? (family_t::h::ALLELE1) : (family_t::h::ALLELE2);
ind.hp(locus)= (!pil) ? (family_t::h::ALLELE(1)) : (family_t::h::ALLELE(2));
ind.hm(locus)= (!mil) ? (family_t::h::ALLELE(1)) : (family_t::h::ALLELE(2));
typename family_t::g real_g= ind.real_g(locus);
if ( pil == mil ) {
if ( !pil ) {
real_g= family_t::g::HOMO1;
real_g.set_alleles(1, 1);
} else {
real_g= family_t::g::HOMO2;
real_g.set_alleles(2, 2);
}
} else {
real_g= family_t::g::HETER12;
real_g.set_alleles(1, 2);
}
if (real_g != ind.real_g(locus)) {
// Not-genotyped loci cannot have errors
Expand Down
23 changes: 20 additions & 3 deletions src/lib/assumptions.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,11 +38,12 @@
std::istream&
operator>>(std::istream& in, pedcnf_t::pedvar_t& var) {
ped_var_kind var_kind= ped_var_kind::DUMMY;
size_t i1, i2;
size_t i1, i2, i3;
in >> var_kind;
in >> i1;
in >> i2;
var= boost::make_tuple(var_kind, i1, i2);
in >> i3;
var= boost::make_tuple(var_kind, i1, i2, i3);
return in;
};

Expand Down Expand Up @@ -88,7 +89,7 @@ class assumption_manager_t
// Process line
try {
std::istringstream is(buff);
pedcnf_t::pedvar_t v= boost::make_tuple(ped_var_kind::DUMMY, 0, 0);
pedcnf_t::pedvar_t v= boost::make_tuple(ped_var_kind::DUMMY, 0, 0, 0);
if (!(is >> v)) throw invalid_line_t("unrecognized variable");
bool value;
if (!(is >> value)) throw invalid_line_t("unrecognized boolean value");
Expand Down Expand Up @@ -144,6 +145,22 @@ class assumption_manager_t
lit_t iv= cnf.get_rm(v.get<1>(), v.get<2>());
cnf.add_clause<1>( (lit_t[]){ (value? +1 : -1) * iv } );
++n_assumptions;
} else if (v.get<0>() == ped_var_kind::PM) {
if (!cnf.has_pm(v.get<1>(), v.get<2>(), v.get<3>())) {
L_WARN("Variable '" << v << "' does not exist in the SAT instance."
" Adding anyway...");
}
lit_t iv= cnf.get_pm(v.get<1>(), v.get<2>(), v.get<3>());
cnf.add_clause<1>( (lit_t[]){ (value? +1 : -1) * iv } );
++n_assumptions;
} else if (v.get<0>() == ped_var_kind::MM) {
if (!cnf.has_mm(v.get<1>(), v.get<2>(), v.get<3>())) {
L_WARN("Variable '" << v << "' does not exist in the SAT instance."
" Adding anyway...");
}
lit_t iv= cnf.get_mm(v.get<1>(), v.get<2>(), v.get<3>());
cnf.add_clause<1>( (lit_t[]){ (value? +1 : -1) * iv } );
++n_assumptions;
} else if (v.get<0>() == ped_var_kind::E) {
if (!cnf.has_e(v.get<1>(), v.get<2>())) {
L_WARN("Variable '" << v << "' does not exist in the SAT instance."
Expand Down
25 changes: 0 additions & 25 deletions src/lib/haplotypes_genotypes.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,16 +40,6 @@
const single_multiallelic_genotype_t
single_multiallelic_genotype_t::MISS(0, 0);

const single_multiallelic_genotype_t
single_multiallelic_genotype_t::HOMO1(1, 1);

const single_multiallelic_genotype_t
single_multiallelic_genotype_t::HOMO2(2, 2);

const single_multiallelic_genotype_t
single_multiallelic_genotype_t::HETER12(1, 2);


bool operator==(const single_multiallelic_genotype_t& g1,
const single_multiallelic_genotype_t& g2) {
return (g1.allele1()==g2.allele1()) &&
Expand Down Expand Up @@ -110,21 +100,6 @@ is_heterozygous(const single_multiallelic_genotype_t& g) {
const single_multiallelic_haplotype_t
single_multiallelic_haplotype_t::MISS(0);

const single_multiallelic_haplotype_t
single_multiallelic_haplotype_t::ALLELE1(1);

const single_multiallelic_haplotype_t
single_multiallelic_haplotype_t::ALLELE2(2);

const single_multiallelic_haplotype_t
single_multiallelic_haplotype_t::ALLELE3(3);

const single_multiallelic_haplotype_t
single_multiallelic_haplotype_t::ALLELE4(4);

const single_multiallelic_haplotype_t
single_multiallelic_haplotype_t::ALLELE5(5);


bool operator==(const single_multiallelic_haplotype_t& h1,
const single_multiallelic_haplotype_t& h2) {
Expand Down
Loading

0 comments on commit c0994f9

Please sign in to comment.