8 //#define MAQ_LONGREADS
11 # define MAX_READLEN 128
13 # define MAX_READLEN 64
16 #define MAX_NAMELEN 36
17 #define MAQMAP_FORMAT_OLD 0
18 #define MAQMAP_FORMAT_NEW -1
20 #define PAIRFLAG_FF 0x01
21 #define PAIRFLAG_FR 0x02
22 #define PAIRFLAG_RF 0x04
23 #define PAIRFLAG_RR 0x08
24 #define PAIRFLAG_PAIRED 0x10
25 #define PAIRFLAG_DIFFCHR 0x20
26 #define PAIRFLAG_NOMATCH 0x40
27 #define PAIRFLAG_SW 0x80
31 uint8_t seq[MAX_READLEN]; /* the last base is the single-end mapping quality. */
32 uint8_t size, map_qual, info1, info2, c[2], flag, alt_qual;
35 char name[MAX_NAMELEN];
42 uint64_t n_mapped_reads;
43 maqmap1_t *mapped_reads;
46 maqmap_t *maq_new_maqmap()
48 maqmap_t *mm = (maqmap_t*)calloc(1, sizeof(maqmap_t));
49 mm->format = MAQMAP_FORMAT_NEW;
52 void maq_delete_maqmap(maqmap_t *mm)
56 for (i = 0; i < mm->n_ref; ++i)
57 free(mm->ref_name[i]);
59 free(mm->mapped_reads);
62 maqmap_t *maqmap_read_header(gzFile fp)
66 mm = maq_new_maqmap();
67 gzread(fp, &mm->format, sizeof(int));
68 if (mm->format != MAQMAP_FORMAT_NEW) {
70 fprintf(stderr, "** Obsolete map format is detected. Please use 'mapass2maq' command to convert the format.\n");
73 assert(mm->format == MAQMAP_FORMAT_NEW);
75 gzread(fp, &mm->n_ref, sizeof(int));
76 mm->ref_name = (char**)calloc(mm->n_ref, sizeof(char*));
77 for (k = 0; k != mm->n_ref; ++k) {
78 gzread(fp, &len, sizeof(int));
79 mm->ref_name[k] = (char*)malloc(len * sizeof(char));
80 gzread(fp, mm->ref_name[k], len);
82 /* read number of mapped reads */
83 gzread(fp, &mm->n_mapped_reads, sizeof(uint64_t));
87 void maq2tam_core(gzFile fp)
93 mm = maqmap_read_header(fp);
94 while ((ret = gzread(fp, m1, sizeof(maqmap1_t))) == sizeof(maqmap1_t)) {
96 if (m1->flag) flag |= 1;
97 if ((m1->flag&PAIRFLAG_PAIRED) || ((m1->flag&PAIRFLAG_SW) && m1->flag != 192)) flag |= 2;
98 if (m1->flag == 192) flag |= 4;
99 if (m1->flag == 64) flag |= 8;
100 if (m1->pos&1) flag |= 0x10;
101 if ((flag&1) && m1->dist != 0) {
104 if (m1->flag&(PAIRFLAG_FF|PAIRFLAG_RF)) c = 0;
105 else if (m1->flag&(PAIRFLAG_FR|PAIRFLAG_RR)) c = 1;
108 if (m1->flag&(PAIRFLAG_FF|PAIRFLAG_FR)) c = 0;
109 else if (m1->flag&(PAIRFLAG_RF|PAIRFLAG_RR)) c = 1;
115 int l = strlen(m1->name);
116 if (m1->name[l-2] == '/') {
117 flag |= (m1->name[l-1] == '1')? 0x40 : 0x80;
118 m1->name[l-2] = '\0';
121 printf("%s\t%d\t", m1->name, flag);
122 printf("%s\t%d\t", mm->ref_name[m1->seqid], (m1->pos>>1)+1);
123 if (m1->flag == 130) {
124 int c = (int8_t)m1->seq[MAX_READLEN-1];
125 printf("%d\t", m1->alt_qual);
126 if (c == 0) printf("%dM\t", m1->size);
128 if (c > 0) printf("%dM%dI%dM\t", m1->map_qual, c, m1->size - m1->map_qual - c);
129 else printf("%dM%dD%dM\t", m1->map_qual, -c, m1->size - m1->map_qual);
132 if (flag&4) printf("0\t*\t");
133 else printf("%d\t%dM\t", m1->map_qual, m1->size);
135 printf("*\t0\t%d\t", m1->dist);
136 for (j = 0; j != m1->size; ++j) {
137 if (m1->seq[j] == 0) putchar('N');
138 else putchar("ACGT"[m1->seq[j]>>6&3]);
141 for (j = 0; j != m1->size; ++j)
142 putchar((m1->seq[j]&0x3f) + 33);
145 printf("MF:i:%d\n", m1->flag);
147 printf("MF:i:%d\t", m1->flag);
148 if (m1->flag) printf("Aq:i:%d\t", m1->alt_qual);
149 printf("NM:i:%d\tUQ:i:%d\tH0:i:%d\tH1:i:%d\n", m1->info1&0xf, m1->info2, m1->c[0], m1->c[1]);
153 fprintf(stderr, "Truncated! Continue anyway.\n");
154 maq_delete_maqmap(mm);
157 int main(int argc, char *argv[])
161 fprintf(stderr, "Usage: maq2tam <in.map>\n");
164 fp = strcmp(argv[1], "-")? gzopen(argv[1], "r") : gzdopen(fileno(stdin), "r");