[Bio] / Kmers2 / rna_seq.c Repository:
ViewVC logotype

View of /Kmers2/rna_seq.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (download) (as text) (annotate)
Wed Apr 9 16:58:43 2014 UTC (5 years, 7 months ago) by overbeek
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +5 -2 lines
fixed initialization

#include <unistd.h>
#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include  <sys/types.h>
#include  <sys/ipc.h>
#include  <stdio.h>
#include  <string.h>
#include <fcntl.h>
#include <sys/stat.h>
#include <sys/wait.h>
#include <sys/mman.h>
#include <errno.h>   /* errno */

#define READSZ 4*BUFSIZ

int K = 24;  /* cannot exceed 39 */
int NUM_BUFFERS = 20;

#define WU_PER_BUFFER 20000
#define MAX_READ_SIZE 2000
#define SEQ_CHUNKS_PER_BUFFER 5000
#define MAX_SZ_IDS_IN_BUFFER 20000

char trans[256] = {4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,0,4,1,4,4,4,2,4,4,4,4,4,4,4,4,4,4,4,4,3,3,4,4,4,4,4,4,4,4,4,4,4,0,4,1,4,4,4,2,4,4,4,4,4,4,4,4,4,4,4,4,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4};

char transC[256] = {4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,3,4,2,4,4,4,1,4,4,4,4,4,4,4,4,4,4,4,4,0,0,4,4,4,4,4,4,4,4,4,4,4,3,4,2,4,4,4,1,4,4,4,4,4,4,4,4,4,4,4,4,0,0,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4};

typedef struct raw_output {
  int h1;
  int h2;
  int bad_input;
} raw_t;

typedef struct fasta_mapping {
  char *id;                   /* ID array */
  int   fastaB;               /* -1 -> empty */
  int   bufferB;
  int   ln;
} fasta_mapping_t;

typedef struct buffer {
  char seq[WU_PER_BUFFER+40];   /* we include a K-1 character pad  */
  int  data_in_buffer;
  raw_t raw[WU_PER_BUFFER];
  fasta_mapping_t where_seq[SEQ_CHUNKS_PER_BUFFER];
  char ids_in_buffer[MAX_SZ_IDS_IN_BUFFER];
} buffer_t;

#define HASH2 15485867
#define HSIZE 100000007

#define ID_TAB_MAX_ENTRIES 100000
#define ID_TAB_SZ 1000000

typedef struct hash_entry {
  int kmer_hash2;
  int id_index;             
  int kmer_position;   /* strand: positive for plus; -(p+(K-1)) for minus */
} hash_entry_t;

typedef struct handle {
  char buffer[READSZ];
  char *buffP;
  int state;
  int left;
} handle_t;
  
int fill (handle_t *h) {
  h->left = read (STDIN_FILENO, h->buffer, READSZ);
  h->buffP = h->buffer;
  return h->left;
}
	
handle_t *get_handle() {
  handle_t *h = (handle_t *) malloc(sizeof(handle_t));

  char c;
  if ((read(STDIN_FILENO,&c,1) < 1) || (c != '>')) {
	fprintf (stderr, "Input does not start with >\n");
	exit(0);
  }

  h->left = 0;
  h->state    = 1;

  return h;
}

int get_id(handle_t *h, char *id) {
  char c;
  int got = 0;

  while (!got) {
    while (h->left  && ((c = *(h->buffP)) != ' ') && (c != '\t')  && (c != '\n') ) {
      *id++ = c;
      h->left--;
      h->buffP++;
    }	

    if (h->left == 0) {
      if (!fill(h)) {
	return 0;
      }
    } else {
      int got2 = 0;
      while (! got2) {
	while (h->left  && ((c = *(h->buffP)) != '\n') ) {
	  h->left--;
	  h->buffP++;
	}
	if (h->left == 0) {
	  if (!fill(h)) {
	    return 0;
	  } 
	} else {
	  h->left--;
	  h->buffP++;
	  got2 = 1;
	}
      }
      *id = '\0';
      got = 1;	
    }
  }
  return 1;
}

int get_seq(handle_t *h, char *where, int want, int *how_many) {

  *how_many = 0;
  int got = 0;
  char c;
  while (! got) {
    while (want && h->left  && ((c = *(h->buffP)) != '>')) {
	if ((c != ' ') && (c != '\t') && (c != '\n')) {
	  *where++ = *(h->buffP);
	  want--;
	  (*how_many)++;
	}
	h->left--;
	h->buffP++;
    }
    if (want == 0) {
      return 0;
    }
    else {
      if (h->left == 0) {
	if (! fill(h)) {
	  return 1;
	} 
      }
      else {
	h->left--;
	h->buffP++;
	return 1;
      }
    }
  }
  return 0;
}    

int dna_char(c)
char c;
{
    switch (c)
    {
      case 'a':
      case 'A':
	return 0;

      case 'c':
      case 'C':
	return 1;

      case 'g':
      case 'G':
	return 2;

      case 't':
      case 'u':
      case 'T':
      case 'U':
	return 3;

      default:
	return 4;;
    }
}

char compl(c)
char c;
{
    switch (c)
    {
      case 'a':
	return 't';
      case 'A':
	return 'T';

      case 'c':
	return 'g';
      case 'C':
	return 'G';

      case 'g':
	return 'c';
      case 'G':
	return 'C';

      case 't':
      case 'u':
	return 'a';
      case 'T':
      case 'U':
	return 'A';

      case 'm':
	return 'k';
      case 'M':
	return 'K';

      case 'r':
	return 'y';
      case 'R':
	return 'Y';

      case 'w':
	return 'w';
      case 'W':
	return 'W';

      case 's':
	return 'S';
      case 'S':
	return 'S';

      case 'y':
	return 'r';
      case 'Y':
	return 'R';

      case 'k':
	return 'm';
      case 'K':
	return 'M';

      case 'b':
	return 'v';
      case 'B':
	return 'V';

      case 'd':
	return 'h';
      case 'D':
	return 'H';

      case 'h':
	return 'd';
      case 'H':
	return 'D';

      case 'v':
	return 'b';
      case 'V':
	return 'B';

      case 'n':
	return 'n';
      case 'N':
	return 'N';

      default:
	return c;
    }
}

raw_t encode_plus(char *kmer) {
  long long v = 0;
  int i;
  raw_t r;

  r.bad_input = 0;
  for (i = 0; i < K; i++) {
    int c = trans[*(kmer+i)];
    if (c == 4) {
      r.bad_input = 1;;
    }
    v = (v << 2) | c;
  }
  r.h1 = v % HSIZE;
  r.h2 = v % HASH2;
  return r;
}

raw_t encode_minus(char *kmer) {
  long long v = 0;
  int i;
  raw_t r;

  r.bad_input = 0;
  for (i = K-1; i >= 0; i--) {
    int c = transC[*(kmer+i)];
    if (c == 4) {
      r.bad_input = 1;;
    }
    v = (v << 2) | c;

  }
  r.h1 = v % HSIZE;
  r.h2 = v % HASH2;
  return r;
}

