dparasail source code

1 module dparasail;
2 
3 public import dparasail.alignment;
4 public import dparasail.result;
5 
6 
7 unittest
8 {
9     import std.stdio;
10     import std.math : approxEqual;
11 
12     auto p = Parasail("dnafull", 3, 2);
13     auto q = p.sw_striped("GATTA", "GACTA");
14 
15     assert(q.queryPosition == 0);
16     assert(q.position == 0);
17     assert(q.score == 16);
18     assert(q.cigar.alignedLength == 5);
19     assert(q.cigar.toString == "2=1X2=");
20     q.writeAlignment;
21     assert(approxEqual(q.identity, 0.8));
22     assert(approxEqual(q.similarity, 1.0));
23     assert(approxEqual(q.gapRate, 0.0));
24     
25 
26     auto q2 = p.sw_striped("GGCTTCTGATCAGGCTTCT", "GACTA");
27 
28     assert(q2.queryPosition == 7);
29     assert(q2.position == 0);
30     assert(q2.score == 14);
31     assert(q2.cigar.alignedLength == 5);
32     assert(q2.cigar.toString == "7S2=1D1=1I1=7S");
33     q2.writeAlignment;
34     assert(approxEqual(q2.identity, 0.667));
35     assert(approxEqual(q2.similarity, 0.667));
36     assert(approxEqual(q2.gapRate, 0.333));
37 
38     auto q3 = p.sw_striped("GACTAA", "GGCTTCTGATCAGGCTTCT");
39 
40     assert(q3.queryPosition == 0);
41     assert(q3.position == 7);
42     assert(q3.score == 14);
43     assert(q3.cigar.alignedLength == 5);
44     assert(q3.cigar.toString == "2=1D1=1I1=1S");
45     q3.writeAlignment;
46     assert(approxEqual(q3.identity, 0.667));
47     assert(approxEqual(q3.similarity, 0.667));
48     assert(approxEqual(q3.gapRate, 0.333));
49 
50     auto q4 = p.nw_scan("GATTA", "GACTA");
51 
52     assert(q4.queryPosition == 0);
53     assert(q4.position == 0);
54     assert(q4.score == 16);
55     assert(q4.cigar.alignedLength == 5);
56     assert(q4.cigar.toString == "2=1X2=");
57     q4.writeAlignment;
58     assert(approxEqual(q4.identity, 0.8));
59     assert(approxEqual(q4.similarity, 1.0));
60     assert(approxEqual(q4.gapRate, 0.0));
61 
62     auto q5 = p.nw_scan("GGCTTCTGATCAGGCTTCT", "GACTA");
63 
64     assert(q5.queryPosition == 0);
65     assert(q5.position == 0);
66     assert(q5.score == -14);
67     assert(q5.cigar.alignedLength == 5);
68     assert(q5.cigar.toString == "1=1X2=4I1=10S");
69     q5.writeAlignment;
70     writeln(q5.identity);
71     writeln(q5.similarity);
72     writeln(q5.gapRate);
73     assert(approxEqual(q5.identity, 0.444));
74     assert(approxEqual(q5.similarity, 0.556));
75     assert(approxEqual(q5.gapRate, 0.444));
76 
77     // TODO: reconcile issues where cigar manipulation for 
78     // use with bams changes output statistics
79 
80     // auto q52 = p.aligner!("nw","stats")("GGCTTCTGATCAGGCTTCT", "GACTA");
81 
82     // writeln(q52.queryPosition);
83     // writeln(q52.position);
84     // writeln(q52.score);
85     // assert(q52.queryPosition == 0);
86     // assert(q52.position == 0);
87     // assert(q52.score == -14);
88     // q5.writeAlignment;
89     // writeln(q52.identity);
90     // writeln(q52.similarity);
91     // writeln(q52.gapRate);
92     // assert(approxEqual(q5.identity, 0.444));
93     // assert(approxEqual(q5.similarity, 0.556));
94     // assert(approxEqual(q5.gapRate, 0.444));
95     // assert(approxEqual(q5.identity, 0.211));
96     // assert(approxEqual(q5.similarity, 0.263));
97     // assert(approxEqual(q5.gapRate, 0.737));
98 
99     auto q6 = p.nw_scan("GACTAA", "GGCTTCTGATCAGGCTTCT");
100     assert(q6.queryPosition == 0);
101     assert(q6.position == 0);
102     assert(q6.score == -8);
103     assert(q6.cigar.alignedLength == 12);
104     assert(q6.cigar.toString == "1=1X2=4D1=2D1=");
105     q6.writeAlignment;
106     writeln(q6.identity);
107     writeln(q6.similarity);
108     writeln(q6.gapRate);
109     assert(approxEqual(q6.identity, 0.4167));
110     assert(approxEqual(q6.similarity, 0.5));
111     assert(approxEqual(q6.gapRate, 0.5));
112 
113     // TODO: report potiential bug as 
114     // parasail similarity == identity 
115     // when it shoudl not be
116 
117     // auto q62 = p.aligner!("nw","stats")("GACTAA", "GGCTTCTGATCAGGCTTCT");
118     // writeln(q62.queryPosition);
119     // writeln(q62.position);
120     // writeln(q62.score);
121     // writeln(q62.identity);
122     // writeln(q62.similarity);
123     // writeln(q62.gapRate);
124     // assert(q6.queryPosition == 0);
125     // assert(q6.position == 0);
126     // assert(q6.score == -8);
127     // assert(q6.cigar.alignedLength == 12);
128     // assert(q6.cigar.toString == "1=1X2=4D1=2D1=");
129     // q6.writeAlignment;
130     // writeln(q6.identity);
131     // writeln(q6.similarity);
132     // writeln(q6.gapRate);
133     // assert(approxEqual(q6.identity, 0.4167));
134     // assert(approxEqual(q6.similarity, 0.5));
135     // assert(approxEqual(q6.gapRate, 0.5));
136     // assert(approxEqual(q6.identity, 0.263));
137     // assert(approxEqual(q6.similarity, 0.316));
138     // assert(approxEqual(q6.gapRate, 0.684));
139     
140 }