code_c/Maasha/src/repeat-O-matic.c

   1 /*
   2     Copyright (C) 2008, Martin A. Hansen
   3
   4     This program determines the repetiveness of a genome by determining
   5     the number of identical 15-mers for each position in the genome.
   6
   7     The output is a fixedStep file ala the phastCons files from the UCSC
   8     Genome browser.
   9
  10     It is very fast and efficient using less than 8 Gb of memory to
  11     complete the human genome in roughly 30 minutes.
  12 */
  13
  14 #include "common.h"
  15 #include "mem.h"
  16 #include "filesys.h"
  17 #include "seq.h"
  18 #include "fasta.h"
  19
  20 // #define OLIGO_SIZE 15
  21 #define OLIGO_SIZE 5
  22 #define SIZE ( 1 << ( OLIGO_SIZE * 2 ) )
  23
  24 #define UINT_BITS 32
  25 #define T    3                            /* 11 on the rightmost two bits of bin. */
  26 #define C    1                            /* 01 on the rightmost two bits of bin. */
  27 #define G    2                            /* 10 on the rightmost two bits of bin. */
  28
  29 uint *oligo_count( char *path );
  30 uint  mask_create( int oligo_size );
  31 void  oligo_count_output( char *path, uint *array );
  32 void  fixedstep_put_entry( char *chr, int beg, int step_size, uint *block_array, int block_size );
  33
  34 int main( int argc, char *argv[] )
  35 {
  36     char *path  = NULL;
  37     uint *array = NULL;
  38
  39     path  = argv[ 1 ];
  40
  41     array = oligo_count( path );
  42
  43     //oligo_count_output( path, array );
  44
  45     return 0;
  46 }
  47
  48
  49 uint *oligo_count( char *path )
  50 {
  51     /* Martin A. Hansen, June 2008 */
  52
  53     /* Count the occurence of all oligos of a fixed size in a FASTA file. */
  54
  55     uint       *array   = NULL;
  56     uint        i       = 0;
  57     uint        mask    = 0;
  58     uint        bin     = 0;
  59     uint        bin_rc1 = 0;
  60     uint        bin_rc2 = 0;
  61     uint        j       = 0;
  62     uint        A_rc    = ( 3 << ( UINT_BITS - 2 ) );   /* 11 on the leftmost two bits an uint. */
  63     uint        G_rc    = ( 2 << ( UINT_BITS - 2 ) );   /* 10 on the leftmost two bits an uint. */
  64     uint        C_rc    = ( 1 << ( UINT_BITS - 2 ) );   /* 01 on the leftmost two bits an uint. */
  65     seq_entry  *entry   = NULL;
  66     FILE       *fp      = NULL;
  67
  68     array = mem_get_zero( sizeof( uint ) * SIZE );
  69
  70     mask  = mask_create( OLIGO_SIZE );
  71
  72     entry = mem_get( sizeof( entry ) );
  73
  74     fp = read_open( path );
  75
  76     while ( ( fasta_get_entry( fp, &entry ) ) )
  77     {
  78         fprintf( stderr, "Counting oligos in: %s ... ", entry->seq_name );
  79
  80         bin     = 0;
  81         bin_rc1 = 0;
  82         j       = 0;
  83
  84         for ( i = 0; entry->seq[ i ]; i++ )
  85         {
  86             bin     <<= 2;
  87             bin_rc1 >>= 2;
  88
  89             switch( entry->seq[ i ] )
  90             {
  91                 case 'A': case 'a':           bin_rc1 |= A_rc; j++; break;
  92                 case 'T': case 't': bin |= T;                  j++; break;
  93                 case 'C': case 'c': bin |= C; bin_rc1 |= G_rc; j++; break;
  94                 case 'G': case 'g': bin |= G; bin_rc1 |= C_rc; j++; break;
  95                 default: bin = 0; bin_rc1 = 0; j = 0; break;
  96             }
  97
  98             if ( j >= OLIGO_SIZE )
  99             {
 100                 array[ ( bin & mask ) ]++;
 101
 102                 bin_rc2 = bin_rc1;
 103
 104                 bin_rc2 >>= ( UINT_BITS - ( OLIGO_SIZE * 2 ) );
 105
 106                 array[ ( bin_rc2 ) ]++;
 107 /*
 108                 printf( "\n" );
 109                 printf( "mask          : %s\n", bits2string( mask ) );
 110                 printf( "bin           : %s\n", bits2string( bin ) );
 111                 printf( "bin     & mask: %s\n", bits2string( bin & mask ) );
 112                 printf( "bin_rc1       : %s\n", bits2string( bin_rc1 ) );
 113                 printf( "bin_rc2       : %s\n", bits2string( bin_rc2 ) );
 114 */
 115             }
 116         }
 117
 118         fprintf( stderr, "done.\n" );
 119     }
 120
 121     close_stream( fp );
 122
 123     free( entry->seq_name );
 124     free( entry->seq );
 125     entry = NULL;
 126
 127     return array;
 128 }
 129
 130
 131 uint mask_create( int oligo_size )
 132 {
 133     /* Martin A. Hansen, June 2008 */
 134
 135     /* Create a bit mask for binary encode oligos less than sizeof( uint ). */
 136
 137     uint i;
 138     uint mask;
 139
 140     mask = 0;
 141
 142     for ( i = 0; i < oligo_size; i++ )
 143     {
 144         mask <<= 2;
 145
 146         mask |= 3;
 147     }
 148
 149     return mask;
 150 }
 151
 152
 153 //void oligo_count_output( char *path, uint *array )
 154 //{
 155 //    /* Martin A. Hansen, June 2008 */
 156 //
 157 //    /* Output oligo count for each sequence position. */
 158 //
 159 //    struct seq_entry *entry;
 160 //    FILE             *fp;
 161 //    uint              mask;
 162 //    uint              i;
 163 //    uint              j;
 164 //    uint              bin;
 165 //    int               count;
 166 //    uint             *block;
 167 //    uint              block_pos;
 168 //    uint              block_beg;
 169 //    uint              block_size;
 170 //    uint              chr_pos;
 171 //    file_buffer      *buffer;
 172 //
 173 //    mask = mask_create( OLIGO_SIZE );
 174 //
 175 //    entry = mem_get( sizeof( entry ) );
 176 //
 177 //    fp = read_open( path );
 178 //
 179 //    while ( ( fasta_get_entry( buffer, &entry ) ) )
 180 //    {
 181 //        fprintf( stderr, "Writing results for: %s ... ", entry->seq_name );
 182 //
 183 //        bin        = 0;
 184 //        j          = 0;
 185 //        block_pos  = 0;
 186 //        block_size = sizeof( uint ) * ( entry->seq_len + OLIGO_SIZE );
 187 //        block      = mem_get_zero( block_size );
 188 //
 189 //        for ( i = 0; entry->seq[ i ]; i++ )
 190 //        {
 191 //            bin <<= 2;
 192 //
 193 //            switch( entry->seq[ i ] )
 194 //            {
 195 //                case 'A': case 'a':           j++; break;
 196 //                case 'T': case 't': bin |= T; j++; break;
 197 //                case 'C': case 'c': bin |= C; j++; break;
 198 //                case 'G': case 'g': bin |= G; j++; break;
 199 //                default: bin = 0; j = 0; break;
 200 //            }
 201 //
 202 //            if ( j >= OLIGO_SIZE )
 203 //            {
 204 //                count   = array[ ( bin & mask ) ];
 205 //
 206 //                if ( count > 1 )
 207 //                {
 208 //                    chr_pos = i - OLIGO_SIZE + 1;
 209 //
 210 //                    if ( block_pos == 0 )
 211 //                    {
 212 //                        memset( block, '\0', block_size );
 213 //
 214 //                        block_beg = chr_pos;
 215 //
 216 //                        block[ block_pos ] = count;
 217 //
 218 //                        block_pos++;
 219 //                    }
 220 //                    else
 221 //                    {
 222 //                        if ( chr_pos > block_beg + block_pos )
 223 //                        {
 224 //                            fixedstep_put_entry( entry->seq_name, block_beg, 1, block, block_pos );
 225 //
 226 //                            block_pos = 0;
 227 //                        }
 228 //                        else
 229 //                        {
 230 //                            block[ block_pos ] = count;
 231 //
 232 //                            block_pos++;
 233 //                        }
 234 //                    }
 235 //                }
 236 //            }
 237 //        }
 238 //
 239 //        if ( block_pos > 0 )
 240 //        {
 241 //            fixedstep_put_entry( entry->seq_name, block_beg, 1, block, block_pos );
 242 //
 243 //            mem_free( ( void * ) &block );
 244 //        }
 245 //
 246 //        fprintf( stderr, "done.\n" );
 247 //    }
 248 //
 249 //    close_stream( fp );
 250 //
 251 //    fasta_free_entry( entry );
 252 //}
 253
 254
 255 void fixedstep_put_entry( char *chr, int beg, int step_size, uint *block_array, int block_size )
 256 {
 257     /* Martin A. Hansen, June 2008 */
 258
 259     /* Outputs a block of fixedStep values. */
 260
 261     int i;
 262
 263     if ( block_size > 0 )
 264     {
 265         beg += 1; /* fixedStep format is 1 based. */
 266
 267         printf( "fixedStep chrom=%s start=%d step=%d\n", chr, beg, step_size );
 268
 269         for ( i = 0; i < block_size; i++ ) {
 270             printf( "%u\n", block_array[ i ] );
 271         }
 272     }
 273 }