Исходный код программы поиска теоретических ПЦР-фрагментов по найденным праймерам с распараллеленным циклом на Си

 

/*fragparr.c*/

 

#include <stdio.h>

#include <string.h>

#include <signal.h>

#include <malloc.h>

#include <stdlib.h>

#include <sys/stat.h>

#include <sys/types.h>

#include <fcntl.h>

#include <sys/time.h>

#include <omp.h>

 

#define FILENAME1 "sequence.raw"

#define FILENAME2 "primersort0.raw"

#define FILENAME3 "optprimer.raw"

 

#define MIN(a,b) ((a) < (b)? (a): (b))

#define MAX(a,b) ((a) > (b)? (a): (b))

#define XSIZE 4200

#define SIGMA 256

#define WORD 32

 

/*void prof(char * str, int i, struct timeval * B){

struct itimerval A;

struct timeval tmp;

if (getitimer(ITIMER_PROF, &A)){

perror("getitimer");

exit(1);

}

tmp.tv_usec=B->tv_usec-A.it_value.tv_usec;

tmp.tv_sec=B->tv_sec-A.it_value.tv_sec;

if (tmp.tv_usec<0){

tmp.tv_usec+=1000000;

tmp.tv_sec--;

}

fprintf(stderr, "%d %s: %d.%d - %d.%d\n", i, str, A.it_value.tv_sec,

A.it_value.tv_usec,

tmp.tv_sec, tmp.tv_usec);

*B=A.it_value;

}*/

 

/*Инвертирование строки*/

void reverse(char p1[], char p2[]){

char * tmp=p1;

while (*tmp++);

for (tmp--; --tmp>=p1; *p2++=*tmp);

*p2='\0';

}

 

/*Инвертирование символов*/

void inverse(char p1[], char p2[], int strl){

int i;

for (i=0; i<strl; i++)

switch (p1[i]){

case 'A':

p2[i]='T';

break;

case 'T':

p2[i]='A';

break;

case 'C':

p2[i]='G';

break;

case 'G':

p2[i]='C';

break;

}

p2[i]='\0';

}

 

/*[1]-----Примитивный алгоритм*/

/*

char* hospool(char *y, char *x, int m){

char* yy;

int i, j;

yy=y;

i=0;

while (y[i]!= '\0'){

for (j=0; j < m; ++j){

if (y[i+j]!= x[j])

break;

}

if (j==m){

yy=y+i;

return (yy);

}

i++;

}

return NULL;

}

*/

 

/*[2]-----Knuth Morris-Pratt-----/kmp/-----1977*/

/*

char* kmpstr(char *y, char *x, int m){

char* yy;

int i, j, n;

yy=y;

int *protr =(int*)malloc(m*sizeof(int));

protr[0]=0;

// Построение таблицы смещений

for(i=1, j=0; i<m; i++){

while(j>0 && x[j]!=x[i])

j = protr[j-1];

if(x[j]==x[i])

j++;

protr[i]=j;

}

// Поиск

i=0;

j=0;

while (y[i]!= '\0'){

while(j>0 && x[j]!=y[i])

j=protr[j-1];

if(x[j]==y[i])

j++;

if (j==m){

free (protr);

yy=y+i-j+1;

return yy;

}

i++;

}

free (protr);

return NULL;

}

*/

 

/*[3]-----Boyer Moore-----/bm/-----1977*/

