|
6
|
1 #include <stdio.h>
|
|
|
2 #include <stdlib.h>
|
|
|
3 #include <assert.h>
|
|
|
4 #include <unistd.h>
|
|
|
5 #include "const.h"
|
|
|
6 #include "maqmap.h"
|
|
|
7
|
|
|
8 maqmap_t *maq_new_maqmap()
|
|
|
9 {
|
|
|
10 maqmap_t *mm = (maqmap_t*)calloc(1, sizeof(maqmap_t));
|
|
|
11 mm->format = MAQMAP_FORMAT_NEW;
|
|
|
12 return mm;
|
|
|
13 }
|
|
|
14 void maq_delete_maqmap(maqmap_t *mm)
|
|
|
15 {
|
|
|
16 int i;
|
|
|
17 if (mm == 0) return;
|
|
|
18 for (i = 0; i < mm->n_ref; ++i)
|
|
|
19 free(mm->ref_name[i]);
|
|
|
20 free(mm->ref_name);
|
|
|
21 free(mm->mapped_reads);
|
|
|
22 free(mm);
|
|
|
23 }
|
|
|
24 void maqmap_write_header(gzFile fp, const maqmap_t *mm)
|
|
|
25 {
|
|
|
26 int i, len;
|
|
|
27 gzwrite(fp, &mm->format, sizeof(int));
|
|
|
28 gzwrite(fp, &mm->n_ref, sizeof(int));
|
|
|
29 for (i = 0; i != mm->n_ref; ++i) {
|
|
|
30 len = strlen(mm->ref_name[i]) + 1;
|
|
|
31 gzwrite(fp, &len, sizeof(int));
|
|
|
32 gzwrite(fp, mm->ref_name[i], len);
|
|
|
33 }
|
|
|
34 gzwrite(fp, &mm->n_mapped_reads, sizeof(bit64_t));
|
|
|
35 }
|
|
|
36 maqmap_t *maqmap_read_header(gzFile fp)
|
|
|
37 {
|
|
|
38 maqmap_t *mm;
|
|
|
39 int k, len;
|
|
|
40 mm = maq_new_maqmap();
|
|
|
41 gzread(fp, &mm->format, sizeof(int));
|
|
|
42 if (mm->format != MAQMAP_FORMAT_NEW) {
|
|
|
43 if (mm->format > 0) {
|
|
|
44 fprintf(stderr, "** Obsolete map format is detected. Please use 'mapass2maq' command to convert the format.\n");
|
|
|
45 exit(3);
|
|
|
46 }
|
|
|
47 assert(mm->format == MAQMAP_FORMAT_NEW);
|
|
|
48 }
|
|
|
49 gzread(fp, &mm->n_ref, sizeof(int));
|
|
|
50 mm->ref_name = (char**)calloc(mm->n_ref, sizeof(char*));
|
|
|
51 for (k = 0; k != mm->n_ref; ++k) {
|
|
|
52 gzread(fp, &len, sizeof(int));
|
|
|
53 mm->ref_name[k] = (char*)malloc(len * sizeof(char));
|
|
|
54 gzread(fp, mm->ref_name[k], len);
|
|
|
55 }
|
|
|
56 /* read number of mapped reads */
|
|
|
57 gzread(fp, &mm->n_mapped_reads, sizeof(bit64_t));
|
|
|
58 return mm;
|
|
|
59 }
|
|
|
60
|
|
|
61 /* mapvalidate */
|
|
|
62
|
|
|
63 static void mapvalidate_core(gzFile fpin)
|
|
|
64 {
|
|
|
65 maqmap_t *m = maqmap_read_header(fpin);
|
|
|
66 maqmap1_t *m1, mm1;
|
|
|
67 bit64_t n = 0;
|
|
|
68 int i, l;
|
|
|
69 bit64_t *cnt;
|
|
|
70 m1 = &mm1;
|
|
|
71 cnt = (bit64_t*)calloc(m->n_ref, 8);
|
|
|
72 printf("[message] number of reference sequences: %d\n", m->n_ref);
|
|
|
73 while ((l = maqmap_read1(fpin, m1)) != 0) {
|
|
|
74 if (l != sizeof(maqmap1_t)) {
|
|
|
75 printf("[fatal error] truncated map file.\n");
|
|
|
76 break;
|
|
|
77 }
|
|
|
78 ++n;
|
|
|
79 if ((int)m1->seqid >= m->n_ref) {
|
|
|
80 printf("[fatal error] maqmap1_t::seqid is invalid (%d >= %d).\n", m1->seqid, m->n_ref);
|
|
|
81 break;
|
|
|
82 }
|
|
|
83 ++cnt[m1->seqid];
|
|
|
84 if (m1->size >= MAX_READLEN - 1) {
|
|
|
85 printf("[faltal error] maqmap1_t::size is invalid (%d >= %d).\n", m1->size, MAX_READLEN - 1);
|
|
|
86 break;
|
|
|
87 }
|
|
|
88 }
|
|
|
89 if (m->n_mapped_reads != 0) {
|
|
|
90 if (m->n_mapped_reads != n) {
|
|
|
91 printf("[warning] maqmap1_t::n_mapped_reads is set, but not equals the real number (%llu != %llu).\n",
|
|
|
92 m->n_mapped_reads, n);
|
|
|
93 }
|
|
|
94 }
|
|
|
95 for (i = 0; i != m->n_ref; ++i)
|
|
|
96 printf("[message] %s : %llu\n", m->ref_name[i], cnt[i]);
|
|
|
97 free(cnt);
|
|
|
98 maq_delete_maqmap(m);
|
|
|
99 }
|
|
|
100
|
|
|
101 /* mapview */
|
|
|
102
|
|
|
103 static void mapview_core(FILE *fpout, gzFile fpin, int is_verbose, int is_mm)
|
|
|
104 {
|
|
|
105 bit32_t j;
|
|
|
106 maqmap_t *m = maqmap_read_header(fpin);
|
|
|
107 maqmap1_t *m1, mm1;
|
|
|
108 m1 = &mm1;
|
|
|
109 while (maqmap_read1(fpin, m1)) {
|
|
|
110 fprintf(fpout, "%s\t%s\t%d\t%c\t%d\t%u\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d",
|
|
|
111 m1->name, m->ref_name[m1->seqid], (m1->pos>>1) + 1,
|
|
|
112 (m1->pos&1)? '-' : '+', m1->dist, m1->flag, m1->map_qual, (signed char)m1->seq[MAX_READLEN-1],
|
|
|
113 m1->alt_qual, m1->info1&0xf, m1->info2, m1->c[0], m1->c[1], m1->size);
|
|
|
114 if (is_verbose) {
|
|
|
115 fputc('\t', fpout);
|
|
|
116 for (j = 0; j != m1->size; ++j) {
|
|
|
117 if (m1->seq[j] == 0) fputc('n', fpout);
|
|
|
118 else if ((m1->seq[j]&0x3f) < 27) fputc("acgt"[m1->seq[j]>>6&3], fpout);
|
|
|
119 else fputc("ACGT"[m1->seq[j]>>6&3], fpout);
|
|
|
120 }
|
|
|
121 fputc('\t', fpout);
|
|
|
122 for (j = 0; j != m1->size; ++j)
|
|
|
123 fputc((m1->seq[j]&0x3f) + 33, fpout);
|
|
|
124 }
|
|
|
125 if (is_mm) {
|
|
|
126 bit64_t *p = (bit64_t*)(m1->seq + 55);
|
|
|
127 fprintf(fpout, "\t%llx", *p);
|
|
|
128 }
|
|
|
129 fputc('\n', fpout);
|
|
|
130 }
|
|
|
131 maq_delete_maqmap(m);
|
|
|
132 }
|
|
|
133
|
|
|
134 int ma_mapview(int argc, char *argv[])
|
|
|
135 {
|
|
|
136 int c, is_verbose = 1, is_mm = 0;
|
|
|
137 while ((c = getopt(argc, argv, "bN")) >= 0) {
|
|
|
138 switch (c) {
|
|
|
139 case 'b': is_verbose = 0; break;
|
|
|
140 case 'N': is_mm = 1; break;
|
|
|
141 }
|
|
|
142 }
|
|
|
143 if (argc == optind) {
|
|
|
144 fprintf(stderr, "Usage: maq mapview [-bN] <in.map>\n");
|
|
|
145 return 1;
|
|
|
146 }
|
|
|
147 gzFile fp = (strcmp(argv[optind], "-") == 0)? gzdopen(STDIN_FILENO, "r") : gzopen(argv[optind], "r");
|
|
|
148 mapview_core(stdout, fp, is_verbose, is_mm);
|
|
|
149 gzclose(fp);
|
|
|
150 return 0;
|
|
|
151 }
|
|
|
152
|
|
|
153 int ma_mapvalidate(int argc, char *argv[])
|
|
|
154 {
|
|
|
155 gzFile fp;
|
|
|
156 if (argc < 2) {
|
|
|
157 fprintf(stderr, "Usage: maq mapvalidate <in.map>\n");
|
|
|
158 return 1;
|
|
|
159 }
|
|
|
160 fp = (strcmp(argv[optind], "-") == 0)? gzdopen(STDIN_FILENO, "r") : gzopen(argv[1], "r");
|
|
|
161 mapvalidate_core(fp);
|
|
|
162 gzclose(fp);
|
|
|
163 return 0;
|
|
|
164 }
|
