-
Notifications
You must be signed in to change notification settings - Fork 1
/
fl.cpp
104 lines (81 loc) · 3.66 KB
/
fl.cpp
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
#include "fl.hpp"
#include "adapters.hpp"
#include "utils.hpp"
#include "polya.hpp"
#include <algorithm>
#include <future>
#include <mutex>
read_set_t get_fl_reads(const read_set_t &reads, std::string a5, std::string a3, bool is_fq, int n_threads) {
read_set_t fl_reads;
std::mutex mu;
std::vector<std::future<void>> tasks;
for (int t = 0; t < n_threads; ++t) {
tasks.emplace_back(std::async(std::launch::async, [t, &mu, &reads, n_threads, &fl_reads, a5, a3, is_fq] {
for (int i = t; i < reads.size(); i += n_threads) {
auto read = reads[i];
// TODO: analyze small reads as well
if (read.seq.size() < 151) {
continue;
}
std::string seq = "";
std::string qt = "";
int pos5 = adapter_pos_start(read.seq, a5);
int pos3 = adapter_pos_end(read.seq, a3);
if (pos5 != -1 && pos3 != -1) {
seq = read.seq.substr(pos5, read.seq.size()-150+pos3-pos5);
if (is_fq) {
qt = read.quality.substr(pos5, read.quality.size()-150+pos3-pos5);
}
} else {
pos5 = adapter_pos_start(read.seq, reverse_complement(a5));
pos3 = adapter_pos_end(read.seq, reverse_complement(a3));
if (pos5 != -1 && pos3 != -1) {
seq = read.seq.substr(pos5, read.seq.size()-150+pos3-pos5);
if (is_fq) {
qt = read.quality.substr(pos5, read.quality.size()-150+pos3-pos5);
}
}
}
if (seq != "") {
auto polya_info = get_polya(seq);
if (polya_info.tail_pos > -1) {
seq = seq.substr(0, polya_info.tail_pos);
read_t fl_read;
fl_read.header = read.header;
fl_read.seq = seq;
if (is_fq) {
qt = qt.substr(0, polya_info.tail_pos);
fl_read.quality = qt;
fl_read.ann = read.ann;
}
std::lock_guard<std::mutex> lock(mu);
fl_reads.push_back(fl_read);
} else {
seq = reverse_complement(seq);
if (is_fq) {
std::reverse(qt.begin(), qt.end());
}
polya_info = get_polya(seq);
if (polya_info.tail_pos > -1) {
seq = seq.substr(0, polya_info.tail_pos);
read_t fl_read;
fl_read.header = read.header;
fl_read.seq = seq;
if (is_fq) {
qt = qt.substr(0, polya_info.tail_pos);
fl_read.quality = qt;
fl_read.ann = read.ann;
}
std::lock_guard<std::mutex> lock(mu);
fl_reads.push_back(fl_read);
}
}
}
}
}));
}
for (auto &&task : tasks) {
task.get();
}
return fl_reads;
}