samtools/misc/seqtk.c.pysam.c

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