void load_kmers(char *filename, hash_entry_t **hash_tab, char ***id_index_ar) {

  *hash_tab = malloc(HSIZE * sizeof(hash_entry_t));
  int i;
  hash_entry_t *p1;
  for (i=0; (i < HSIZE); i++) {
    p1 = *hash_tab + i;
    p1->kmer_position = 0;
  }
    
  *id_index_ar    = malloc(ID_TAB_MAX_ENTRIES * sizeof(char *));
  char **p_index  = *id_index_ar;

  char *vals      = malloc(ID_TAB_SZ);
  char *p         = vals;
  FILE *ifp      = fopen(filename,"r");
  if (ifp == NULL) { 
    printf("could not open %s\n",filename);
    exit(1);
  }

  int sz = 0;
  int j;
  while ((fscanf(ifp,"%d\t",&j) == 1) && fgets(p,1000,ifp)) {
    if (sz != j) {
      fprintf(stderr,"Your index must be dense and in order (see line %ld, should be %d)\n",p-vals,sz);
      exit(1);
    }
    p_index[sz] = p;
    
    p += strlen(p) -1;
    *(p++) = '\0';
    if ((sz >= ID_TAB_MAX_ENTRIES) || ((p-vals) > (ID_TAB_SZ - 1000))) {
      fprintf(stderr,"Your contig id array is too small; bump ID_TAB_MAX_ENTRIES and ID_TAB_SZ\n");
      exit(1);
    }
    sz += 1;
  }

  int c;
  while ((c = getc(ifp)) != '\n') {}

  char kmer_string[K+1];
  int id;
  int pos;
  while (fscanf(ifp,"%s\t%d\t%d",kmer_string,&id,&pos) == 3) {
    raw_t v = encode_plus(kmer_string);
    if (! v.bad_input) {
      hash_entry_t *hpos = *hash_tab + v.h1;
      if (hpos->kmer_position == 0) {
	hpos->id_index      = id;
	hpos->kmer_position = pos;
	hpos->kmer_hash2  = v.h2;
      }
      else {
	/* fprintf(stderr,"clash: losing %s\n",kmer_string); */
      }
      v = encode_minus(kmer_string);
      hpos = *hash_tab + v.h1;
      if (hpos->kmer_position == 0) {
	hpos->id_index      = id;
	hpos->kmer_position = -pos; 
	hpos->kmer_hash2  = v.h2;
      }
      else {
	/* fprintf(stderr,"clash: losing comlement of %s\n",kmer_string); */
      }
    }
  }
}

buffer_t *buffer_loc;

void allocate_buffer() {

   size_t sz = sizeof(buffer_t);
   buffer_loc = (buffer_t *) malloc(sz);
}

void print_usage();

char *kmers;
char *input;
char *output;

hash_entry_t *ref_hash;
char         **ref_id_index_ar;

int read_input(handle_t *h) {
  buffer_t *bp = buffer_loc;
  int fastaII = 0;
  bp->data_in_buffer = 0;
  char *next_id = bp->ids_in_buffer;

  int got = 0;
  while ((bp->data_in_buffer < (WU_PER_BUFFER - MAX_READ_SIZE)) &&
	 get_id(h,next_id) &&
	 (get_seq(h,bp->seq+bp->data_in_buffer,
		  (WU_PER_BUFFER+(K-1) - bp->data_in_buffer),
		  &(bp->where_seq[fastaII].ln)))) {

    got++;
    bp->where_seq[fastaII].id = next_id;
    bp->where_seq[fastaII].fastaB = 0;
    bp->where_seq[fastaII].bufferB = bp->data_in_buffer;
    bp->where_seq[fastaII+1].fastaB = -1;
    next_id += strlen(next_id)+1;
    bp->data_in_buffer += bp->where_seq[fastaII].ln;
    if (bp->where_seq[fastaII].ln > K) {
	bp->where_seq[fastaII].ln -= (K-1);
    }
    else {
      fprintf(stderr,"BUG\n");
    }
    
    fastaII++;
    bp->where_seq[fastaII].fastaB = -1;
  }
  return got;
}

void work() {
  buffer_t *bp = buffer_loc;
  int i;
  for (i=0; (i < bp->data_in_buffer);i++) {
    bp->raw[i] = encode_plus(bp->seq+i);
  }
}

/* ==============  GLOBAL: hit recording =================== */ 
/* We consider a string of hits to begin with two against the same
   contig, with a close offset.
*/

char id_of_read[300];

typedef struct hit {
  int contig_id_index;
  int contig_offset;
  int strand;
  char * read_id;
  int read_offset;
} hit_t;

hit_t first;
hit_t last;
hit_t possible_challenger;

int num_hits;

/* ========= end hit global ========== */

void init_hit_list(char *read_id) {
  if (read_id == NULL) { fprintf(stderr,"BAD\n"); }
  strcpy(id_of_read,read_id);
  possible_challenger.contig_id_index = -1;
  num_hits = 0;
}