/*

//Построение таблицы для сдвига по функции плохого символа

void preBmBc(unsigned char *x, int m, int bmBc[]) {

int i;

for (i = 0; i < SIGMA; ++i)

bmBc[i] = m;

for (i = 0; i < m - 1; ++i)

bmBc[x[i]] = m - i - 1;

}

 

void suffixes(unsigned char *x, int m, int *suff) {

int f, g, i;

suff[m - 1] = m;

g = m - 1;

for (i = m - 2; i >= 0; --i) {

if (i > g && suff[i + m - 1 - f] < i - g)

suff[i] = suff[i + m - 1 - f];

else {

if (i < g)

g = i;

f = i;

while (g >= 0 && x[g] == x[g + m - 1 - f])

--g;

suff[i] = f - g;

}

}

}

 

//Построение таблицы для сдвига по функции хорошего суффикса

void preBmGs(unsigned char *x, int m, int bmGs[]) {

int i, j, suff[XSIZE];

suffixes(x, m, suff);

for (i = 0; i < m; ++i)

bmGs[i] = m;

j = 0;

for (i = m - 1; i >= 0; --i)

if (suff[i] == i + 1)

for (; j < m - 1 - i; ++j)

if (bmGs[j] == m)

bmGs[j] = m - 1 - i;

for (i = 0; i <= m - 2; ++i)

bmGs[m - 1 - suff[i]] = m - 1 - i;

}

 

char * search(unsigned char *y, unsigned char *x, int m) {

int i, j, bmGs[XSIZE], bmBc[SIGMA], n;

char * yy;

yy=y;

n=strlen(y);

 

// Подготовительные вычисления

preBmGs(x, m, bmGs);

preBmBc(x, m, bmBc);

 

// Поиск

j = 0;

while (j <= n - m) {

for (i = m - 1; i >= 0 && x[i] == y[i + j]; --i);

if (i < 0)

return yy+j;

else

j += MAX(bmGs[i], bmBc[y[i + j]] - m + 1 + i);

}

return NULL;

}

*/

 

/*[4]-----Turbo Boyer Moore-----/tbm/-----1994*/

/*

//Построение таблицы для сдвига по функции плохого символа

void preBmBc(unsigned char *x, int m, int bmBc[]) {

int i;

for (i = 0; i < SIGMA; ++i)

bmBc[i] = m;

for (i = 0; i < m - 1; ++i)

bmBc[x[i]] = m - i - 1;

}

 

void suffixes(unsigned char *x, int m, int *suff) {

int f, g, i;

suff[m - 1] = m;

g = m - 1;

for (i = m - 2; i >= 0; --i) {

if (i > g && suff[i + m - 1 - f] < i - g)

suff[i] = suff[i + m - 1 - f];

else {

if (i < g)

g = i;

f = i;

while (g >= 0 && x[g] == x[g + m - 1 - f])

--g;

suff[i] = f - g;

}

}

}

 

//Построение таблицы для сдвига по функции хорошего суффикса

void preBmGs(unsigned char *x, int m, int bmGs[]) {

int i, j, suff[XSIZE];

suffixes(x, m, suff);

for (i = 0; i < m; ++i)

bmGs[i] = m;

j = 0;

for (i = m - 1; i >= 0; --i)

if (suff[i] == i + 1)

for (; j < m - 1 - i; ++j)

if (bmGs[j] == m)

bmGs[j] = m - 1 - i;

for (i = 0; i <= m - 2; ++i)

bmGs[m - 1 - suff[i]] = m - 1 - i;

}

 

char * search(unsigned char *y, unsigned char *x, int m) {

int bcShift, i, j, shift, u, v, turboShift, bmGs[XSIZE], bmBc[SIGMA], n;

char * yy;

yy=y;

n=strlen(y);

 

// Подготовительные вычисления

preBmGs(x, m, bmGs);

preBmBc(x, m, bmBc);

 

// Поиск

j = u = 0;

shift = m;

while (j <= n - m) {

i = m - 1;

while (i >= 0 && x[i] == y[i + j]) {

--i;

if (u!= 0 && i == m - 1 - shift)

i -= u;

}

if (i < 0)

return yy+j;

else {

v = m - 1 - i;

turboShift = u - v;

bcShift = bmBc[y[i + j]] - m + 1 + i;

shift = MAX(turboShift, bcShift);

shift = MAX(shift, bmGs[i]);

if (shift == bmGs[i])

u = MIN(m - shift, v);

else {

if (turboShift < bcShift)

shift = MAX(shift, u + 1);

u = 0;

}

}

j += shift;

}

return NULL;

}

*/

 

