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 "fasta.h"
  18
  19 // #define OLIGO_SIZE 15
  20 #define OLIGO_SIZE 5
  21 #define SIZE ( 1 << ( OLIGO_SIZE * 2 ) )
  22
  23 #define UINT_BITS 32
  24 #define T    3                            /* 11 on the rightmost two bits of bin. */
  25 #define C    1                            /* 01 on the rightmost two bits of bin. */
  26 #define G    2                            /* 10 on the rightmost two bits of bin. */
  27
  28 uint *oligo_count( char *path );
  29 uint  mask_create( int oligo_size );
  30 void  oligo_count_output( char *path, uint *array );
  31 void  fixedstep_put_entry( char *chr, int beg, int step_size, uint *block_array, int block_size );
  32
  33 int main( int argc, char *argv[] )
  34 {
  35     char *path  = NULL;
  36     uint *array = NULL;
  37
  38     path  = argv[ 1 ];
  39
  40     array = oligo_count( path );
  41
  42     //oligo_count_output( path, array );
  43
  44     return 0;
  45 }
  46
  47
  48 uint *oligo_count( char *path )
  49 {
  50     /* Martin A. Hansen, June 2008 */
  51
  52     /* Count the occurence of all oligos of a fixed size in a FASTA file. */
  53
  54     uint       *array   = NULL;
  55     uint        i       = 0;
  56     uint        mask    = 0;
  57     uint        bin     = 0;
  58     uint        bin_rc1 = 0;
  59     uint        bin_rc2 = 0;
  60     uint        j       = 0;
  61     uint        A_rc    = ( 3 << ( UINT_BITS - 2 ) );   /* 11 on the leftmost two bits an uint. */
  62     uint        G_rc    = ( 2 << ( UINT_BITS - 2 ) );   /* 10 on the leftmost two bits an uint. */
  63     uint        C_rc    = ( 1 << ( UINT_BITS - 2 ) );   /* 01 on the leftmost two bits an uint. */
  64     seq_entry  *entry   = NULL;
  65     FILE       *fp      = NULL;
  66     file_buffer *buffer = 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( &buffer, &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     fasta_free_entry( entry );
 124
 125     return array;
 126 }
 127
 128
 129 uint mask_create( int oligo_size )
 130 {
 131     /* Martin A. Hansen, June 2008 */
 132
 133     /* Create a bit mask for binary encode oligos less than sizeof( uint ). */
 134
 135     uint i;
 136     uint mask;
 137
 138     mask = 0;
 139
 140     for ( i = 0; i < oligo_size; i++ )
 141     {
 142         mask <<= 2;
 143
 144         mask |= 3;
 145     }
 146
 147     return mask;
 148 }
 149
 150
 151 //void oligo_count_output( char *path, uint *array )
 152 //{
 153 //    /* Martin A. Hansen, June 2008 */
 154 //
 155 //    /* Output oligo count for each sequence position. */
 156 //
 157 //    struct seq_entry *entry;
 158 //    FILE             *fp;
 159 //    uint              mask;
 160 //    uint              i;
 161 //    uint              j;
 162 //    uint              bin;
 163 //    int               count;
 164 //    uint             *block;
 165 //    uint              block_pos;
 166 //    uint              block_beg;
 167 //    uint              block_size;
 168 //    uint              chr_pos;
 169 //    file_buffer      *buffer;
 170 //
 171 //    mask = mask_create( OLIGO_SIZE );
 172 //
 173 //    entry = mem_get( sizeof( entry ) );
 174 //
 175 //    fp = read_open( path );
 176 //
 177 //    while ( ( fasta_get_entry( buffer, &entry ) ) )
 178 //    {
 179 //        fprintf( stderr, "Writing results for: %s ... ", entry->seq_name );
 180 //
 181 //        bin        = 0;
 182 //        j          = 0;
 183 //        block_pos  = 0;
 184 //        block_size = sizeof( uint ) * ( entry->seq_len + OLIGO_SIZE );
 185 //        block      = mem_get_zero( block_size );
 186 //
 187 //        for ( i = 0; entry->seq[ i ]; i++ )
 188 //        {
 189 //            bin <<= 2;
 190 //
 191 //            switch( entry->seq[ i ] )
 192 //            {
 193 //                case 'A': case 'a':           j++; break;
 194 //                case 'T': case 't': bin |= T; j++; break;
 195 //                case 'C': case 'c': bin |= C; j++; break;
 196 //                case 'G': case 'g': bin |= G; j++; break;
 197 //                default: bin = 0; j = 0; break;
 198 //            }
 199 //
 200 //            if ( j >= OLIGO_SIZE )
 201 //            {
 202 //                count   = array[ ( bin & mask ) ];
 203 //
 204 //                if ( count > 1 )
 205 //                {
 206 //                    chr_pos = i - OLIGO_SIZE + 1;
 207 //
 208 //                    if ( block_pos == 0 )
 209 //                    {
 210 //                        memset( block, '\0', block_size );
 211 //
 212 //                        block_beg = chr_pos;
 213 //
 214 //                        block[ block_pos ] = count;
 215 //
 216 //                        block_pos++;
 217 //                    }
 218 //                    else
 219 //                    {
 220 //                        if ( chr_pos > block_beg + block_pos )
 221 //                        {
 222 //                            fixedstep_put_entry( entry->seq_name, block_beg, 1, block, block_pos );
 223 //
 224 //                            block_pos = 0;
 225 //                        }
 226 //                        else
 227 //                        {
 228 //                            block[ block_pos ] = count;
 229 //
 230 //                            block_pos++;
 231 //                        }
 232 //                    }
 233 //                }
 234 //            }
 235 //        }
 236 //
 237 //        if ( block_pos > 0 )
 238 //        {
 239 //            fixedstep_put_entry( entry->seq_name, block_beg, 1, block, block_pos );
 240 //
 241 //            mem_free( ( void * ) &block );
 242 //        }
 243 //
 244 //        fprintf( stderr, "done.\n" );
 245 //    }
 246 //
 247 //    close_stream( fp );
 248 //
 249 //    fasta_free_entry( entry );
 250 //}
 251
 252
 253 void fixedstep_put_entry( char *chr, int beg, int step_size, uint *block_array, int block_size )
 254 {
 255     /* Martin A. Hansen, June 2008 */
 256
 257     /* Outputs a block of fixedStep values. */
 258
 259     int i;
 260
 261     if ( block_size > 0 )
 262     {
 263         beg += 1; /* fixedStep format is 1 based. */
 264
 265         printf( "fixedStep chrom=%s start=%d step=%d\n", chr, beg, step_size );
 266
 267         for ( i = 0; i < block_size; i++ ) {
 268             printf( "%u\n", block_array[ i ] );
 269         }
 270     }
 271 }