-
Notifications
You must be signed in to change notification settings - Fork 0
/
filter_dup_reads.c
80 lines (61 loc) · 1.74 KB
/
filter_dup_reads.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include <zlib.h>
#include <stdio.h>
#include <glib.h>
#include "util.h"
#include "memutil.h"
#define MAX_LINE 1024
int main(int argc, char **argv) {
gzFile gzf;
char *input_filename, *seq_str, *count_ptr;
char line1[MAX_LINE], line2[MAX_LINE], line3[MAX_LINE], line4[MAX_LINE];
GHashTable *hashtab;
long *val, count;
if(argc != 2) {
fprintf(stderr, "usage: %s my_fastq_file.txt.gz | "
"gzip > filtered_fastq.txt.gz\n", argv[0]);
exit(2);
}
input_filename = argv[1];
gzf = util_must_gzopen(input_filename, "rb");
hashtab = g_hash_table_new_full(g_str_hash, g_str_equal,
free, free);
count = 0;
while(TRUE) {
/* read all 4 fastq lines at once time */
if((gzgets(gzf, line1, MAX_LINE) == NULL) ||
(gzgets(gzf, line2, MAX_LINE) == NULL) ||
(gzgets(gzf, line3, MAX_LINE) == NULL) ||
(gzgets(gzf, line4, MAX_LINE) == NULL)) {
break;
}
count += 1;
if(count == 100000) {
count = 0;
fprintf(stderr, ".");
}
/* the second line is the sequence, we want to check if we've seen
* the same one before
*/
seq_str = &line2[1];
val = g_hash_table_lookup(hashtab, seq_str);
if(val == NULL) {
/* we have not seen a read with this sequence before */
/* record the read in the hash table */
val = my_new(long, 1);
*val = 1;
g_hash_table_insert(hashtab, util_str_dup(seq_str), val);
/* print the read to stdout in fastq format */
fprintf(stdout, "%s%s%s%s", line1, line2, line3, line4);
} else {
/* we have seen this sequence before, increment count and discard */
*val += 1;
/* fprintf(stderr, "discarding duplicate read (seen %d times) %s \n",
* seq_str, *val);
*/
}
}
fprintf(stderr, "\n");
gzclose(gzf);
g_hash_table_destroy(hashtab);
return 0;
}