Skip to content
Snippets Groups Projects
Commit 46cd37eb authored by Heng Li's avatar Heng Li
Browse files

r87: trim to a fixed length (#9)

parent 7d809104
Branches
Tags
No related merge requests found
...@@ -280,13 +280,14 @@ int stk_trimfq(int argc, char *argv[]) ...@@ -280,13 +280,14 @@ int stk_trimfq(int argc, char *argv[])
gzFile fp; gzFile fp;
kseq_t *seq; kseq_t *seq;
double param = 0.05, q_int2real[128]; double param = 0.05, q_int2real[128];
int i, c, min_len = 30, left = 0, right = 0; int i, c, min_len = 30, left = 0, right = 0, fixed_len = -1;
while ((c = getopt(argc, argv, "l:q:b:e:")) >= 0) { while ((c = getopt(argc, argv, "l:q:b:e:L:")) >= 0) {
switch (c) { switch (c) {
case 'q': param = atof(optarg); break; case 'q': param = atof(optarg); break;
case 'l': min_len = atoi(optarg); break; case 'l': min_len = atoi(optarg); break;
case 'b': left = atoi(optarg); break; case 'b': left = atoi(optarg); break;
case 'e': right = atoi(optarg); break; case 'e': right = atoi(optarg); break;
case 'L': fixed_len = atoi(optarg); break;
} }
} }
if (optind == argc) { if (optind == argc) {
...@@ -296,6 +297,7 @@ int stk_trimfq(int argc, char *argv[]) ...@@ -296,6 +297,7 @@ int stk_trimfq(int argc, char *argv[])
fprintf(stderr, " -l INT maximally trim down to INT bp (disabled by -b/-e) [%d]\n", min_len); fprintf(stderr, " -l INT maximally trim down to INT bp (disabled by -b/-e) [%d]\n", min_len);
fprintf(stderr, " -b INT trim INT bp from left (non-zero to disable -q/-l) [0]\n"); fprintf(stderr, " -b INT trim INT bp from left (non-zero to disable -q/-l) [0]\n");
fprintf(stderr, " -e INT trim INT bp from right (non-zero to disable -q/-l) [0]\n"); fprintf(stderr, " -e INT trim INT bp from right (non-zero to disable -q/-l) [0]\n");
fprintf(stderr, " -L INT retain at most INT bp from the 5'-end (non-zero to disable -q/-l) [0]\n");
fprintf(stderr, "\n"); fprintf(stderr, "\n");
return 1; return 1;
} }
...@@ -306,9 +308,10 @@ int stk_trimfq(int argc, char *argv[]) ...@@ -306,9 +308,10 @@ int stk_trimfq(int argc, char *argv[])
while (kseq_read(seq) >= 0) { while (kseq_read(seq) >= 0) {
int beg, tmp, end; int beg, tmp, end;
double s, max; double s, max;
if (left || right) { if (left || right || fixed_len > 0) {
beg = left; end = seq->seq.l - right; beg = left; end = seq->seq.l - right;
if (beg >= end) beg = end = 0; if (beg >= end) beg = end = 0;
if (end - beg > fixed_len) end = beg + fixed_len;
} else if (seq->qual.l > min_len) { } else if (seq->qual.l > min_len) {
for (i = 0, beg = tmp = 0, end = seq->qual.l, s = max = 0.; i < seq->qual.l; ++i) { for (i = 0, beg = tmp = 0, end = seq->qual.l, s = max = 0.; i < seq->qual.l; ++i) {
int q = seq->qual.s[i]; int q = seq->qual.s[i];
...@@ -1555,7 +1558,7 @@ static int usage() ...@@ -1555,7 +1558,7 @@ static int usage()
{ {
fprintf(stderr, "\n"); fprintf(stderr, "\n");
fprintf(stderr, "Usage: seqtk <command> <arguments>\n"); fprintf(stderr, "Usage: seqtk <command> <arguments>\n");
fprintf(stderr, "Version: 1.0-r86-dirty\n\n"); fprintf(stderr, "Version: 1.0-r87-dirty\n\n");
fprintf(stderr, "Command: seq common transformation of FASTA/Q\n"); fprintf(stderr, "Command: seq common transformation of FASTA/Q\n");
fprintf(stderr, " comp get the nucleotide composition of FASTA/Q\n"); fprintf(stderr, " comp get the nucleotide composition of FASTA/Q\n");
fprintf(stderr, " sample subsample sequences\n"); fprintf(stderr, " sample subsample sequences\n");
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment