misc/maq2sam.c

   1 #include <string.h>
   2 #include <zlib.h>
   3 #include <stdio.h>
   4 #include <inttypes.h>
   5 #include <stdlib.h>
   6 #include <assert.h>
   7
   8 //#define MAQ_LONGREADS
   9
  10 #ifdef MAQ_LONGREADS
  11 #  define MAX_READLEN 128
  12 #else
  13 #  define MAX_READLEN 64
  14 #endif
  15
  16 #define MAX_NAMELEN 36
  17 #define MAQMAP_FORMAT_OLD 0
  18 #define MAQMAP_FORMAT_NEW -1
  19
  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
  28
  29 typedef struct
  30 {
  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;
  33         uint32_t seqid, pos;
  34         int dist;
  35         char name[MAX_NAMELEN];
  36 } maqmap1_t;
  37
  38 typedef struct
  39 {
  40         int format, n_ref;
  41         char **ref_name;
  42         uint64_t n_mapped_reads;
  43         maqmap1_t *mapped_reads;
  44 } maqmap_t;
  45
  46 maqmap_t *maq_new_maqmap()
  47 {
  48         maqmap_t *mm = (maqmap_t*)calloc(1, sizeof(maqmap_t));
  49         mm->format = MAQMAP_FORMAT_NEW;
  50         return mm;
  51 }
  52 void maq_delete_maqmap(maqmap_t *mm)
  53 {
  54         int i;
  55         if (mm == 0) return;
  56         for (i = 0; i < mm->n_ref; ++i)
  57                 free(mm->ref_name[i]);
  58         free(mm->ref_name);
  59         free(mm->mapped_reads);
  60         free(mm);
  61 }
  62 maqmap_t *maqmap_read_header(gzFile fp)
  63 {
  64         maqmap_t *mm;
  65         int k, len;
  66         mm = maq_new_maqmap();
  67         gzread(fp, &mm->format, sizeof(int));
  68         if (mm->format != MAQMAP_FORMAT_NEW) {
  69                 if (mm->format > 0) {
  70                         fprintf(stderr, "** Obsolete map format is detected. Please use 'mapass2maq' command to convert the format.\n");
  71                         exit(3);
  72                 }
  73                 assert(mm->format == MAQMAP_FORMAT_NEW);
  74         }
  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);
  81         }
  82         /* read number of mapped reads */
  83         gzread(fp, &mm->n_mapped_reads, sizeof(uint64_t));
  84         return mm;
  85 }
  86
  87 void maq2tam_core(gzFile fp)
  88 {
  89         maqmap_t *mm;
  90         maqmap1_t mm1, *m1;
  91         int ret;
  92         m1 = &mm1;
  93         mm = maqmap_read_header(fp);
  94         while ((ret = gzread(fp, m1, sizeof(maqmap1_t))) == sizeof(maqmap1_t)) {
  95                 int j, flag = 0;
  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) {
 102                         int c;
 103                         if (m1->dist > 0) {
 104                                 if (m1->flag&(PAIRFLAG_FF|PAIRFLAG_RF)) c = 0;
 105                                 else if (m1->flag&(PAIRFLAG_FR|PAIRFLAG_RR)) c = 1;
 106                                 else c = m1->pos&1;
 107                         } else {
 108                                 if (m1->flag&(PAIRFLAG_FF|PAIRFLAG_FR)) c = 0;
 109                                 else if (m1->flag&(PAIRFLAG_RF|PAIRFLAG_RR)) c = 1;
 110                                 else c = m1->pos&1;
 111                         }
 112                         flag |= c;
 113                 }
 114                 if (flag) {
 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';
 119                         }
 120                 }
 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);
 127                         else {
 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);
 130                         }
 131                 } else {
 132                         if (flag&4) printf("0\t*\t");
 133                         else printf("%d\t%dM\t", m1->map_qual, m1->size);
 134                 }
 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]);
 139                 }
 140                 putchar('\t');
 141                 for (j = 0; j != m1->size; ++j)
 142                         putchar((m1->seq[j]&0x3f) + 33);
 143                 putchar('\t');
 144                 if (flag&4) {
 145                         printf("MF:i:%d\n", m1->flag);
 146                 } else {
 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]);
 150                 }
 151         }
 152         if (ret > 0)
 153                 fprintf(stderr, "Truncated! Continue anyway.\n");
 154         maq_delete_maqmap(mm);
 155 }
 156
 157 int main(int argc, char *argv[])
 158 {
 159         gzFile fp;
 160         if (argc == 1) {
 161                 fprintf(stderr, "Usage: maq2tam <in.map>\n");
 162                 return 1;
 163         }
 164         fp = strcmp(argv[1], "-")? gzopen(argv[1], "r") : gzdopen(fileno(stdin), "r");
 165         maq2tam_core(fp);
 166         gzclose(fp);
 167         return 0;
 168 }