/*[5]-----Zhu Takaoka-----/zt/-----1987*/

void suffixes(unsigned char *x, int m, int *suff) {

int f, g, i;

suff[m - 1] = m;

g = m - 1;

for (i = m - 2; i >= 0; --i) {

if (i > g && suff[i + m - 1 - f] < i - g)

suff[i] = suff[i + m - 1 - f];

else {

if (i < g)

g = i;

f = i;

while (g >= 0 && x[g] == x[g + m - 1 - f])

--g;

suff[i] = f - g;

}

}

}

 

//Построение таблицы для сдвига по функции хорошего суффикса

void preBmGs(unsigned char *x, int m, int bmGs[]) {

int i, j, suff[XSIZE];

suffixes(x, m, suff);

for (i = 0; i < m; ++i)

bmGs[i] = m;

j = 0;

for (i = m - 1; i >= 0; --i)

if (suff[i] == i + 1)

for (; j < m - 1 - i; ++j)

if (bmGs[j] == m)

bmGs[j] = m - 1 - i;

for (i = 0; i <= m - 2; ++i)

bmGs[m - 1 - suff[i]] = m - 1 - i;

}

 

//Построение таблицы для сдвига по функции плохого символа

void preZtBc(unsigned char *x, int m, int ztBc[SIGMA][SIGMA]) {

int i, j;

for (i = 0; i < SIGMA; ++i)

for (j = 0; j < SIGMA; ++j)

ztBc[i][j] = m;

for (i = 0; i < SIGMA; ++i)

ztBc[i][x[0]] = m - 1;

for (i = 1; i < m - 1; ++i)

ztBc[x[i - 1]][x[i]] = m - 1 - i;

}

 

char * search(unsigned char *y, unsigned char *x, int m) {

int i, j, ztBc[SIGMA][SIGMA], bmGs[XSIZE], n;

char * yy;

yy=y;

n=strlen(y);

 

// Подготовительные вычисления

preZtBc(x, m, ztBc);

preBmGs(x, m, bmGs);

 

// Поиск

j = 0;

while (j <= n - m) {

i = m - 1;

while (i >=0 && x[i] == y[i + j])

--i;

if (i < 0)

return yy+j;

else

j += MAX(bmGs[i],ztBc[y[j + m - 2]][y[j + m - 1]]);

}

return NULL;

}

 

/*[6]-----Genomic Rapid Algofor String Pattern Matching-----/graspm/-----2009*/

/*

typedef struct GRASPmList {

int k;

struct GRASPmList *next;

} GList;

 

void ADD_LIST(GList **l, int e) {

GList *t = (GList *)malloc(sizeof(GList));

t->k = e;

t->next = *l;

*l = t;

}

 

char * search(unsigned char *t, unsigned char *p, int m) {

GList *pos, *z[SIGMA];

int i, j, k, first, hbc[SIGMA], n;

char * yy;

yy=t;

n=strlen(t);

 

// Подготовительные вычисления

for(i=0; i<SIGMA; i++)

z[i]=NULL;

if (p[0] == p[m-1])

for(i=0; i<SIGMA; i++)

ADD_LIST(&z[i], 0);

for(i=0; i<m-1;i++)

if (p[i+1] == p[m-1])

ADD_LIST(&z[p[i]],(i+1));

for(i=0;i<SIGMA;i++)

hbc[i]=m;

for(i=0;i<m;i++)

hbc[p[i]]=m-i-1;

 

// Поиск

j = m-1;

while (j<n) {

while (k=hbc[t[j]]) j += k; {

pos = z[t[j-1]];

while(pos!=NULL) {

k = pos->k;

i = 0;

first = j-k;

while(i<m && p[i]==t[first+i])

i++;

if(i==m && first<=n-m)

return yy+first;

pos = pos->next;

}

}

j+=m;

}

return NULL;

}

*/

 

