This code is meant to compute the reverse complement for a given sequence of bases. The code is tested versus some test cases and seems to work so far.
The primary goal is that it should be as fast as possible (e.g. recovery from "file not found" or format errors in the input file is not required).
The code is pretty fast already. But I would like to ask if someone sees any obvious opportunities for a worthwhile optimisation or obvious mistakes of course.
#include <iostream>
#include <chrono>
#include <cstdio>
#include <sys/stat.h>
#include <thread>
#include <queue>
#include <tuple>
#include <algorithm>
#define INPUTFILENAME "~/70000000.txt"
#define OUTPUTFILENAME "~/output.txt"
namespace helpers {
size_t getFileSize(void)
{
struct stat stat_buf;
int rc = stat(INPUTFILENAME, &stat_buf);
return rc == 0 ? stat_buf.st_size : -1;
}
}
char* static_front = 0;
char* static_back = 0;
std::mutex mutex;
std::condition_variable cv;
std::queue<std::tuple<char*, char*>> queue;
static inline void staticTransformation(char* front, char *back) {
constexpr char complement[] = {
['A']='T', ['a']='T',
['C']='G', ['c']='G',
['G']='C', ['g']='C',
['T']='A', ['t']='A',
['U']='A', ['u']='A',
['M']='K', ['m']='K',
['R']='Y', ['r']='Y',
['W']='W', ['w']='W',
['S']='S', ['s']='S',
['Y']='R', ['y']='R',
['K']='M', ['k']='M',
['V']='B', ['v']='B',
['H']='D', ['h']='D',
['D']='H', ['d']='H',
['B']='V', ['b']='V',
['N']='N', ['n']='N'};
char help;
while(front < back) {
if(*front=='\n') {
front++;
}
if(*back=='\n') {
back--;
}
// swap after complement
help = complement[(*front)];
*front = complement[(*back)];
*back = help;
front++;
back--;
}
}
void staticTransformationWorker(int name) {
char* _front;
char* _back;
auto tuple = std::tie(_front, _back);
while(1) {
{
std::unique_lock<std::mutex> lock(mutex);
cv.wait(lock, []{ return !queue.empty(); });
tuple = queue.front();
queue.pop();
}
staticTransformation(_front, _back);
}
}
int main(int argc, const char * argv[]) {
size_t buffer_size = 0;
auto start = std::chrono::steady_clock::now();
unsigned concurentThreadsSupported = std::thread::hardware_concurrency();
unsigned numberOfThreads = std::min(concurentThreadsSupported,8u);
std::thread threads[8];
for(unsigned i=0; i<numberOfThreads; i++)
threads[i] = std::thread( staticTransformationWorker,i);
buffer_size = helpers::getFileSize();
char* buffer = static_cast<char*>(malloc(buffer_size));
// Read the file and close
FILE *file = fopen(INPUTFILENAME,"r");
fread(buffer, sizeof(char), buffer_size, file);
fclose(file);
char* sequence_start = 0; // start of the plain base codes, excl. header and first nl
char* sequence_end = 0; // end of the plain base codes, eccl. the last nl
auto tuple = std::tie(sequence_start, sequence_end);
for(char* cursor = buffer; cursor < buffer+buffer_size; cursor++) {
if(*cursor == '>') { // found header of a sequence
if(sequence_end > sequence_start) { // Sequenz fertig gecodet
{
std::lock_guard<std::mutex> lock(mutex);
queue.push(tuple);
}
cv.notify_all();
}
while(*cursor != '\n') cursor++; // skip the header
cursor++; // skip the newline at the end of the header
sequence_start = cursor; // store address of first letter code
}
if(*cursor == '\n') { // potential end of a sequence
sequence_end = cursor-1;
if(*(cursor+1)!='>') // when not a sequence header
while(*cursor != '\n') cursor++; // skip the line
}
}
// Reverse complement transformation in the while-loop is only triggered when a new sequence starts
// the last transformation won´t be triggered. We can run ist synchronous now, since this task has nothing to do.
staticTransformation(sequence_start,sequence_end);
file = fopen(OUTPUTFILENAME,"w");
fwrite(buffer, sizeof(char), buffer_size, file);
fclose(file);
for(unsigned i=0; i<numberOfThreads; i++)
threads[i].detach();
free(buffer);
auto duration = std::chrono::duration_cast< std::chrono::milliseconds>
(std::chrono::steady_clock::now() - start);
std::cout << "\n\nZeit: " << duration.count() << " ms" << std::endl;
return 0;
}
a typical input file would be something like: (much larger of course)
>ONE Homo sapiens alu
GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGGGAGGCCGAGGCGGGCGGA
CGGGCGTGGTGGCGCGCGCCTGTAATCCCAGCTACTCGGGAGGCTGAGGCAGGAGAATCG
CTTGAACCCGGGAGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCCAGCCTG
GGCGACAGAGCGAGACTCCG
One annoying thing is, for example, that running it several times directly in a row gets slower and slower results.
fwrite is about 40.5% and fread is 16% of overall execution time.
Here is the testfile I´m using: http://privatfrickler.de/70000000.txt
