code_c/Maasha/src/tetra_count.c

   1 #include "common.h"
   2 #include "filesys.h"
   3 #include "mem.h"
   4 #include "seq.h"
   5 #include "fasta.h"
   6
   7 #define BLOCK32     15
   8 #define ARRAY_SIZE ( 1 << ( BLOCK32 * 2 ) )
   9
  10 #define BLOCK_SIZE   4
  11 #define MAX_SPACE  256
  12 #define T            3  /* 11 on the rightmost two bits of bin. */
  13 #define C            1  /* 01 on the rightmost two bits of bin. */
  14 #define G            2  /* 10 on the rightmost two bits of bin. */
  15
  16
  17 uint mask_create( int size );
  18 void block_count( char *path, uint **array_ppt, seq_entry **entry_ppt, uint mask );
  19
  20
  21 int main( int argc, char *argv[] )
  22 {
  23     uint       mask  = 0;
  24     int        i     = 0;
  25     uint      *array = NULL;
  26     seq_entry *entry = NULL;
  27
  28     mask  = mask_create( BLOCK_SIZE );
  29
  30     array = mem_get_zero( sizeof( uint ) * ARRAY_SIZE );
  31
  32     seq_new( &entry, MAX_SEQ_NAME, MAX_SEQ );
  33
  34     printf( "mask: %d\n", mask );
  35
  36     for ( i = 1; i < argc; i++ ) {
  37         block_count( argv[ i ], &array, &entry, mask );
  38     }
  39
  40     return EXIT_SUCCESS;
  41 }
  42
  43
  44 uint mask_create( int size )
  45 {
  46     /* Martin A. Hansen, June 2008 */
  47
  48     /* Create a bit mask. */
  49
  50     uint i    = 0;
  51     uint mask = 0;
  52
  53     for ( i = 0; i < size; i++ )
  54     {
  55         mask <<= 2;
  56
  57         mask |= 3; /* add 11 to mask */
  58     }
  59
  60     return mask;
  61 }
  62
  63
  64 void block_count( char *path, uint **array_ppt, seq_entry **entry_ppt, uint mask )
  65 {
  66     /* Martin A. Hansen, August 2008 */
  67
  68     /* Scan a FASTA file for ... */
  69
  70     uint      *array = *array_ppt;
  71     seq_entry *entry = *entry_ppt;
  72     FILE      *fp    = NULL;
  73     uint       i     = 0;
  74     uint       j     = 0;
  75     uint       bin1  = 0;
  76     uint       bin2  = 0;
  77
  78     fp = read_open( path );
  79
  80     fprintf( stderr, "Processing: %s\n", path );
  81
  82     while ( ( fasta_get_entry( fp, &entry ) ) != 0 )
  83     {
  84         fprintf( stderr, "   Counting blocks in: %s ... ", entry->seq_name );
  85
  86         bin1 = 0;
  87         j    = 0;
  88
  89         for ( i = 0; entry->seq[ i ]; i++ )
  90         {
  91             bin <<= 2;
  92
  93             switch( entry->seq[ i ] )
  94             {
  95                 case 'A': case 'a':           j++; break;
  96                 case 'T': case 't': bin |= T; j++; break;
  97                 case 'C': case 'c': bin |= C; j++; break;
  98                 case 'G': case 'g': bin |= G; j++; break;
  99                 default: bin = 0; j = 0; break;
 100             }
 101
 102             if ( j >= BLOCK32 )
 103             {
 104                 array[ ( bin & mask ) ]++;
 105 /*
 106                 printf( "\n" );
 107                 printf( "mask          : %s\n", bits2string( mask ) );
 108                 printf( "bin           : %s\n", bits2string( bin ) );
 109                 printf( "bin     & mask: %s\n", bits2string( bin & mask ) );
 110 */
 111             }
 112         }
 113
 114         fprintf( stderr, "done.\n" );
 115     }
 116
 117     fprintf( stderr, "done.\n" );
 118
 119     close_stream( fp );
 120 }