BGI/SOAPdenovo2/standardPregraph/kmer.c

   1 /*
   2  * kmer.c
   3  *
   4  * Copyright (c) 2008-2012 BGI-Shenzhen <soap at genomics dot org dot cn>.
   5  *
   6  * This file is part of SOAPdenovo.
   7  *
   8  * SOAPdenovo is free software: you can redistribute it and/or modify
   9  * it under the terms of the GNU General Public License as published by
  10  * the Free Software Foundation, either version 3 of the License, or
  11  * (at your option) any later version.
  12  *
  13  * SOAPdenovo is distributed in the hope that it will be useful,
  14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  * GNU General Public License for more details.
  17  *
  18  * You should have received a copy of the GNU General Public License
  19  * along with SOAPdenovo.  If not, see <http://www.gnu.org/licenses/>.
  20  *
  21  */
  22
  23 #include "stdinc.h"
  24 #include "newhash.h"
  25 #include "kmerhash.h"
  26 #include "extfunc.h"
  27 #include "extvab.h"
  28
  29 #ifdef MER127
  30 void PrintKmer ( FILE * fp, Kmer kmer )
  31 {
  32         fprintf ( fp, "%llx %llx %llx %llx", kmer.high1, kmer.low1, kmer.high2, kmer.low2 );
  33 }
  34
  35 /*************************************************
  36 Function:
  37     KmerSmaller
  38 Description:
  39     Checks if kmer1 is smaller than kmer2. if no, returns 0.
  40 Input:
  41     1. kmer1:       the first kmer
  42     2. kmer2:       the second kmer
  43 Output:
  44     None.
  45 Return:
  46     1 if kmer1 is smaller than kmer2.
  47 *************************************************/
  48 boolean KmerSmaller ( Kmer kmer1, Kmer kmer2 )
  49 {
  50         if ( kmer1.high1 != kmer2.high1 )
  51         {
  52                 return ( kmer1.high1 < kmer2.high1 );
  53         }
  54         else
  55         {
  56                 if ( kmer1.low1 != kmer2.low1 )
  57                 {
  58                         return ( kmer1.low1 < kmer2.low1 );
  59                 }
  60                 else
  61                 {
  62                         if ( kmer1.high2 != kmer2.high2 )
  63                         {
  64                                 return ( kmer1.high2 < kmer2.high2 );
  65                         }
  66                         else
  67                         {
  68                                 return ( kmer1.low2 < kmer2.low2 );
  69                         }
  70                 }
  71         }
  72 }
  73
  74 /*************************************************
  75 Function:
  76     KmerLarger
  77 Description:
  78     Checks if kmer1 is larger than kmer2. if no, returns 0.
  79 Input:
  80     1. kmer1:       the first kmer
  81     2. kmer2:       the second kmer
  82 Output:
  83     None.
  84 Return:
  85     1 if kmer1 is larger than kmer2.
  86 *************************************************/
  87 boolean KmerLarger ( Kmer kmer1, Kmer kmer2 )
  88 {
  89         if ( kmer1.high1 != kmer2.high1 )
  90         {
  91                 return ( kmer1.high1 > kmer2.high1 );
  92         }
  93         else
  94         {
  95                 if ( kmer1.low1 != kmer2.low1 )
  96                 {
  97                         return ( kmer1.low1 > kmer2.low1 );
  98                 }
  99                 else
 100                 {
 101                         if ( kmer1.high2 != kmer2.high2 )
 102                         {
 103                                 return ( kmer1.high2 > kmer2.high2 );
 104                         }
 105                         else
 106                         {
 107                                 return ( kmer1.low2 > kmer2.low2 );
 108                         }
 109                 }
 110         }
 111 }
 112
 113 /*************************************************
 114 Function:
 115     KmerEqual
 116 Description:
 117     Checks if kmer1 is same as kmer2. if no, returns 0.
 118 Input:
 119     1. kmer1:       the first kmer
 120     2. kmer2:       the second kmer
 121 Output:
 122     None.
 123 Return:
 124     1 if if kmer1 is same as kmer2.
 125 *************************************************/
 126 boolean KmerEqual ( Kmer kmer1, Kmer kmer2 )
 127 {
 128         if ( kmer1.low2 != kmer2.low2 || kmer1.high2 != kmer2.high2 || kmer1.low1 != kmer2.low1 || kmer1.high1 != kmer2.high1 )
 129         {
 130                 return 0;
 131         }
 132         else
 133         {
 134                 return 1;
 135         }
 136 }
 137
 138 /*************************************************
 139 Function:
 140     KmerAnd
 141 Description:
 142     kmer1 & kmer2
 143 Input:
 144     1. kmer1:       the first kmer
 145     2. kmer2:       the second kmer
 146 Output:
 147     None.
 148 Return:
 149     kmer1 & kmer2.
 150 *************************************************/
 151 Kmer KmerAnd ( Kmer kmer1, Kmer kmer2 )
 152 {
 153         kmer1.high1 &= kmer2.high1;
 154         kmer1.low1 &= kmer2.low1;
 155         kmer1.high2 &= kmer2.high2;
 156         kmer1.low2 &= kmer2.low2;
 157         return kmer1;
 158 }
 159
 160 /*************************************************
 161 Function:
 162     KmerLeftBitMoveBy2
 163 Description:
 164     word <<= 2.
 165 Input:
 166     1. word:        the kmer
 167 Output:
 168     None.
 169 Return:
 170     word <<= 2.
 171 *************************************************/
 172 Kmer KmerLeftBitMoveBy2 ( Kmer word )
 173 {
 174         word.high1 = ( word.high1 << 2 ) | ( word.low1 >> 62 );
 175         word.low1 = ( word.low1 << 2 ) | ( word.high2 >> 62 );
 176         word.high2 = ( word.high2 << 2 ) | ( word.low2 >> 62 );
 177         word.low2 <<= 2;
 178         return word;
 179 }
 180
 181 /*************************************************
 182 Function:
 183     KmerRightBitMoveBy2
 184 Description:
 185     word >>= 2.
 186 Input:
 187     1. word:        the kmer
 188 Output:
 189     None.
 190 Return:
 191     word >>= 2.
 192 *************************************************/
 193 Kmer KmerRightBitMoveBy2 ( Kmer word )
 194 {
 195         word.low2 = ( word.low2 >> 2 ) | ( word.high2 & 0x3 ) << 62;
 196         word.high2 = ( word.high2 >> 2 ) | ( word.low1 & 0x3 ) << 62;
 197         word.low1 = ( word.low1 >> 2 ) | ( word.high1 & 0x3 ) << 62;
 198         word.high1 >>= 2;
 199         return word;
 200 }
 201
 202 /*************************************************
 203 Function:
 204     KmerPlus
 205 Description:
 206     Appends a base to kmer to get (k+1)mer.
 207 Input:
 208     1. prev:        the kmer
 209     2. ch:      the downstream
 210 Output:
 211     None.
 212 Return:
 213     (k+1)mer.
 214 *************************************************/
 215 Kmer KmerPlus ( Kmer prev, char ch )
 216 {
 217         Kmer word = KmerLeftBitMoveBy2 ( prev );
 218         word.low2 |= ch;
 219         return word;
 220 }
 221
 222 /*************************************************
 223 Function:
 224     nextKmer
 225 Description:
 226     Combines the last (k-1) bases of a kmer with another base to get new kmer.
 227 Input:
 228     1. prev:        the kmer
 229     2. ch:      the appened base
 230 Output:
 231     None.
 232 Return:
 233     The new kmer.
 234 *************************************************/
 235 Kmer nextKmer ( Kmer prev, char ch )
 236 {
 237         Kmer word = KmerLeftBitMoveBy2 ( prev );
 238         word = KmerAnd ( word, WORDFILTER );
 239         word.low2 |= ch;
 240         return word;
 241 }
 242
 243 /*************************************************
 244 Function:
 245     prevKmer
 246 Description:
 247     Combines a base with the first (k-1) bases of a kmer to get new kmer.
 248 Input:
 249     1. next:        the kmer
 250     2. ch:      the first base of new kmer
 251 Output:
 252     None.
 253 Return:
 254     The new kmer.
 255 *************************************************/
 256 Kmer prevKmer ( Kmer next, char ch )
 257 {
 258         Kmer word = KmerRightBitMoveBy2 ( next );
 259
 260         switch ( overlaplen )
 261         {
 262                 case 1 ... 32:
 263                         word.low2 |= ( ( ( ubyte8 ) ch ) << 2 * ( overlaplen - 1 ) );
 264                         break;
 265                 case 33 ... 64:
 266                         word.high2 |= ( ( ubyte8 ) ch ) << ( 2 * ( overlaplen - 1 ) - 64 );
 267                         break;
 268                 case 65 ... 96:
 269                         word.low1 |= ( ( ubyte8 ) ch ) << ( 2 * ( overlaplen - 1 ) - 128 );
 270                         break;
 271                 case 97 ... 128:
 272                         word.high1 |= ( ( ubyte8 ) ch ) << ( 2 * ( overlaplen - 1 ) - 192 );
 273                         break;
 274         }
 275
 276         return word;
 277 }
 278
 279 /*************************************************
 280 Function:
 281     lastCharInKmer
 282 Description:
 283     Gets the last char of a kmer, and returns it.
 284 Input:
 285     1. kmer:        the kmer
 286 Output:
 287     None.
 288 Return:
 289     The last char of kmer.
 290 *************************************************/
 291 char lastCharInKmer ( Kmer kmer )
 292 {
 293         return ( char ) ( kmer.low2 & 0x3 );
 294 }
 295
 296 /*************************************************
 297 Function:
 298     firstCharInKmer
 299 Description:
 300     Gets the first char of a kmer, and returns it.
 301 Input:
 302     1. kmer:        the kmer
 303 Output:
 304     None.
 305 Return:
 306     The first char of kmer.
 307 *************************************************/
 308 char firstCharInKmer ( Kmer kmer )
 309 {
 310         switch ( overlaplen )
 311         {
 312                 case 1 ... 32:
 313                         kmer.low2 >>= 2 * ( overlaplen - 1 );
 314                         return kmer.low2;   // & 3;
 315                 case 33 ... 64:
 316                         kmer.high2 >>= 2 * ( overlaplen - 1 ) - 64;
 317                         return kmer.high2;  // & 3;
 318                 case 65 ... 96:
 319                         kmer.low1 >>= 2 * ( overlaplen - 1 ) - 128;
 320                         return kmer.low1;
 321                 case 97 ... 128:
 322                         kmer.high1 >>= 2 * ( overlaplen - 1 ) - 192;
 323                         return kmer.high1;
 324         }
 325 }
 326
 327 /*************************************************
 328 Function:
 329     createFilter
 330 Description:
 331     Creates filter for masking.
 332 Input:
 333     1. overlaplen:      the kmer value
 334 Output:
 335     None.
 336 Return:
 337     The filter kmer.
 338 *************************************************/
 339 Kmer createFilter ( int overlaplen )
 340 {
 341         Kmer word;
 342         word.high1 = word.low1 = word.high2 = word.low2 = 0;
 343
 344         switch ( overlaplen )
 345         {
 346                 case 1 ... 31:
 347                         word.low2 = ( ( ( ubyte8 ) 1 ) << ( 2 * overlaplen ) ) - 1;
 348                         break;
 349                 case 32 ... 63:
 350                         word.low2 = ~word.low2;
 351                         word.high2 = ( ( ( ubyte8 ) 1 ) << ( 2 * overlaplen - 64 ) ) - 1;
 352                         break;
 353                 case 64 ... 95:
 354                         word.high2 = word.low2 = ~word.low2;
 355                         word.low1 = ( ( ( ubyte8 ) 1 ) << ( 2 * overlaplen - 128 ) ) - 1;
 356                         break;
 357                 case 96 ... 127:
 358                         word.low1 = word.high2 = word.low2 = ~word.low2;
 359                         word.high1 = ( ( ( ubyte8 ) 1 ) << ( 2 * overlaplen - 192 ) ) - 1;
 360                         break;
 361         }
 362
 363         return word;
 364 }
 365
 366 /*************************************************
 367 Function:
 368     KmerRightBitMove
 369 Description:
 370     Makes right shifts to a kmer by specified bits and returns the new kmer.
 371 Input:
 372     1. word:        the kmer
 373     2. dis:     the shifted bits
 374 Output:
 375     None.
 376 Return:
 377     The new kmer.
 378 *************************************************/
 379 Kmer KmerRightBitMove ( Kmer word, int dis )
 380 {
 381         ubyte8 mask;
 382
 383         switch ( dis )
 384         {
 385                 case 0 ... 63:
 386                         mask = ( ( ( ubyte8 ) 1 ) << dis ) - 1;
 387                         word.low2 = ( word.low2 >> dis ) | ( word.high2 & mask ) << ( 64 - dis );
 388                         word.high2 = ( word.high2 >> dis ) | ( word.low1 & mask ) << ( 64 - dis );
 389                         word.low1 = ( word.low1 >> dis ) | ( word.high1 & mask ) << ( 64 - dis );
 390                         word.high1 >>= dis;
 391                         return word;
 392                 case 64 ... 127:
 393                         mask = ( ( ( ubyte8 ) 1 ) << ( dis - 64 ) ) - 1;
 394                         word.low2 = word.high2 >> ( dis - 64 ) | ( word.low1 & mask ) << ( 128 - dis );
 395                         word.high2 = word.low1 >> ( dis - 64 ) | ( word.high1 & mask ) << ( 128 - dis );
 396                         word.low1 = word.high1 >> ( dis - 64 );
 397                         word.high1 = 0;
 398                         return word;
 399                 case 128 ... 191:
 400                         mask = ( ( ( ubyte8 ) 1 ) << ( dis - 128 ) ) - 1;
 401                         word.low2 = word.low1 >> ( dis - 128 ) | ( word.high1 & mask ) << ( 192 - dis );
 402                         word.high2 = word.high1 >> ( dis - 128 );
 403                         word.high1 = word.low1 = 0;
 404                         return word;
 405                 case 192 ... 255:
 406                         word.low2 = word.high1 >> ( dis - 192 );
 407                         word.high1 = word.low1 = word.high2 = 0;
 408                         return word;
 409         }
 410 }
 411
 412 /*************************************************
 413 Function:
 414     printKmerSeq
 415 Description:
 416     Prints seq of kmer.
 417 Input:
 418     1. fp:      output file
 419     2. kmer:        the output kmer
 420 Output:
 421     None.
 422 Return:
 423     None.
 424 *************************************************/
 425 void printKmerSeq ( FILE * fp, Kmer kmer )
 426 {
 427         int i, bit1, bit2, bit3, bit4;
 428         bit4 = bit3 = bit2 = bit1 = 0;
 429         char kmerSeq[128];
 430
 431         switch ( overlaplen )
 432         {
 433                 case 1 ... 31:
 434                         bit4 = overlaplen;
 435                         break;
 436                 case 32 ... 63:
 437                         bit4 = 32;
 438                         bit3 = overlaplen - 32;
 439                         break;
 440                 case 64 ... 95:
 441                         bit4 = bit3 = 32;
 442                         bit2 = overlaplen - 64;
 443                         break;
 444                 case 96 ... 127:
 445                         bit4 = bit3 = bit2 = 32;
 446                         bit1 = overlaplen - 96;
 447                         break;
 448         }
 449
 450         for ( i = bit1 - 1; i >= 0; i-- )
 451         {
 452                 kmerSeq[i] = kmer.high1 & 0x3;
 453                 kmer.high1 >>= 2;
 454         }
 455
 456         for ( i = bit2 - 1; i >= 0; i-- )
 457         {
 458                 kmerSeq[i + bit1] = kmer.low1 & 0x3;
 459                 kmer.low1 >>= 2;
 460         }
 461
 462         for ( i = bit3 - 1; i >= 0; i-- )
 463         {
 464                 kmerSeq[i + bit1 + bit2] = kmer.high2 & 0x3;
 465                 kmer.high2 >>= 2;
 466         }
 467
 468         for ( i = bit4 - 1; i >= 0; i-- )
 469         {
 470                 kmerSeq[i + bit1 + bit2 + bit3] = kmer.low2 & 0x3;
 471                 kmer.low2 >>= 2;
 472         }
 473
 474         for ( i = 0; i < overlaplen; i++ )
 475         {
 476                 fprintf ( fp, "%c", int2base ( ( int ) kmerSeq[i] ) );
 477         }
 478 }
 479
 480 void print_kmer ( FILE * fp, Kmer kmer, char c )
 481 {
 482         fprintf ( fp, "%llx %llx %llx %llx", kmer.high1, kmer.low1, kmer.high2, kmer.low2 );
 483         fprintf ( fp, "%c", c );
 484 }
 485
 486 void print_kmer_gz ( gzFile * fp, Kmer kmer, char c )
 487 {
 488         gzprintf ( fp, "%llx %llx %llx %llx", kmer.high1, kmer.low1, kmer.high2, kmer.low2 );
 489         gzprintf ( fp, "%c", c );
 490 }
 491
 492 static const ubyte2 BitReverseTable[65536] =
 493 {
 494 # define R2(n)      n,     n + 1*16384,    n + 2*16384,    n + 3*16384
 495 # define R4(n)   R2(n), R2(n + 1*4096), R2(n + 2*4096), R2(n + 3*4096)
 496 # define R6(n)   R4(n), R4(n + 1*1024), R4(n + 2*1024), R4(n + 3*1024)
 497 # define R8(n)   R6(n), R6(n + 1*256 ), R6(n + 2*256 ), R6(n + 3*256 )
 498 # define R10(n)  R8(n), R8(n + 1*64  ), R8(n + 2*64  ), R8(n + 3*64  )
 499 # define R12(n)  R10(n),R10(n + 1*16), R10(n + 2*16 ), R10(n + 3*16  )
 500 # define R14(n)  R12(n),R12(n + 1*4 ), R12(n + 2*4  ), R12(n + 3*4  )
 501         R14 ( 0 ), R14 ( 1 ), R14 ( 2 ), R14 ( 3 )
 502 };
 503
 504 /*************************************************
 505 Function:
 506     fastReverseComp
 507 Description:
 508     Gets the reverse complement of a kmer.
 509 Input:
 510     1. seq:         the kmer
 511     2. seq_size:        the length of the kmer
 512 Output:
 513     None.
 514 Return:
 515     The reverse complement of kmer.
 516 *************************************************/
 517 static Kmer fastReverseComp ( Kmer seq, char seq_size )
 518 {
 519         seq.low2 ^= 0xAAAAAAAAAAAAAAAALLU;
 520         seq.low2 = ( ( ubyte8 ) BitReverseTable[seq.low2 & 0xffff] << 48 ) |
 521                    ( ( ubyte8 ) BitReverseTable[ ( seq.low2 >> 16 ) & 0xffff] << 32 ) | ( ( ubyte8 ) BitReverseTable[ ( seq.low2 >> 32 ) & 0xffff] << 16 ) | ( ( ubyte8 ) BitReverseTable[ ( seq.low2 >> 48 ) & 0xffff] );
 522
 523         if ( seq_size < 32 )
 524         {
 525                 seq.low2 >>= ( 64 - ( seq_size << 1 ) );
 526                 return seq;
 527         }
 528
 529         seq.high2 ^= 0xAAAAAAAAAAAAAAAALLU;
 530         seq.high2 = ( ( ubyte8 ) BitReverseTable[seq.high2 & 0xffff] << 48 ) |
 531                     ( ( ubyte8 ) BitReverseTable[ ( seq.high2 >> 16 ) & 0xffff] << 32 ) | ( ( ubyte8 ) BitReverseTable[ ( seq.high2 >> 32 ) & 0xffff] << 16 ) | ( ( ubyte8 ) BitReverseTable[ ( seq.high2 >> 48 ) & 0xffff] );
 532
 533         if ( seq_size < 64 )
 534         {
 535                 seq.high2 = seq.high2 ^ seq.low2;
 536                 seq.low2 = seq.high2 ^ seq.low2;
 537                 seq.high2 = seq.high2 ^ seq.low2;
 538                 seq = KmerRightBitMove ( seq, 128 - ( seq_size << 1 ) );
 539                 return seq;
 540         }
 541
 542         seq.low1 ^= 0xAAAAAAAAAAAAAAAALLU;
 543         seq.low1 = ( ( ubyte8 ) BitReverseTable[seq.low1 & 0xffff] << 48 ) |
 544                    ( ( ubyte8 ) BitReverseTable[ ( seq.low1 >> 16 ) & 0xffff] << 32 ) | ( ( ubyte8 ) BitReverseTable[ ( seq.low1 >> 32 ) & 0xffff] << 16 ) | ( ( ubyte8 ) BitReverseTable[ ( seq.low1 >> 48 ) & 0xffff] );
 545
 546         if ( seq_size < 96 )
 547         {
 548                 seq.low1 = seq.low1 ^ seq.low2;
 549                 seq.low2 = seq.low1 ^ seq.low2;
 550                 seq.low1 = seq.low1 ^ seq.low2;
 551                 seq = KmerRightBitMove ( seq, 192 - ( seq_size << 1 ) );
 552                 return seq;
 553         }
 554
 555         seq.high1 ^= 0xAAAAAAAAAAAAAAAALLU;
 556         seq.high1 = ( ( ubyte8 ) BitReverseTable[seq.high1 & 0xffff] << 48 ) |
 557                     ( ( ubyte8 ) BitReverseTable[ ( seq.high1 >> 16 ) & 0xffff] << 32 ) | ( ( ubyte8 ) BitReverseTable[ ( seq.high1 >> 32 ) & 0xffff] << 16 ) | ( ( ubyte8 ) BitReverseTable[ ( seq.high1 >> 48 ) & 0xffff] );
 558         seq.low1 = seq.low1 ^ seq.high2;
 559         seq.high2 = seq.low1 ^ seq.high2;
 560         seq.low1 = seq.low1 ^ seq.high2;
 561         seq.low2 = seq.low2 ^ seq.high1;
 562         seq.high1 = seq.low2 ^ seq.high1;
 563         seq.low2 = seq.low2 ^ seq.high1;
 564         seq = KmerRightBitMove ( seq, 256 - ( seq_size << 1 ) );
 565         return seq;
 566 }
 567
 568 Kmer reverseComplementVerbose ( Kmer word, int overlap )
 569 {
 570         return fastReverseComp ( word, overlap );
 571 }
 572
 573 Kmer reverseComplement ( Kmer word, int overlap )
 574 {
 575         return fastReverseComp ( word, overlap );
 576 }
 577
 578 #else
 579 void PrintKmer ( FILE * fp, Kmer kmer )
 580 {
 581         fprintf ( fp, "%llx %llx", kmer.high, kmer.low );
 582 }
 583
 584
 585 __uint128_t Kmer2int128 ( Kmer seq )
 586 {
 587         __uint128_t temp;
 588         temp = seq.high;
 589         temp <<= 64;
 590         temp |= seq.low;
 591         return temp;
 592 }
 593 boolean KmerSmaller ( Kmer kmer1, Kmer kmer2 )
 594 {
 595         if ( kmer1.high < kmer2.high )
 596         {
 597                 return 1;
 598         }
 599         else if ( kmer1.high == kmer2.high )
 600         {
 601                 if ( kmer1.low < kmer2.low )
 602                 {
 603                         return 1;
 604                 }
 605                 else
 606                 {
 607                         return 0;
 608                 }
 609         }
 610         else
 611         {
 612                 return 0;
 613         }
 614 }
 615
 616 boolean KmerLarger ( Kmer kmer1, Kmer kmer2 )
 617 {
 618         if ( kmer1.high > kmer2.high )
 619         {
 620                 return 1;
 621         }
 622         else if ( kmer1.high == kmer2.high )
 623         {
 624                 if ( kmer1.low > kmer2.low )
 625                 {
 626                         return 1;
 627                 }
 628                 else
 629                 {
 630                         return 0;
 631                 }
 632         }
 633         else
 634         {
 635                 return 0;
 636         }
 637 }
 638 boolean KmerEqual ( Kmer kmer1, Kmer kmer2 )
 639 {
 640         if ( kmer1.high == kmer2.high && kmer1.low == kmer2.low )
 641         {
 642                 return 1;
 643         }
 644         else
 645         {
 646                 return 0;
 647         }
 648 }
 649
 650 Kmer KmerAnd ( Kmer kmer1, Kmer kmer2 )
 651 {
 652         kmer1.high &= kmer2.high;
 653         kmer1.low &= kmer2.low;
 654         return kmer1;
 655 }
 656
 657 Kmer KmerLeftBitMoveBy2 ( Kmer word )
 658 {
 659         ubyte8 temp = word.low >> 62;
 660         word.high <<= 2;
 661         word.high |= temp;
 662         word.low <<= 2;
 663         return word;
 664 }
 665
 666 Kmer KmerRightBitMoveBy2 ( Kmer word )
 667 {
 668         ubyte8 temp = ( word.high & 0x3 ) << 62;
 669         word.high >>= 2;
 670         word.low >>= 2;
 671         word.low |= temp;
 672         return word;
 673 }
 674
 675 Kmer KmerPlus ( Kmer prev, char ch )
 676 {
 677         Kmer word = KmerLeftBitMoveBy2 ( prev );
 678         word.low |= ch;
 679         return word;
 680 }
 681 Kmer nextKmer ( Kmer prev, char ch )
 682 {
 683         Kmer word = KmerLeftBitMoveBy2 ( prev );
 684         word = KmerAnd ( word, WORDFILTER );
 685         word.low |= ch;
 686         return word;
 687 }
 688
 689 Kmer prevKmer ( Kmer next, char ch )
 690 {
 691         Kmer word = KmerRightBitMoveBy2 ( next );
 692
 693         if ( 2 * ( overlaplen - 1 ) < 64 )
 694         {
 695                 word.low |= ( ( ( ubyte8 ) ch ) << 2 * ( overlaplen - 1 ) );
 696         }
 697         else
 698         {
 699                 word.high |= ( ( ubyte8 ) ch ) << ( 2 * ( overlaplen - 1 ) - 64 );
 700         }
 701
 702         return word;
 703 }
 704
 705 char lastCharInKmer ( Kmer kmer )
 706 {
 707         return ( char ) ( kmer.low & 0x3 );
 708 }
 709 char firstCharInKmer ( Kmer kmer )
 710 {
 711         if ( 2 * ( overlaplen - 1 ) < 64 )
 712         {
 713                 kmer.low >>= 2 * ( overlaplen - 1 );
 714                 return kmer.low;    // & 3;
 715         }
 716         else
 717         {
 718                 kmer.high >>= 2 * ( overlaplen - 1 ) - 64;
 719                 return kmer.high;   // & 3;
 720         }
 721 }
 722
 723 Kmer createFilter ( int overlaplen )
 724 {
 725         Kmer word;
 726         word.high = word.low = 0;
 727
 728         if ( 2 * overlaplen < 64 )
 729         {
 730                 word.low = ( ( ( ubyte8 ) 1 ) << ( 2 * overlaplen ) ) - 1;
 731         }
 732         else
 733         {
 734                 word.low = ~word.low;
 735
 736                 if ( 2 * overlaplen > 64 )
 737                 {
 738                         word.high = ( ( ( ubyte8 ) 1 ) << ( 2 * overlaplen - 64 ) ) - 1;
 739                 }
 740         }
 741
 742         return word;
 743 }
 744
 745 Kmer KmerRightBitMove ( Kmer word, int dis )
 746 {
 747         if ( dis < 64 )
 748         {
 749                 ubyte8 mask = ( ( ( ubyte8 ) 1 ) << dis ) - 1;
 750                 ubyte8 temp = ( word.high & mask ) << ( 64 - dis );
 751                 word.high >>= dis;
 752                 word.low >>= dis;
 753                 word.low |= temp;
 754                 return word;
 755         }
 756
 757         word.high >>= ( dis - 64 );
 758         word.low = word.high;
 759         word.high = 0;
 760         return word;
 761 }
 762
 763
 764 void printKmerSeq ( FILE * fp, Kmer kmer )
 765 {
 766         int i, bit1, bit2;
 767         char ch;
 768         char kmerSeq[64];
 769         bit2 = overlaplen > 32 ? 32 : overlaplen;
 770         bit1 = overlaplen > 32 ? overlaplen - 32 : 0;
 771
 772         for ( i = bit1 - 1; i >= 0; i-- )
 773         {
 774                 ch = kmer.high & 0x3;
 775                 kmer.high >>= 2;
 776                 kmerSeq[i] = ch;
 777         }
 778
 779         for ( i = bit2 - 1; i >= 0; i-- )
 780         {
 781                 ch = kmer.low & 0x3;
 782                 kmer.low >>= 2;
 783                 kmerSeq[i + bit1] = ch;
 784         }
 785
 786         for ( i = 0; i < overlaplen; i++ )
 787         {
 788                 fprintf ( fp, "%c", int2base ( ( int ) kmerSeq[i] ) );
 789         }
 790 }
 791
 792 void print_kmer ( FILE * fp, Kmer kmer, char c )
 793 {
 794         fprintf ( fp, "%llx %llx", kmer.high, kmer.low );
 795         fprintf ( fp, "%c", c );
 796 }
 797
 798 void print_kmer_gz ( gzFile * fp, Kmer kmer, char c )
 799 {
 800         gzprintf ( fp, "%llx %llx", kmer.high, kmer.low );
 801         gzprintf ( fp, "%c", c );
 802 }
 803
 804 static Kmer fastReverseComp ( Kmer seq, char seq_size )
 805 {
 806         seq.low ^= 0xAAAAAAAAAAAAAAAALLU;
 807         seq.low = ( ( seq.low & 0x3333333333333333LLU ) << 2 ) | ( ( seq.low & 0xCCCCCCCCCCCCCCCCLLU ) >> 2 );
 808         seq.low = ( ( seq.low & 0x0F0F0F0F0F0F0F0FLLU ) << 4 ) | ( ( seq.low & 0xF0F0F0F0F0F0F0F0LLU ) >> 4 );
 809         seq.low = ( ( seq.low & 0x00FF00FF00FF00FFLLU ) << 8 ) | ( ( seq.low & 0xFF00FF00FF00FF00LLU ) >> 8 );
 810         seq.low = ( ( seq.low & 0x0000FFFF0000FFFFLLU ) << 16 ) | ( ( seq.low & 0xFFFF0000FFFF0000LLU ) >> 16 );
 811         seq.low = ( ( seq.low & 0x00000000FFFFFFFFLLU ) << 32 ) | ( ( seq.low & 0xFFFFFFFF00000000LLU ) >> 32 );
 812
 813         if ( seq_size < 32 )
 814         {
 815                 seq.low >>= ( 64 - ( seq_size << 1 ) );
 816                 return seq;
 817         }
 818
 819         seq.high ^= 0xAAAAAAAAAAAAAAAALLU;
 820         seq.high = ( ( seq.high & 0x3333333333333333LLU ) << 2 ) | ( ( seq.high & 0xCCCCCCCCCCCCCCCCLLU ) >> 2 );
 821         seq.high = ( ( seq.high & 0x0F0F0F0F0F0F0F0FLLU ) << 4 ) | ( ( seq.high & 0xF0F0F0F0F0F0F0F0LLU ) >> 4 );
 822         seq.high = ( ( seq.high & 0x00FF00FF00FF00FFLLU ) << 8 ) | ( ( seq.high & 0xFF00FF00FF00FF00LLU ) >> 8 );
 823         seq.high = ( ( seq.high & 0x0000FFFF0000FFFFLLU ) << 16 ) | ( ( seq.high & 0xFFFF0000FFFF0000LLU ) >> 16 );
 824         seq.high = ( ( seq.high & 0x00000000FFFFFFFFLLU ) << 32 ) | ( ( seq.high & 0xFFFFFFFF00000000LLU ) >> 32 );
 825         ubyte8 temp = seq.high;
 826         seq.high = seq.low;
 827         seq.low = temp;
 828         seq = KmerRightBitMove ( seq, 128 - ( seq_size << 1 ) );
 829         return seq;
 830 }
 831
 832 Kmer reverseComplementVerbose ( Kmer word, int overlap )
 833 {
 834         return fastReverseComp ( word, overlap );
 835 }
 836
 837 Kmer reverseComplement ( Kmer word, int overlap )
 838 {
 839         return fastReverseComp ( word, overlap );
 840 }
 841
 842 #endif