-
Notifications
You must be signed in to change notification settings - Fork 4
Expand file tree
/
Copy pathdotplot.cpp
More file actions
93 lines (80 loc) · 3.1 KB
/
dotplot.cpp
File metadata and controls
93 lines (80 loc) · 3.1 KB
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
#include <Rcpp.h>
using namespace Rcpp;
//' Generate a dotplot matrix by comparing sequence 1 and sequence 2 using the
//' given window size and percent match
//'
//' @param seq1 The first sequence
//' @param seq2 The second sequence
//' @param wsize The window of the size to use
//' @param wstep the size in chars for the steps of the moving window. Use wstep == wsize for non-overlapping windows.
//' @param nmatch if the number of match per window is greater than or equal to nmatch then a dot is produced.
//' @return Returns a matrix
//' @export
// [[Rcpp::export]]
NumericMatrix mkDotPlotMatrix(std::string seq1, std::string seq2, int wsize, int wstep, int nmatch ) {
// Check arguments:
if(wsize < 1) stop("non allowed value for wsize");
if(wstep < 1) stop("non allowed value for wstep");
if(nmatch < 1) stop("non allowed value for nmatch");
if(nmatch > wsize) stop("nmatch > wsize is not allowed");
int s1l = seq1.length();
int s2l = seq2.length();
int xl = 1 + (s1l - wsize) / wstep;
int yl = 1 + (s2l - wsize) / wstep;
NumericMatrix X(xl, yl);
std::fill(X.begin(), X.end(), 1.0);
// Brute force, inefficient with overlapping comparisons
for(int i=0; i < (s1l - wsize); i += wstep) {
if( (i / wstep) == 1000) { Rcpp::checkUserInterrupt(); }
for(int j=0; j < (s2l - wsize); j += wstep) {
if( (j / wstep) == 1000) { Rcpp::checkUserInterrupt(); }
int nm =0;
for(int k=0; k < wsize; k++) {
if (seq1[i+k] == seq2[j+k]) nm++;
}
if (nm >= nmatch) {
X(i, j) = 0.0;
}
}
}
return X;
}
//' Generate a dotplot matrix by comparing sequence 1 and sequence 2 using the
//' given window size and percent match
//'
//' @param seq1 The first sequence
//' @param seq2 The second sequence
//' @param wsize The window of the size to use
//' @param wstep the size in chars for the steps of the moving window. Use wstep == wsize for non-overlapping windows.
//' @param nmatch if the number of match per window is greater than or equal to nmatch then a dot is produced.
//' @return Returns a matrix
//' @export
// [[Rcpp::export]]
DataFrame mkDotPlotDataFrame(std::string seq1, std::string seq2, int wsize, int wstep, int nmatch ) {
// Check arguments:
if(wsize < 1) stop("non allowed value for wsize");
if(wstep < 1) stop("non allowed value for wstep");
if(nmatch < 1) stop("non allowed value for nmatch");
if(nmatch > wsize) stop("nmatch > wsize is not allowed");
int s1l = seq1.length();
int s2l = seq2.length();
std::vector<double> X;
std::vector<double> Y;
// Brute force, inefficient with overlapping comparisons
for(int i=0; i < (s1l - wsize); i += wstep) {
if( (i / wstep) == 1000) { Rcpp::checkUserInterrupt(); }
for(int j=0; j < (s2l - wsize); j += wstep) {
if( (j / wstep) == 1000) { Rcpp::checkUserInterrupt(); }
int nm =0;
for(int k=0; k < wsize; k++) {
if (seq1[i+k] == seq2[j+k]) { nm++; }
}
if (nm >= nmatch) {
X.push_back(static_cast<double>(i));
Y.push_back(static_cast<double>(j));
}
}
}
return DataFrame::create(_["x"] = X,
_["y"] = Y);
}