/*[7]-----Quick Search-----/qs/-----1990*/

/*

//Построение таблицы смещений

void preQsBc(unsigned char *P, int m, int qbc[]) {

int i;

for(i=0;i<SIGMA;i++)

qbc[i]=m+1;

for(i=0;i<m;i++)

qbc[P[i]]=m-i;

}

 

char * search(unsigned char *T, unsigned char *P, int m) {

int i, s, qsbc[SIGMA], n;

char * yy;

yy=T;

n=strlen(T);

 

// Подготовительные вычисления

preQsBc(P, m, qsbc);

 

// Поиск

s = 0;

while(s<=n-m) {

i=0;

while(i<m && P[i]==T[s+i])

i++;

if(i==m)

return yy+s;

s+=qsbc[T[s+m]];

}

return NULL;

}

*/

 

/*[8]-----Shift Or-----/so/-----1992*/

/*

int preSo(unsigned char *x, int m, unsigned int S[]) {

unsigned int j, lim;

int i;

for (i = 0; i < SIGMA; ++i)

S[i] = ~0;

for (lim = i = 0, j = 1; i < m; ++i, j <<= 1) {

S[x[i]] &= ~j;

lim |= j;

}

lim = ~(lim>>1);

return(lim);

}

 

char * search_large(unsigned char *y, unsigned char *x, int m) {

unsigned int lim, D, k, h, p_len;

unsigned int S[SIGMA];

int j, n;

char *yy;

yy=y;

n=strlen(y);

p_len = m;

m=32;

// Подготовительные вычисления

lim = preSo(x, m, S);

// Поиск

for (D = ~0, j = 0; j < n; ++j) {

D = (D<<1) | S[y[j]];

if (D < lim) {

k = 0;

h = j-m+1;

while(k<p_len && x[k]==y[h+k])

k++;

if(k==p_len)

return yy+j-m+1;

}

}

return NULL;

}

 

char * search(unsigned char *y, unsigned char *x, int m) {

unsigned int lim, D;

unsigned int S[SIGMA];

int j, n;

if (m > WORD) return search_large(y, x, m);

char *yy;

yy=y;

n=strlen(y);

// Подготовительные вычисления

lim = preSo(x, m, S);

// Поиск

for (D = ~0, j = 0; j < n; ++j) {

D = (D<<1) | S[y[j]];

if (D < lim)

return yy+j-m+1;

}

return NULL;

}

*/

 

/*[9]-----Karp Rabin-----/kr/-----1987*/

/*

#define REHASH(a, b, h) ((((h) - (a)*d) << 1) + (b))

 

char * search(char *y, char *x, int m) {

int d, hx, hy, i, j, n;

char *yy;

yy=y;

n=strlen(y);

// Подготовительные вычисления

for (d = i = 1; i < m; ++i)

d = (d<<1);

for (hy = hx = i = 0; i < m; ++i) {

hx = ((hx<<1) + x[i]);

hy = ((hy<<1) + y[i]);

}

// Поиск

j = 0;

while (j <= n-m) {

if (hx == hy && memcmp(x, y + j, m) == 0)

return yy+j;

hy = REHASH(y[j], y[j + m], hy);

++j;

}

return NULL;

}

*/

 

struct write_rec{

int offset;

short len;

char lenp;

char num;

char prim[24];

};

 

/*Для ассемблерной реализации*/

/*char* mystrstr(char *x, char *y);*/

 

