misc/seqtk.c

   1 #include <stdio.h>
   2 #include <ctype.h>
   3 #include <stdlib.h>
   4 #include <stdint.h>
   5 #include <zlib.h>
   6 #include <string.h>
   7 #include <unistd.h>
   8 #include <limits.h>
   9
  10 #include "kseq.h"
  11 KSEQ_INIT(gzFile, gzread)
  12
  13 typedef struct {
  14         int n, m;
  15         uint64_t *a;
  16 } reglist_t;
  17
  18 #include "khash.h"
  19 KHASH_MAP_INIT_STR(reg, reglist_t)
  20
  21 typedef kh_reg_t reghash_t;
  22
  23 reghash_t *stk_reg_read(const char *fn)
  24 {
  25         reghash_t *h = kh_init(reg);
  26         gzFile fp;
  27         kstream_t *ks;
  28         int dret;
  29         kstring_t *str;
  30         // read the list
  31         str = calloc(1, sizeof(kstring_t));
  32         fp = strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
  33         ks = ks_init(fp);
  34         while (ks_getuntil(ks, 0, str, &dret) >= 0) {
  35                 int beg = -1, end = -1;
  36                 reglist_t *p;
  37                 khint_t k = kh_get(reg, h, str->s);
  38                 if (k == kh_end(h)) {
  39                         int ret;
  40                         char *s = strdup(str->s);
  41                         k = kh_put(reg, h, s, &ret);
  42                         memset(&kh_val(h, k), 0, sizeof(reglist_t));
  43                 }
  44                 p = &kh_val(h, k);
  45                 if (dret != '\n') {
  46                         if (ks_getuntil(ks, 0, str, &dret) > 0 && isdigit(str->s[0])) {
  47                                 beg = atoi(str->s);
  48                                 if (dret != '\n') {
  49                                         if (ks_getuntil(ks, 0, str, &dret) > 0 && isdigit(str->s[0])) {
  50                                                 end = atoi(str->s);
  51                                                 if (end < 0) end = -1;
  52                                         }
  53                                 }
  54                         }
  55                 }
  56                 // skip the rest of the line
  57                 if (dret != '\n') while ((dret = ks_getc(ks)) > 0 && dret != '\n');
  58                 if (end < 0 && beg > 0) end = beg, beg = beg - 1; // if there is only one column
  59                 if (beg < 0) beg = 0, end = INT_MAX;
  60                 if (p->n == p->m) {
  61                         p->m = p->m? p->m<<1 : 4;
  62                         p->a = realloc(p->a, p->m * 8);
  63                 }
  64                 p->a[p->n++] = (uint64_t)beg<<32 | end;
  65         }
  66         ks_destroy(ks);
  67         gzclose(fp);
  68         free(str->s); free(str);
  69         return h;
  70 }
  71
  72 void stk_reg_destroy(reghash_t *h)
  73 {
  74         khint_t k;
  75         for (k = 0; k < kh_end(h); ++k) {
  76                 if (kh_exist(h, k)) {
  77                         free(kh_val(h, k).a);
  78                         free((char*)kh_key(h, k));
  79                 }
  80         }
  81         kh_destroy(reg, h);
  82 }
  83
  84 /* constant table */
  85
  86 unsigned char seq_nt16_table[256] = {
  87         15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
  88         15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
  89         15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15 /*'-'*/,15,15,
  90         15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
  91         15, 1,14, 2, 13,15,15, 4, 11,15,15,12, 15, 3,15,15,
  92         15,15, 5, 6,  8,15, 7, 9,  0,10,15,15, 15,15,15,15,
  93         15, 1,14, 2, 13,15,15, 4, 11,15,15,12, 15, 3,15,15,
  94         15,15, 5, 6,  8,15, 7, 9,  0,10,15,15, 15,15,15,15,
  95         15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
  96         15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
  97         15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
  98         15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
  99         15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
 100         15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
 101         15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
 102         15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15
 103 };
 104
 105 char *seq_nt16_rev_table = "XACMGRSVTWYHKDBN";
 106 unsigned char seq_nt16to4_table[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 };
 107 int bitcnt_table[] = { 4, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4 };
 108
 109 /* composition */
 110 int stk_comp(int argc, char *argv[])
 111 {
 112         gzFile fp;
 113         kseq_t *seq;
 114         int l, c, upper_only = 0;
 115         reghash_t *h = 0;
 116         reglist_t dummy;
 117         while ((c = getopt(argc, argv, "ur:")) >= 0) {
 118                 switch (c) {
 119                         case 'u': upper_only = 1; break;
 120                         case 'r': h = stk_reg_read(optarg); break;
 121                 }
 122         }
 123         if (argc == optind) {
 124                 fprintf(stderr, "Usage:  seqtk comp [-u] [-r in.bed] <in.fa>\n\n");
 125                 fprintf(stderr, "Output format: chr, length, #A, #C, #G, #T, #2, #3, #4, #CpG, #tv, #ts, #CpG-ts\n");
 126                 return 1;
 127         }
 128         fp = (strcmp(argv[optind], "-") == 0)? gzdopen(fileno(stdin), "r") : gzopen(argv[optind], "r");
 129         seq = kseq_init(fp);
 130         dummy.n= dummy.m = 1; dummy.a = calloc(1, 8);
 131         while ((l = kseq_read(seq)) >= 0) {
 132                 int i, k;
 133                 reglist_t *p = 0;
 134                 if (h) {
 135                         khint_t k = kh_get(reg, h, seq->name.s);
 136                         if (k != kh_end(h)) p = &kh_val(h, k);
 137                 } else {
 138                         p = &dummy;
 139                         dummy.a[0] = l;
 140                 }
 141                 for (k = 0; p && k < p->n; ++k) {
 142                         int beg = p->a[k]>>32, end = p->a[k]&0xffffffff;
 143                         int la, lb, lc, na, nb, nc, cnt[11];
 144                         if (beg > 0) la = seq->seq.s[beg-1], lb = seq_nt16_table[la], lc = bitcnt_table[lb];
 145                         else la = 'a', lb = -1, lc = 0;
 146                         na = seq->seq.s[beg]; nb = seq_nt16_table[na]; nc = bitcnt_table[nb];
 147                         memset(cnt, 0, 11 * sizeof(int));
 148                         for (i = beg; i < end; ++i) {
 149                                 int is_CpG = 0, a, b, c;
 150                                 a = na; b = nb; c = nc;
 151                                 na = seq->seq.s[i+1]; nb = seq_nt16_table[na]; nc = bitcnt_table[nb];
 152                                 if (b == 2 || b == 10) { // C or Y
 153                                         if (nb == 4 || nb == 5) is_CpG = 1;
 154                                 } else if (b == 4 || b == 5) { // G or R
 155                                         if (lb == 2 || lb == 10) is_CpG = 1;
 156                                 }
 157                                 if (upper_only == 0 || isupper(a)) {
 158                                         if (c > 1) ++cnt[c+2];
 159                                         if (c == 1) ++cnt[seq_nt16to4_table[b]];
 160                                         if (b == 10 || b == 5) ++cnt[9];
 161                                         else if (c == 2) {
 162                                                 ++cnt[8];
 163                                         }
 164                                         if (is_CpG) {
 165                                                 ++cnt[7];
 166                                                 if (b == 10 || b == 5) ++cnt[10];
 167                                         }
 168                                 }
 169                                 la = a; lb = b; lc = c;
 170                         }
 171                         if (h) printf("%s\t%d\t%d", seq->name.s, beg, end);
 172                         else printf("%s\t%d", seq->name.s, l);
 173                         for (i = 0; i < 11; ++i) printf("\t%d", cnt[i]);
 174                         putchar('\n');
 175                 }
 176                 fflush(stdout);
 177         }
 178         free(dummy.a);
 179         kseq_destroy(seq);
 180         gzclose(fp);
 181         return 0;
 182 }
 183
 184 int stk_randbase(int argc, char *argv[])
 185 {
 186         gzFile fp;
 187         kseq_t *seq;
 188         int l;
 189         if (argc == 1) {
 190                 fprintf(stderr, "Usage: seqtk randbase <in.fa>\n");
 191                 return 1;
 192         }
 193         fp = (strcmp(argv[1], "-") == 0)? gzdopen(fileno(stdin), "r") : gzopen(argv[1], "r");
 194         seq = kseq_init(fp);
 195         while ((l = kseq_read(seq)) >= 0) {
 196                 int i;
 197                 printf(">%s", seq->name.s);
 198                 for (i = 0; i < l; ++i) {
 199                         int c, b, a, j, k, m;
 200                         b = seq->seq.s[i];
 201                         c = seq_nt16_table[b];
 202                         a = bitcnt_table[c];
 203                         if (a == 2) {
 204                                 m = (drand48() < 0.5);
 205                                 for (j = k = 0; j < 4; ++j) {
 206                                         if ((1<<j & c) == 0) continue;
 207                                         if (k == m) break;
 208                                         ++k;
 209                                 }
 210                                 seq->seq.s[i] = islower(b)? "acgt"[j] : "ACGT"[j];
 211                         }
 212                         if (i%60 == 0) putchar('\n');
 213                         putchar(seq->seq.s[i]);
 214                 }
 215                 putchar('\n');
 216         }
 217         kseq_destroy(seq);
 218         gzclose(fp);
 219         return 0;
 220 }
 221
 222 int stk_hety(int argc, char *argv[])
 223 {
 224         gzFile fp;
 225         kseq_t *seq;
 226         int l, c, win_size = 50000, n_start = 5, win_step, is_lower_mask = 0;
 227         char *buf;
 228         uint32_t cnt[3];
 229         if (argc == 1) {
 230                 fprintf(stderr, "\n");
 231                 fprintf(stderr, "Usage:   seqtk hety [options] <in.fa>\n\n");
 232                 fprintf(stderr, "Options: -w INT   window size [%d]\n", win_size);
 233                 fprintf(stderr, "         -t INT   # start positions in a window [%d]\n", n_start);
 234                 fprintf(stderr, "         -m       treat lowercases as masked\n");
 235                 fprintf(stderr, "\n");
 236                 return 1;
 237         }
 238         while ((c = getopt(argc, argv, "w:t:m")) >= 0) {
 239                 switch (c) {
 240                 case 'w': win_size = atoi(optarg); break;
 241                 case 't': n_start = atoi(optarg); break;
 242                 case 'm': is_lower_mask = 1; break;
 243                 }
 244         }
 245         fp = (strcmp(argv[optind], "-") == 0)? gzdopen(fileno(stdin), "r") : gzopen(argv[optind], "r");
 246         seq = kseq_init(fp);
 247         win_step = win_size / n_start;
 248         buf = calloc(win_size, 1);
 249         while ((l = kseq_read(seq)) >= 0) {
 250                 int x, i, y, z, next = 0;
 251                 cnt[0] = cnt[1] = cnt[2] = 0;
 252                 for (i = 0; i <= l; ++i) {
 253                         if ((i >= win_size && i % win_step == 0) || i == l) {
 254                                 if (i == l && l >= win_size) {
 255                                         for (y = l - win_size; y < next; ++y) --cnt[(int)buf[y % win_size]];
 256                                 }
 257                                 if (cnt[1] + cnt[2] > 0)
 258                                         printf("%s\t%d\t%d\t%.2lf\t%d\t%d\n", seq->name.s, next, i,
 259                                                    (double)cnt[2] / (cnt[1] + cnt[2]) * win_size, cnt[1] + cnt[2], cnt[2]);
 260                                 next = i;
 261                         }
 262                         if (i < l) {
 263                                 y = i % win_size;
 264                                 c = seq->seq.s[i];
 265                                 if (is_lower_mask && islower(c)) c = 'N';
 266                                 c = seq_nt16_table[c];
 267                                 x = bitcnt_table[c];
 268                                 if (i >= win_size) --cnt[(int)buf[y]];
 269                                 buf[y] = z = x > 2? 0 : x == 2? 2 : 1;
 270                                 ++cnt[z];
 271                         }
 272                 }
 273         }
 274         free(buf);
 275         kseq_destroy(seq);
 276         gzclose(fp);
 277         return 0;
 278 }
 279
 280 /* fq2fa */
 281 int stk_fq2fa(int argc, char *argv[])
 282 {
 283         gzFile fp;
 284         kseq_t *seq;
 285         char *buf;
 286         int l, i, c, qual_thres = 0, linelen = 60;
 287         while ((c = getopt(argc, argv, "q:l:")) >= 0) {
 288                 switch (c) {
 289                         case 'q': qual_thres = atoi(optarg); break;
 290                         case 'l': linelen = atoi(optarg); break;
 291                 }
 292         }
 293         if (argc == optind) {
 294                 fprintf(stderr, "Usage: seqtk fq2fa [-q qualThres=0] [-l lineLen=60] <in.fq>\n");
 295                 return 1;
 296         }
 297         buf = linelen > 0? malloc(linelen + 1) : 0;
 298         fp = strcmp(argv[optind], "-")? gzopen(argv[optind], "r") : gzdopen(fileno(stdin), "r");
 299         seq = kseq_init(fp);
 300         while ((l = kseq_read(seq)) >= 0) {
 301                 if (seq->qual.l && qual_thres > 0) {
 302                         for (i = 0; i < l; ++i)
 303                                 if (seq->qual.s[i] - 33 < qual_thres)
 304                                         seq->seq.s[i] = tolower(seq->seq.s[i]);
 305                 }
 306                 putchar('>');
 307                 if (seq->comment.l) {
 308                         fputs(seq->name.s, stdout);
 309                         putchar(' ');
 310                         puts(seq->comment.s);
 311                 } else puts(seq->name.s);
 312                 if (buf) { // multi-line
 313                         for (i = 0; i < l; i += linelen) {
 314                                 int x = i + linelen < l? linelen : l - i;
 315                                 memcpy(buf, seq->seq.s + i, x);
 316                                 buf[x] = 0;
 317                                 puts(buf);
 318                         }
 319                 } else puts(seq->seq.s);
 320         }
 321         free(buf);
 322         kseq_destroy(seq);
 323         gzclose(fp);
 324         return 0;
 325 }
 326
 327 int stk_maskseq(int argc, char *argv[])
 328 {
 329         khash_t(reg) *h = kh_init(reg);
 330         gzFile fp;
 331         kseq_t *seq;
 332         int l, i, j, c, is_complement = 0, is_lower = 0;
 333         khint_t k;
 334         while ((c = getopt(argc, argv, "cl")) >= 0) {
 335                 switch (c) {
 336                 case 'c': is_complement = 1; break;
 337                 case 'l': is_lower = 1; break;
 338                 }
 339         }
 340         if (argc - optind < 2) {
 341                 fprintf(stderr, "Usage:   seqtk maskseq [-cl] <in.fa> <in.bed>\n\n");
 342                 fprintf(stderr, "Options: -c     mask the complement regions\n");
 343                 fprintf(stderr, "         -l     soft mask (to lower cases)\n");
 344                 return 1;
 345         }
 346         h = stk_reg_read(argv[optind+1]);
 347         // maskseq
 348         fp = strcmp(argv[optind], "-")? gzopen(argv[optind], "r") : gzdopen(fileno(stdin), "r");
 349         seq = kseq_init(fp);
 350         while ((l = kseq_read(seq)) >= 0) {
 351                 k = kh_get(reg, h, seq->name.s);
 352                 if (k == kh_end(h)) { // not found in the hash table
 353                         if (is_complement) {
 354                                 for (j = 0; j < l; ++j)
 355                                         seq->seq.s[j] = is_lower? tolower(seq->seq.s[j]) : 'N';
 356                         }
 357                 } else {
 358                         reglist_t *p = &kh_val(h, k);
 359                         if (!is_complement) {
 360                                 for (i = 0; i < p->n; ++i) {
 361                                         int beg = p->a[i]>>32, end = p->a[i];
 362                                         if (beg >= seq->seq.l) {
 363                                                 fprintf(stderr, "[maskseq] start position >= the sequence length.\n");
 364                                                 continue;
 365                                         }
 366                                         if (end >= seq->seq.l) end = seq->seq.l;
 367                                         if (is_lower) for (j = beg; j < end; ++j) seq->seq.s[j] = tolower(seq->seq.s[j]);
 368                                         else for (j = beg; j < end; ++j) seq->seq.s[j] = 'N';
 369                                 }
 370                         } else {
 371                                 int8_t *mask = calloc(seq->seq.l, 1);
 372                                 for (i = 0; i < p->n; ++i) {
 373                                         int beg = p->a[i]>>32, end = p->a[i];
 374                                         if (end >= seq->seq.l) end = seq->seq.l;
 375                                         for (j = beg; j < end; ++j) mask[j] = 1;
 376                                 }
 377                                 for (j = 0; j < l; ++j)
 378                                         if (mask[j] == 0) seq->seq.s[j] = is_lower? tolower(seq->seq.s[j]) : 'N';
 379                                 free(mask);
 380                         }
 381                 }
 382                 printf(">%s", seq->name.s);
 383                 for (j = 0; j < seq->seq.l; ++j) {
 384                         if (j%60 == 0) putchar('\n');
 385                         putchar(seq->seq.s[j]);
 386                 }
 387                 putchar('\n');
 388         }
 389         // free
 390         kseq_destroy(seq);
 391         gzclose(fp);
 392         stk_reg_destroy(h);
 393         return 0;
 394 }
 395
 396 /* subseq */
 397
 398 int stk_subseq(int argc, char *argv[])
 399 {
 400         khash_t(reg) *h = kh_init(reg);
 401         gzFile fp;
 402         kseq_t *seq;
 403         int l, i, j, c, is_tab = 0;
 404         khint_t k;
 405         while ((c = getopt(argc, argv, "t")) >= 0) {
 406                 switch (c) {
 407                 case 't': is_tab = 1; break;
 408                 }
 409         }
 410         if (optind + 2 > argc) {
 411                 fprintf(stderr, "Usage: seqtk subseq [-t] <in.fa> <in.bed>\n\n");
 412                 fprintf(stderr, "Note: Use 'samtools faidx' if only a few regions are intended.\n");
 413                 return 1;
 414         }
 415         h = stk_reg_read(argv[optind+1]);
 416         // subseq
 417         fp = strcmp(argv[optind], "-")? gzopen(argv[optind], "r") : gzdopen(fileno(stdin), "r");
 418         seq = kseq_init(fp);
 419         while ((l = kseq_read(seq)) >= 0) {
 420                 reglist_t *p;
 421                 k = kh_get(reg, h, seq->name.s);
 422                 if (k == kh_end(h)) continue;
 423                 p = &kh_val(h, k);
 424                 for (i = 0; i < p->n; ++i) {
 425                         int beg = p->a[i]>>32, end = p->a[i];
 426                         if (beg >= seq->seq.l) {
 427                                 fprintf(stderr, "[subseq] %s: %d >= %ld\n", seq->name.s, beg, seq->seq.l);
 428                                 continue;
 429                         }
 430                         if (end > seq->seq.l) end = seq->seq.l;
 431                         if (is_tab == 0) {
 432                                 printf("%c%s", seq->qual.l == seq->seq.l? '@' : '>', seq->name.s);
 433                                 if (end == INT_MAX) {
 434                                         if (beg) printf(":%d", beg+1);
 435                                 } else printf(":%d-%d", beg+1, end);
 436                         } else printf("%s\t%d\t", seq->name.s, beg + 1);
 437                         if (end > seq->seq.l) end = seq->seq.l;
 438                         for (j = 0; j < end - beg; ++j) {
 439                                 if (is_tab == 0 && j % 60 == 0) putchar('\n');
 440                                 putchar(seq->seq.s[j + beg]);
 441                         }
 442                         putchar('\n');
 443                         if (seq->qual.l != seq->seq.l || is_tab) continue;
 444                         printf("+");
 445                         for (j = 0; j < end - beg; ++j) {
 446                                 if (j % 60 == 0) putchar('\n');
 447                                 putchar(seq->qual.s[j + beg]);
 448                         }
 449                         putchar('\n');
 450                 }
 451         }
 452         // free
 453         kseq_destroy(seq);
 454         gzclose(fp);
 455         stk_reg_destroy(h);
 456         return 0;
 457 }
 458
 459 /* mergefa */
 460 int stk_mergefa(int argc, char *argv[])
 461 {
 462         gzFile fp[2];
 463         kseq_t *seq[2];
 464         int i, l, c, is_intersect = 0, is_haploid = 0, qual = 0, is_mask = 0;
 465         while ((c = getopt(argc, argv, "himq:")) >= 0) {
 466                 switch (c) {
 467                         case 'i': is_intersect = 1; break;
 468                         case 'h': is_haploid = 1; break;
 469                         case 'm': is_mask = 1; break;
 470                         case 'q': qual = atoi(optarg); break;
 471                 }
 472         }
 473         if (is_mask && is_intersect) {
 474                 fprintf(stderr, "[%s] `-i' and `-h' cannot be applied at the same time.\n", __func__);
 475                 return 1;
 476         }
 477         if (optind + 2 > argc) {
 478                 fprintf(stderr, "\nUsage: seqtk mergefa [options] <in1.fa> <in2.fa>\n\n");
 479                 fprintf(stderr, "Options: -q INT   quality threshold [0]\n");
 480                 fprintf(stderr, "         -i       take intersection\n");
 481                 fprintf(stderr, "         -m       convert to lowercase when one of the input base is N.\n");
 482                 fprintf(stderr, "         -h       suppress hets in the input\n\n");
 483                 return 1;
 484         }
 485         for (i = 0; i < 2; ++i) {
 486                 fp[i] = strcmp(argv[optind+i], "-")? gzopen(argv[optind+i], "r") : gzdopen(fileno(stdin), "r");
 487                 seq[i] = kseq_init(fp[i]);
 488         }
 489         while (kseq_read(seq[0]) >= 0) {
 490                 int min_l, c[2], is_upper;
 491                 kseq_read(seq[1]);
 492                 if (strcmp(seq[0]->name.s, seq[1]->name.s))
 493                         fprintf(stderr, "[%s] Different sequence names: %s != %s\n", __func__, seq[0]->name.s, seq[1]->name.s);
 494                 if (seq[0]->seq.l != seq[1]->seq.l)
 495                         fprintf(stderr, "[%s] Unequal sequence length: %ld != %ld\n", __func__, seq[0]->seq.l, seq[1]->seq.l);
 496                 min_l = seq[0]->seq.l < seq[1]->seq.l? seq[0]->seq.l : seq[1]->seq.l;
 497                 printf(">%s", seq[0]->name.s);
 498                 for (l = 0; l < min_l; ++l) {
 499                         c[0] = seq[0]->seq.s[l]; c[1] = seq[1]->seq.s[l];
 500                         if (seq[0]->qual.l && seq[0]->qual.s[l] - 33 < qual) c[0] = tolower(c[0]);
 501                         if (seq[1]->qual.l && seq[1]->qual.s[l] - 33 < qual) c[1] = tolower(c[1]);
 502                         if (is_intersect) is_upper = (isupper(c[0]) || isupper(c[1]))? 1 : 0;
 503                         else if (is_mask) is_upper = (isupper(c[0]) || isupper(c[1]))? 1 : 0;
 504                         else is_upper = (isupper(c[0]) && isupper(c[1]))? 1 : 0;
 505                         c[0] = seq_nt16_table[c[0]]; c[1] = seq_nt16_table[c[1]];
 506                         if (c[0] == 0) c[0] = 15;
 507                         if (c[1] == 0) c[1] = 15;
 508                         if (is_haploid && (bitcnt_table[c[0]] > 1 || bitcnt_table[c[1]] > 1)) is_upper = 0;
 509                         if (is_intersect) {
 510                                 c[0] = c[0] & c[1];
 511                                 if (c[0] == 0) is_upper = 0;
 512                         } else if (is_mask) {
 513                                 if (c[0] == 15 || c[1] == 15) is_upper = 0;
 514                                 c[0] = c[0] & c[1];
 515                                 if (c[0] == 0) is_upper = 0;
 516                         } else c[0] = c[0] | c[1];
 517                         c[0] = seq_nt16_rev_table[c[0]];
 518                         if (!is_upper) c[0] = tolower(c[0]);
 519                         if (l%60 == 0) putchar('\n');
 520                         putchar(c[0]);
 521                 }
 522                 putchar('\n');
 523         }
 524         return 0;
 525 }
 526
 527 int stk_famask(int argc, char *argv[])
 528 {
 529         gzFile fp[2];
 530         kseq_t *seq[2];
 531         int i, l;
 532         if (argc < 3) {
 533                 fprintf(stderr, "Usage: seqtk famask <src.fa> <mask.fa>\n");
 534                 return 1;
 535         }
 536         for (i = 0; i < 2; ++i) {
 537                 fp[i] = strcmp(argv[optind+i], "-")? gzopen(argv[optind+i], "r") : gzdopen(fileno(stdin), "r");
 538                 seq[i] = kseq_init(fp[i]);
 539         }
 540         while (kseq_read(seq[0]) >= 0) {
 541                 int min_l, c[2];
 542                 kseq_read(seq[1]);
 543                 if (strcmp(seq[0]->name.s, seq[1]->name.s))
 544                         fprintf(stderr, "[%s] Different sequence names: %s != %s\n", __func__, seq[0]->name.s, seq[1]->name.s);
 545                 if (seq[0]->seq.l != seq[1]->seq.l)
 546                         fprintf(stderr, "[%s] Unequal sequence length: %ld != %ld\n", __func__, seq[0]->seq.l, seq[1]->seq.l);
 547                 min_l = seq[0]->seq.l < seq[1]->seq.l? seq[0]->seq.l : seq[1]->seq.l;
 548                 printf(">%s", seq[0]->name.s);
 549                 for (l = 0; l < min_l; ++l) {
 550                         c[0] = seq[0]->seq.s[l]; c[1] = seq[1]->seq.s[l];
 551                         if (c[1] == 'x') c[0] = tolower(c[0]);
 552                         else if (c[1] != 'X') c[0] = c[1];
 553                         if (l%60 == 0) putchar('\n');
 554                         putchar(c[0]);
 555                 }
 556                 putchar('\n');
 557         }
 558         return 0;
 559 }
 560
 561 int stk_mutfa(int argc, char *argv[])
 562 {
 563         khash_t(reg) *h = kh_init(reg);
 564         gzFile fp;
 565         kseq_t *seq;
 566         kstream_t *ks;
 567         int l, i, dret;
 568         kstring_t *str;
 569         khint_t k;
 570         if (argc < 3) {
 571                 fprintf(stderr, "Usage: seqtk mutfa <in.fa> <in.snp>\n\n");
 572                 fprintf(stderr, "Note: <in.snp> contains at least four columns per line which are:\n");
 573                 fprintf(stderr, "      'chr  1-based-pos  any  base-changed-to'.\n");
 574                 return 1;
 575         }
 576         // read the list
 577         str = calloc(1, sizeof(kstring_t));
 578         fp = strcmp(argv[2], "-")? gzopen(argv[2], "r") : gzdopen(fileno(stdin), "r");
 579         ks = ks_init(fp);
 580         while (ks_getuntil(ks, 0, str, &dret) >= 0) {
 581                 char *s = strdup(str->s);
 582                 int beg = 0, ret;
 583                 reglist_t *p;
 584                 k = kh_get(reg, h, s);
 585                 if (k == kh_end(h)) {
 586                         k = kh_put(reg, h, s, &ret);
 587                         memset(&kh_val(h, k), 0, sizeof(reglist_t));
 588                 }
 589                 p = &kh_val(h, k);
 590                 if (ks_getuntil(ks, 0, str, &dret) > 0) beg = atol(str->s) - 1; // 2nd col
 591                 ks_getuntil(ks, 0, str, &dret); // 3rd col
 592                 ks_getuntil(ks, 0, str, &dret); // 4th col
 593                 // skip the rest of the line
 594                 if (dret != '\n') while ((dret = ks_getc(ks)) > 0 && dret != '\n');
 595                 if (isalpha(str->s[0]) && str->l == 1) {
 596                         if (p->n == p->m) {
 597                                 p->m = p->m? p->m<<1 : 4;
 598                                 p->a = realloc(p->a, p->m * 8);
 599                         }
 600                         p->a[p->n++] = (uint64_t)beg<<32 | str->s[0];
 601                 }
 602         }
 603         ks_destroy(ks);
 604         gzclose(fp);
 605         free(str->s); free(str);
 606         // mutfa
 607         fp = strcmp(argv[1], "-")? gzopen(argv[1], "r") : gzdopen(fileno(stdin), "r");
 608         seq = kseq_init(fp);
 609         while ((l = kseq_read(seq)) >= 0) {
 610                 reglist_t *p;
 611                 k = kh_get(reg, h, seq->name.s);
 612                 if (k != kh_end(h)) {
 613                         p = &kh_val(h, k);
 614                         for (i = 0; i < p->n; ++i) {
 615                                 int beg = p->a[i]>>32;
 616                                 if (beg < seq->seq.l)
 617                                         seq->seq.s[beg] = (int)p->a[i];
 618                         }
 619                 }
 620                 printf(">%s", seq->name.s);
 621                 for (i = 0; i < l; ++i) {
 622                         if (i%60 == 0) putchar('\n');
 623                         putchar(seq->seq.s[i]);
 624                 }
 625                 putchar('\n');
 626         }
 627         // free
 628         kseq_destroy(seq);
 629         gzclose(fp);
 630         for (k = 0; k < kh_end(h); ++k) {
 631                 if (kh_exist(h, k)) {
 632                         free(kh_val(h, k).a);
 633                         free((char*)kh_key(h, k));
 634                 }
 635         }
 636         kh_destroy(reg, h);
 637         return 0;
 638 }
 639
 640 int stk_listhet(int argc, char *argv[])
 641 {
 642         gzFile fp;
 643         kseq_t *seq;
 644         int i, l;
 645         if (argc == 1) {
 646                 fprintf(stderr, "Usage: seqtk listhet <in.fa>\n");
 647                 return 1;
 648         }
 649         fp = (strcmp(argv[1], "-") == 0)? gzdopen(fileno(stdin), "r") : gzopen(argv[1], "r");
 650         seq = kseq_init(fp);
 651         while ((l = kseq_read(seq)) >= 0) {
 652                 for (i = 0; i < l; ++i) {
 653                         int b = seq->seq.s[i];
 654                         if (bitcnt_table[seq_nt16_table[b]] == 2)
 655                                 printf("%s\t%d\t%c\n", seq->name.s, i+1, b);
 656                 }
 657         }
 658         kseq_destroy(seq);
 659         gzclose(fp);
 660         return 0;
 661 }
 662
 663 /* cutN */
 664 static int cutN_min_N_tract = 1000;
 665 static int cutN_nonN_penalty = 10;
 666
 667 static int find_next_cut(const kseq_t *ks, int k, int *begin, int *end)
 668 {
 669         int i, b, e;
 670         while (k < ks->seq.l) {
 671                 if (seq_nt16_table[(int)ks->seq.s[k]] == 15) {
 672                         int score, max;
 673                         score = 0; e = max = -1;
 674                         for (i = k; i < ks->seq.l && score >= 0; ++i) { /* forward */
 675                                 if (seq_nt16_table[(int)ks->seq.s[i]] == 15) ++score;
 676                                 else score -= cutN_nonN_penalty;
 677                                 if (score > max) max = score, e = i;
 678                         }
 679                         score = 0; b = max = -1;
 680                         for (i = e; i >= 0 && score >= 0; --i) { /* backward */
 681                                 if (seq_nt16_table[(int)ks->seq.s[i]] == 15) ++score;
 682                                 else score -= cutN_nonN_penalty;
 683                                 if (score > max) max = score, b = i;
 684                         }
 685                         if (e + 1 - b >= cutN_min_N_tract) {
 686                                 *begin = b;
 687                                 *end = e + 1;
 688                                 return *end;
 689                         }
 690                         k = e + 1;
 691                 } else ++k;
 692         }
 693         return -1;
 694 }
 695 static void print_seq(FILE *fpout, const kseq_t *ks, int begin, int end)
 696 {
 697         int i;
 698         if (begin >= end) return; // FIXME: why may this happen? Understand it!
 699         fprintf(fpout, ">%s:%d-%d", ks->name.s, begin+1, end);
 700         for (i = begin; i < end && i < ks->seq.l; ++i) {
 701                 if ((i - begin)%60 == 0) fputc('\n', fpout);
 702                 fputc(ks->seq.s[i], fpout);
 703         }
 704         fputc('\n', fpout);
 705 }
 706 int stk_cutN(int argc, char *argv[])
 707 {
 708         int c, l, gap_only = 0;
 709         gzFile fp;
 710         kseq_t *ks;
 711         while ((c = getopt(argc, argv, "n:p:g")) >= 0) {
 712                 switch (c) {
 713                 case 'n': cutN_min_N_tract = atoi(optarg); break;
 714                 case 'p': cutN_nonN_penalty = atoi(optarg); break;
 715                 case 'g': gap_only = 1; break;
 716                 default: return 1;
 717                 }
 718         }
 719         if (argc == optind) {
 720                 fprintf(stderr, "\n");
 721                 fprintf(stderr, "Usage:   seqtk cutN [options] <in.fa>\n\n");
 722                 fprintf(stderr, "Options: -n INT    min size of N tract [%d]\n", cutN_min_N_tract);
 723                 fprintf(stderr, "         -p INT    penalty for a non-N [%d]\n", cutN_nonN_penalty);
 724                 fprintf(stderr, "         -g        print gaps only, no sequence\n\n");
 725                 return 1;
 726         }
 727         fp = (strcmp(argv[optind], "-") == 0)? gzdopen(fileno(stdin), "r") : gzopen(argv[optind], "r");
 728         ks = kseq_init(fp);
 729         while ((l = kseq_read(ks)) >= 0) {
 730                 int k = 0, begin = 0, end = 0;
 731                 while (find_next_cut(ks, k, &begin, &end) >= 0) {
 732                         if (begin != 0) {
 733                                 if (gap_only) printf("%s\t%d\t%d\n", ks->name.s, begin, end);
 734                                 else print_seq(stdout, ks, k, begin);
 735                         }
 736                         k = end;
 737                 }
 738                 if (!gap_only) print_seq(stdout, ks, k, l);
 739         }
 740         kseq_destroy(ks);
 741         gzclose(fp);
 742         return 0;
 743 }
 744
 745 /* main function */
 746 static int usage()
 747 {
 748         fprintf(stderr, "\n");
 749         fprintf(stderr, "Usage:   seqtk <command> <arguments>\n\n");
 750         fprintf(stderr, "Command: comp      get the nucleotide composite of FASTA/Q\n");
 751         fprintf(stderr, "         hety      regional heterozygosity\n");
 752         fprintf(stderr, "         fq2fa     convert FASTQ to FASTA\n");
 753         fprintf(stderr, "         subseq    extract subsequences from FASTA/Q\n");
 754         fprintf(stderr, "         maskseq   mask sequences\n");
 755         fprintf(stderr, "         mutfa     point mutate FASTA at specified positions\n");
 756         fprintf(stderr, "         mergefa   merge two FASTA/Q files\n");
 757         fprintf(stderr, "         randbase  choose a random base from hets\n");
 758         fprintf(stderr, "         cutN      cut sequence at long N\n");
 759         fprintf(stderr, "         listhet   extract the position of each het\n");
 760         fprintf(stderr, "\n");
 761         return 1;
 762 }
 763
 764 int main(int argc, char *argv[])
 765 {
 766         if (argc == 1) return usage();
 767         if (strcmp(argv[1], "comp") == 0) stk_comp(argc-1, argv+1);
 768         else if (strcmp(argv[1], "hety") == 0) stk_hety(argc-1, argv+1);
 769         else if (strcmp(argv[1], "fq2fa") == 0) stk_fq2fa(argc-1, argv+1);
 770         else if (strcmp(argv[1], "subseq") == 0) stk_subseq(argc-1, argv+1);
 771         else if (strcmp(argv[1], "maskseq") == 0) stk_maskseq(argc-1, argv+1);
 772         else if (strcmp(argv[1], "mutfa") == 0) stk_mutfa(argc-1, argv+1);
 773         else if (strcmp(argv[1], "mergefa") == 0) stk_mergefa(argc-1, argv+1);
 774         else if (strcmp(argv[1], "randbase") == 0) stk_randbase(argc-1, argv+1);
 775         else if (strcmp(argv[1], "cutN") == 0) stk_cutN(argc-1, argv+1);
 776         else if (strcmp(argv[1], "listhet") == 0) stk_listhet(argc-1, argv+1);
 777         else if (strcmp(argv[1], "famask") == 0) stk_famask(argc-1, argv+1);
 778         else {
 779                 fprintf(stderr, "[main] unrecognized commad '%s'. Abort!\n", argv[1]);
 780                 return 1;
 781         }
 782         return 0;
 783 }