~ubuntu-branches/ubuntu/raring/bedtools/raring

« back to all changes in this revision

Viewing changes to src/jaccard/jaccardMain.cpp

  • Committer: Package Import Robot
  • Author(s): Charles Plessy
  • Date: 2012-11-04 17:59:41 UTC
  • mfrom: (1.1.9)
  • Revision ID: package-import@ubuntu.com-20121104175941-hahqzy1w8uy650z0
Tags: 2.17.0-1
bb9012e Imported Upstream version 2.16.2.
9006b23 Imported Upstream version 2.17.0.
9112569 Documented BEDTools license as a whole.
325689c Removed Pre-Depends: dpkg (>= 1.15.6).
a781b14 Conforms with Policy 3.9.4.
84b1167 Use Debhelper 9. 
0bf572d Distribute the test suite.
422cd34 Bash completion for BEDTools.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/*****************************************************************************
 
2
  jaccardMain.cpp
 
3
 
 
4
  (c) 2009 - Aaron Quinlan
 
5
  Hall Laboratory
 
6
  Department of Biochemistry and Molecular Genetics
 
7
  University of Virginia
 
8
  aaronquinlan@gmail.com
 
9
 
 
10
  Licenced under the GNU General Public License 2.0 license.
 
11
******************************************************************************/
 
12
#include "jaccard.h"
 
13
#include "version.h"
 
14
 
 
15
using namespace std;
 
16
 
 
17
// define our program name
 
18
#define PROGRAM_NAME "bedtools jaccard"
 
19
 
 
20
 
 
21
// define our parameter checking macro
 
22
#define PARAMETER_CHECK(param, paramLen, actualLen) (strncmp(argv[i], param, min(actualLen, paramLen))== 0) && (actualLen == paramLen)
 
23
 
 
24
// function declarations
 
25
void jaccard_help(void);
 
26
 
 
27
int jaccard_main(int argc, char* argv[]) {
 
28
 
 
29
    // our configuration variables
 
30
    bool showHelp = false;
 
31
 
 
32
    // input files
 
33
    string bedAFile;
 
34
    string bedBFile;
 
35
 
 
36
    // input arguments
 
37
    float overlapFraction = 1E-9;
 
38
 
 
39
    bool haveBedA           = false;
 
40
    bool haveBedB           = false;
 
41
    bool haveFraction       = false;
 
42
    bool reciprocalFraction = false;
 
43
 
 
44
    // check to see if we should print out some help
 
45
    if(argc <= 1) showHelp = true;
 
46
 
 
47
    for(int i = 1; i < argc; i++) {
 
48
        int parameterLength = (int)strlen(argv[i]);
 
49
 
 
50
        if((PARAMETER_CHECK("-h", 2, parameterLength)) ||
 
51
        (PARAMETER_CHECK("--help", 5, parameterLength))) {
 
52
            showHelp = true;
 
53
        }
 
54
    }
 
55
 
 
56
    if(showHelp) jaccard_help();
 
57
 
 
58
    // do some parsing (all of these parameters require 2 strings)
 
59
    for(int i = 1; i < argc; i++) {
 
60
 
 
61
        int parameterLength = (int)strlen(argv[i]);
 
62
 
 
63
        if(PARAMETER_CHECK("-a", 2, parameterLength)) {
 
64
            if ((i+1) < argc) {
 
65
                haveBedA = true;
 
66
                bedAFile = argv[i + 1];
 
67
                i++;
 
68
            }
 
69
        }
 
70
        else if(PARAMETER_CHECK("-b", 2, parameterLength)) {
 
71
            if ((i+1) < argc) {
 
72
                haveBedB = true;
 
73
                bedBFile = argv[i + 1];
 
74
                i++;
 
75
            }
 
76
        }
 
77
        else if(PARAMETER_CHECK("-f", 2, parameterLength)) {
 
78
            if ((i+1) < argc) {
 
79
                haveFraction = true;
 
80
                overlapFraction = atof(argv[i + 1]);
 
81
                i++;
 
82
            }
 
83
        }
 
84
        else if(PARAMETER_CHECK("-r", 2, parameterLength)) {
 
85
            reciprocalFraction = true;
 
86
        }
 
87
        else {
 
88
            cerr << endl 
 
89
                 << "*****ERROR: Unrecognized parameter: " 
 
90
                 << argv[i] 
 
91
                 << " *****" 
 
92
                 << endl << endl;
 
93
            showHelp = true;
 
94
        }
 
95
    }
 
96
 
 
97
    // make sure we have both input files
 
98
    if (!haveBedA || !haveBedB) {
 
99
        cerr << endl << "*****" << endl << "*****ERROR: Need -a and -b files. " << endl << "*****" << endl;
 
100
        showHelp = true;
 
101
    }
 
102
 
 
103
    if (reciprocalFraction && !haveFraction) {
 
104
        cerr << endl << "*****" << endl << "*****ERROR: If using -r, you need to define -f." << endl << "*****" << endl;
 
105
        showHelp = true;
 
106
    }
 
107
 
 
108
    if (!showHelp) {
 
109
 
 
110
        Jaccard *j = new Jaccard(bedAFile, bedBFile,
 
111
                                 overlapFraction, reciprocalFraction);
 
112
        delete j;
 
113
        return 0;
 
114
    }
 
115
    else {
 
116
        jaccard_help();
 
117
        return 0;
 
118
    }
 
119
}
 
120
 
 
121
void jaccard_help(void) {
 
122
 
 
123
    cerr << "\nTool:    bedtools jaccard (aka jaccard)" << endl;
 
124
    cerr << "Version: " << VERSION << "\n";    
 
125
    cerr << "Summary: Calculate Jaccard statistic b/w two feature files."
 
126
         << endl
 
127
         << "\t Jaccard is the length of the intersection over the union."
 
128
         << endl
 
129
         << "\t Values range from 0 (no intersection) to 1 (self intersection)."
 
130
         << endl << endl;
 
131
 
 
132
    cerr << "Usage:   " << PROGRAM_NAME << " [OPTIONS] -a <bed/gff/vcf> -b <bed/gff/vcf>" << endl << endl;
 
133
 
 
134
    cerr << "Options: " << endl;
 
135
 
 
136
 
 
137
    cerr << "\t-f\t"            << "Minimum overlap required as a fraction of A." << endl;
 
138
    cerr                        << "\t\t- Default is 1E-9 (i.e., 1bp)." << endl;
 
139
    cerr                        << "\t\t- FLOAT (e.g. 0.50)" << endl << endl;
 
140
 
 
141
    cerr << "\t-r\t"            << "Require that the fraction overlap be reciprocal for A and B." << endl;
 
142
    cerr                        << "\t\t- In other words, if -f is 0.90 and -r is used, this requires" << endl;
 
143
    cerr                        << "\t\t  that B overlap 90% of A and A _also_ overlaps 90% of B." << endl << endl;
 
144
 
 
145
    cerr << "Notes: " << endl;
 
146
    cerr << "\t(1) Input files must be sorted by chrom, then start position." 
 
147
         << endl << endl;
 
148
 
 
149
    // end the program here
 
150
    exit(1);
 
151
 
 
152
}