int main(int args, char ** argv){

int f, k, flag, i, lenf, j, kolvo, cpu;

char * bufdata, * primer, * p, * penter, * pend, * pp, *p2enter, *ptester;

struct stat st;

struct itimerval A;

struct timeval B;

struct write_rec W;

A.it_interval.tv_sec=0x7fffffff;

A.it_interval.tv_usec=0;

A.it_value=A.it_interval;

B=A.it_value;

signal(SIGPROF, SIG_IGN);

if (setitimer(ITIMER_PROF, &A, NULL)){

perror("setittime");

exit(1);

}

 

if ((f=open(FILENAME1, O_RDONLY))==-1){

perror("open");

exit(1);

}

if (fstat(f, &st)==-1){

perror("stat");

exit(1);

}

if ((bufdata=malloc(st.st_size+1))==NULL){

perror("malloc");

exit(1);

}

memset(bufdata, '\0', st.st_size+1);

if (read(f, bufdata, st.st_size)!=st.st_size){

perror("read");

exit(1);

}

close(f);

if ((f=open(FILENAME2, O_RDONLY))==-1){

perror("open");

exit(1);

}

if (fstat(f, &st)==-1){

perror("stat");

exit(1);

}

if ((primer=malloc(st.st_size+1))==NULL){

perror("malloc");

exit(1);

}

memset(primer, '\0', st.st_size+1);

if (read(f, primer, st.st_size)!=st.st_size){

perror("read");

exit(1);

}

lenf=st.st_size;

close(f);

if ((f=open(FILENAME3, O_WRONLY | O_APPEND | O_CREAT, 0600))==-1){

perror("openwrite");

exit(1);

}

// prof("Before out for", 0, &B);

printf("Длина %d\n", lenf);

W.num=1;

cpu=atoi(argv[1]);

printf("Количество потоков %d\n", cpu);

omp_set_num_threads(cpu);

 

#pragma omp parallel for private(penter), private(pend), private(pp), private(p)

for (j=0; j<lenf; j=j+32){

char str[26]="ATG", rev[26], inv[26];

p=primer+j;

flag=0;

strcpy(str+3, p);

// prof("Before strlen", i, &B);

int l=strlen(str);

W.lenp=l;

memset(W.prim, 0, 24);

strncpy(W.prim, str+3, l);

// prof("Before strlen", i, &B);

reverse(str, rev);

inverse(rev, inv, l);

/*Первое вхождение праймера в строку*/

// prof("Before strstr 1", i, &B);

for (p=bufdata; penter=search(p, str, l); p=pend+l){

/*for (p=bufdata; penter=strstr(p, str); p=pend+l){*/

/*for (p=bufdata; penter=mystrstr(p, str); p=pend+l){*/

// prof("Before strstr 2", i, &B);

/*Первое вхождение инвертированного праймера в строку*/

pend=search(penter+l, inv, l);

/*pend=strstr(penter+l, inv);*/

/*pend=mystrstr(penter+l, inv);*/

// prof("After strstr 2", i, &B);

if (pend){

char ctmp;

ctmp=*pend;

*pend='\0';

/*Поиск дополнительных вхождений праймера в строку*/

// prof("Before strstr 3", i, &B);

for (pp=penter; pp; pp=search(pp+1, str, l))

/*for (pp=penter; pp; pp=strstr(pp+1, str))*/

/*for (pp=penter; pp; pp=mystrstr(pp+1, str))*/

penter=pp;

// prof("After strstr 3", i, &B);

*pend=ctmp;

k=pend-penter+l;

// prof("In", i, &B);

if (k>60 && k<3500){

if (flag==0){

printf("Прямой праймер %s\n", str);

printf("Обратный праймер %s\n", rev);

printf("Обратный инвертированный праймер %s\n", inv);

flag=1;

}

W.offset=penter-bufdata+1;

W.len=k;

if (write(f, &W, 32)!=32){

perror("write");

exit(1);

}

printf("%d %d %d\n", penter-bufdata+1, k, l);

}

}

else

break;

}

 

}

free(bufdata);

free(primer);

close(f);

printf("Поиск завершен!\n");

exit(0);

}