signif.cpp 2.77 KB
Newer Older
RomainFeron's avatar
RomainFeron committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
/*
* Copyright (C) 2020 Romain Feron
* This file is part of RADSex.

* RADSex is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.

* RADSex is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.

* You should have received a copy of the GNU General Public License
* along with RADSex.  If not, see <https://www.gnu.org/licenses/>.
*/

19
20
21
22
/*!
 * @file signif.cpp
 * @brief Implements the Signif class.
*/
23

24
#include "signif.h"
25

26

27
void Signif::process_marker(Marker& marker) {
28

29
    if (marker.n_individuals > 0) {
30

31
        ++this->results.n_markers;  // Increment total number of markers for Bonferroni correction
32

33
        marker.p = get_p_association(marker.group_counts[parameters.group1], marker.group_counts[parameters.group2], this->popmap.get_count(this->parameters.group1), this->popmap.get_count(this->parameters.group2));
34

35
36
        // First pass: filter markers with non-corrected p < 0.05
        if (static_cast<float>(marker.p) < this->parameters.signif_threshold) this->results.candidate_markers.push_back(marker);
37
38
39
40
41
42

    }

}


43
44


45

46
void Signif::generate_output() {
47

48
    std::ofstream output_file = open_output(this->parameters.output_file_path);
49

50
    if (not this->parameters.output_fasta) {
51

52
        output_file << "#source:radsex-signif;min_depth:" << parameters.min_depth << ";signif_threshold:" << parameters.signif_threshold <<
RomainFeron's avatar
RomainFeron committed
53
54
                       ";bonferroni:" <<  std::boolalpha << (not parameters.disable_correction) <<
                       ";n_markers:" << this->results.n_markers << "\n";
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
        output_file << print_list(this->markers_table.header.header, "\t") << "\n";

    }

    if (not this->parameters.disable_correction) {

        // Bonferroni correction: divide threshold by total number of tests (i.e. number of retained markers)
        this->parameters.signif_threshold /= this->results.n_markers;

    } else {

        // Set number of markers to 1 so corrected p-values are the same as original p-values
        this->results.n_markers = 1;

    }
70

71
72
    // Second pass: filter markers with threshold corrected p < 0.05
    for (auto marker: this->results.candidate_markers) {
73

74
        if (static_cast<float>(marker.p) < this->parameters.signif_threshold) {
75

76
77
            marker.p_corr = std::min(1.0, marker.p * static_cast<double>(results.n_markers));

78
            this->parameters.output_fasta ? marker.output_as_fasta(output_file, this->parameters.min_depth) : marker.output_as_table(output_file);
79

80
81
82
        }

    }
83
84
85

    output_file.close();

86
}
87