void flush_hit_list() {
  if (num_hits > 1) {
    if (first.contig_id_index == -1) {
      fprintf(stderr,"BAD BAD\n"); exit(1);
    }
    if (first.strand == 0) {
      fprintf(stdout,"%d\t%d\t%s\t%d\t%d\t%d\t%d\t%d\n",
	      num_hits,
	      first.contig_id_index,
	      first.read_id,
	      first.strand,
	      first.contig_offset+1,
	      last.contig_offset+K,
	      first.read_offset+1,
	      last.read_offset+K);
    }
    else {
      fprintf(stdout,"%d\t%d\t%s\t%d\t%d\t%d\t%d\t%d\n",
	      num_hits,
	      first.contig_id_index,
	      first.read_id,
	      first.strand,
	      first.contig_offset+((K-1)+1),
	      last.contig_offset+1,
	      first.read_offset+1,
	      last.read_offset+K);
    }
    num_hits = 0;
  }
}

void record_hit(int off_into_read,int strand,int contig_idI,int off_into_contig) {
  /*  fprintf(stderr,"HIT: %s %d %d %d\n",id_of_read,off_into_read,strand,off_into_contig);*/
  if (num_hits == 0) {
    first.contig_id_index = contig_idI;
    first.contig_offset = last.contig_offset = off_into_contig;
    first.strand = strand;
    first.read_id = id_of_read;
    first.read_offset = off_into_read;
    num_hits = 1;
  }
  else {
    int read_gap = abs(off_into_read - first.read_offset);
    if ((first.contig_id_index == contig_idI) && 
	(first.strand == strand) &&
	(abs(last.contig_offset - off_into_contig) < (read_gap + 20)) &&
	(abs(first.contig_offset - off_into_contig) < (read_gap + 20)) &&
	(first.read_id == id_of_read)) {
      last.contig_offset = off_into_contig;
      last.read_offset   = off_into_read;
      num_hits++;
      possible_challenger.contig_id_index = -1;
    }
    else {
      if (possible_challenger.contig_id_index == -1) {
	possible_challenger.contig_id_index = contig_idI;
	possible_challenger.contig_offset = off_into_contig;
	possible_challenger.strand = strand;
	possible_challenger.read_id = id_of_read;
	possible_challenger.read_offset = off_into_read;
      }
      else {
	flush_hit_list();
	first = possible_challenger;
	num_hits = 1;
	possible_challenger.contig_id_index = -1;
	record_hit(off_into_read,strand,contig_idI,off_into_contig);
      }
    }
  }
}

void process_raw(int p1, int p2, int fasta_off) {
  hash_entry_t *t = ref_hash+p1;
  
  if (p2 == t->kmer_hash2) {
    if (t->kmer_position < 0) {
      record_hit(fasta_off,1,t->id_index,-(t->kmer_position));
    }
    else {
      record_hit(fasta_off,0,t->id_index,t->kmer_position);
    }
  }
}

void write_output() {

  buffer_t *bp = buffer_loc;
  fasta_mapping_t *fm = bp->where_seq;
  while (fm->fastaB != -1) {
    init_hit_list(fm->id);
    raw_t *outP = bp->raw+fm->bufferB;
    int fasta_off = fm->fastaB;
    int entries_to_process = fm->ln;
    while (entries_to_process) {
      if (! outP->bad_input) {
	process_raw(outP->h1,outP->h2,fasta_off);       /* strand: 0 -> +, 1 -> - */
      }
      entries_to_process--;
      outP++;
      fasta_off++;
    }
    flush_hit_list();
    fm++;
  }
}


int main(int argc, char *argv[]) {
    int option = 0;

    /*
	  s - the name of the file containing the contig index table and the
	      signature kmers for the reference genome

	  K - the Kmer size
    */
    while ((option = getopt(argc, argv,"s:K:")) != -1) {
        switch (option) {
	     case 's' : kmers = optarg;  
                 break;
	     case 'K' : K = atoi(optarg);
                 break;
             default: print_usage(); 
                 exit(EXIT_FAILURE);
        }
    }

    allocate_buffer();
    load_kmers(kmers,&ref_hash,&ref_id_index_ar);
    handle_t *h = get_handle();
    while (read_input(h)) {
      work();
      write_output();
    }
    exit(0);
}	

void print_usage() {
    printf("Usage: rna_seq -s signature-kmers -K kmer-size\n");
}

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3