pairwiseld.cc

#include <Sequence/summstats/ld.hpp>
#include <Sequence/VariantMatrix.hpp>

#include <vector>
namespace Sequence
{
    struct PairwiseLDStats
    {
        double i, j, rsq, D, Dprime;
        bool skipped;
    };

    PairwiseLDStats
    pairwise_ld_details(const VariantMatrix& m, const std::size_t i,
                        const std::size_t j, const std::int8_t refstate)
    {
        auto counts = two_locus_haplotype_counts(m, i, j, false);

        //We must have two alleles at each state and the reference
        //allele at each state
        int ref_i = 0, ref_j = 0, nonref_i = 0, nonref_j = 0, num_nonref_i = 0,
            num_nonref_j = 0;
        int num_nonref_nonref = 0;
        for (auto& c : counts)
            {
                if (c.i == refstate)
                    {
                        ref_i += c.n;
                    }
                else if (c.i >= 0)
                    {
                        nonref_i += c.n;
                        ++num_nonref_i;
                    }
                if (c.j == refstate)
                    {
                        ref_j += c.n;
                    }
                else if (c.j >= 0)
                    {
                        nonref_j += c.n;
                        ++num_nonref_j;
                    }
            }
    }

    std::vector<PairwiseLDStats>
    pairwise_ld(const Sequence::VariantMatrix& m, const std::int8_t refstate)
    {
        std::vector<PairwiseLDStats> rv;
        if (m.data.empty() || m.positions.empty() || !m.nsam || !m.nsites
            || m.nsites < 2)
            {
                return rv;
            }
        for (std::size_t i = 0; i < m.nsites - 1; ++i)
            {
                for (std::size_t j = i + 1; j < m.nsites; ++j)
                    {
                    }
            }
        return rv;
    }
} // namespace Sequence