File: coverage.hpp

package info (click to toggle)
spades 3.13.1+dfsg-2
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 22,172 kB
  • sloc: cpp: 136,213; ansic: 48,218; python: 16,809; perl: 4,252; sh: 2,115; java: 890; makefile: 507; pascal: 348; xml: 303
file content (112 lines) | stat: -rw-r--r-- 3,324 bytes parent folder | download | duplicates (4)
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
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
//***************************************************************************
//* Copyright (c) 2015 Saint Petersburg State University
//* Copyright (c) 2011-2014 Saint Petersburg Academic University
//* All Rights Reserved
//* See file LICENSE for details.
//***************************************************************************

/*
 * coverage.hpp
 *
 *  Created on: Jun 21, 2011
 *      Author: sergey
 */

#pragma once

#include "utils/logger/logger.hpp"
#include <iostream>
#include <vector>
#include <algorithm>
#include "math/xmath.h"
#include "action_handlers.hpp"

namespace omnigraph {

//todo save/load raw k-mer coverage
template<class Graph>
class CoverageIndex : public GraphActionHandler<Graph> {
    typedef typename Graph::VertexId VertexId;
    typedef typename Graph::EdgeId EdgeId;

    Graph& g_;

 public:
    CoverageIndex(Graph &g)
            : GraphActionHandler<Graph>(g, "CoverageIndex"), g_(g) {
    }

    /**
     * In NON averaged units
     */
    void SetRawCoverage(EdgeId e, unsigned cov) {
        g_.data(e).set_raw_coverage(cov);
    }

    void IncRawCoverage(EdgeId e, unsigned count) {
        g_.data(e).inc_raw_coverage((int)count);
    }

    void SetAvgCoverage(EdgeId e, double cov) {
        g_.data(e).set_raw_coverage((int) math::round(cov * (double) this->g().length(e)));
    }

    /**
     * Returns average coverage of the edge
     */
    double coverage(EdgeId edge) const {
        return (double) RawCoverage(edge) / (double) this->g().length(edge);
    }

    unsigned RawCoverage(EdgeId edge) const {
        return g_.data(edge).raw_coverage();
    }

    void HandleDelete(EdgeId edge) override {
        SetRawCoverage(edge, 0);
    }

    void HandleMerge(const std::vector<EdgeId>& old_edges, EdgeId new_edge) override {
        unsigned coverage = 0;
        for (auto it = old_edges.begin(); it != old_edges.end(); ++it) {
            coverage += RawCoverage(*it);
        }
        SetRawCoverage(new_edge, coverage);
    }

    void HandleGlue(EdgeId new_edge, EdgeId edge1, EdgeId edge2) override {
        SetRawCoverage(new_edge, RawCoverage(edge1) + RawCoverage(edge2));
    }

    void HandleSplit(EdgeId old_edge, EdgeId new_edge1, EdgeId new_edge2) override {
        double avg_cov = coverage(old_edge);
        if (old_edge == g_.conjugate(old_edge)) {
            int raw1 = std::max(1, (int) math::round(avg_cov * (double) this->g().length(new_edge1)));
            SetRawCoverage(new_edge1, raw1);
            SetRawCoverage(g_.conjugate(new_edge1), raw1);
            SetRawCoverage(new_edge2, std::max(1, (int) math::round(avg_cov * (double) this->g().length(new_edge2))));
        } else {
            SetRawCoverage(new_edge1, std::max(1, (int) math::round(avg_cov * (double) this->g().length(new_edge1))));
            SetRawCoverage(new_edge2, std::max(1, (int) math::round(avg_cov * (double) this->g().length(new_edge2))));
        }
    }

    void Save(EdgeId e, std::ostream& out) const {
        out << fmt::format("{:.6f}", coverage(e));
    }

    void Load(EdgeId e, std::istream& in) {
        double cov;
        in >> cov;
        SetAvgCoverage(e, cov);
    }

    /*
     * Is thread safe if different threads process different edges.
     */
    bool IsThreadSafe() const override {
        return true;
    }
};

}