glf.c

   1 #include <string.h>
   2 #include <stdlib.h>
   3 #include "glf.h"
   4
   5 #ifdef _NO_BGZF
   6 // then alias bgzf_*() functions
   7 #endif
   8
   9 static int glf3_is_BE = 0;
  10
  11 static inline uint32_t bam_swap_endian_4(uint32_t v)
  12 {
  13         v = ((v & 0x0000FFFFU) << 16) | (v >> 16);
  14         return ((v & 0x00FF00FFU) << 8) | ((v & 0xFF00FF00U) >> 8);
  15 }
  16
  17 static inline uint16_t bam_swap_endian_2(uint16_t v)
  18 {
  19         return (uint16_t)(((v & 0x00FF00FFU) << 8) | ((v & 0xFF00FF00U) >> 8));
  20 }
  21
  22 static inline int bam_is_big_endian()
  23 {
  24         long one= 1;
  25         return !(*((char *)(&one)));
  26 }
  27
  28 glf3_header_t *glf3_header_init()
  29 {
  30         glf3_is_BE = bam_is_big_endian();
  31         return (glf3_header_t*)calloc(1, sizeof(glf3_header_t));
  32 }
  33
  34 glf3_header_t *glf3_header_read(glfFile fp)
  35 {
  36         glf3_header_t *h;
  37         char magic[4];
  38         h = glf3_header_init();
  39         bgzf_read(fp, magic, 4);
  40         if (strncmp(magic, "GLF\3", 4)) {
  41                 fprintf(stderr, "[glf3_header_read] invalid magic. Abort.\n");
  42                 exit(1);
  43         }
  44         bgzf_read(fp, &h->l_text, 4);
  45         if (glf3_is_BE) h->l_text = bam_swap_endian_4(h->l_text);
  46         if (h->l_text) {
  47                 h->text = (uint8_t*)calloc(h->l_text + 1, 1);
  48                 bgzf_read(fp, h->text, h->l_text);
  49         }
  50         return h;
  51 }
  52
  53 void glf3_header_write(glfFile fp, const glf3_header_t *h)
  54 {
  55         int32_t x;
  56         bgzf_write(fp, "GLF\3", 4);
  57         x = glf3_is_BE? bam_swap_endian_4(h->l_text) : h->l_text;
  58         bgzf_write(fp, &x, 4);
  59         if (h->l_text) bgzf_write(fp, h->text, h->l_text);
  60 }
  61
  62 void glf3_header_destroy(glf3_header_t *h)
  63 {
  64         free(h->text);
  65         free(h);
  66 }
  67
  68 char *glf3_ref_read(glfFile fp, int *len)
  69 {
  70         int32_t n, x;
  71         char *str;
  72         *len = 0;
  73         if (bgzf_read(fp, &n, 4) != 4) return 0;
  74         if (glf3_is_BE) n = bam_swap_endian_4(n);
  75         if (n < 0) {
  76                 fprintf(stderr, "[glf3_ref_read] invalid reference name length: %d.\n", n);
  77                 return 0;
  78         }
  79         str = (char*)calloc(n + 1, 1); // not necesarily n+1 in fact
  80         x = bgzf_read(fp, str, n);
  81         x += bgzf_read(fp, len, 4);
  82         if (x != n + 4) {
  83                 free(str); *len = -1; return 0; // truncated
  84         }
  85         if (glf3_is_BE) *len = bam_swap_endian_4(*len);
  86         return str;
  87 }
  88
  89 void glf3_ref_write(glfFile fp, const char *str, int len)
  90 {
  91         int32_t m, n = strlen(str) + 1;
  92         m = glf3_is_BE? bam_swap_endian_4(n) : n;
  93         bgzf_write(fp, &m, 4);
  94         bgzf_write(fp, str, n);
  95         if (glf3_is_BE) len = bam_swap_endian_4(len);
  96         bgzf_write(fp, &len, 4);
  97 }
  98
  99 void glf3_view1(const char *ref_name, const glf3_t *g3, int pos)
 100 {
 101         int j;
 102         if (g3->rtype == GLF3_RTYPE_END) return;
 103         printf("%s\t%d\t%c\t%d\t%d\t%d", ref_name, pos + 1,
 104                    g3->rtype == GLF3_RTYPE_INDEL? '*' : "XACMGRSVTWYHKDBN"[g3->ref_base],
 105                    g3->depth, g3->rms_mapQ, g3->min_lk);
 106         if (g3->rtype == GLF3_RTYPE_SUB)
 107                 for (j = 0; j != 10; ++j) printf("\t%d", g3->lk[j]);
 108         else {
 109                 printf("\t%d\t%d\t%d\t%d\t%d\t%s\t%s\t", g3->lk[0], g3->lk[1], g3->lk[2], g3->indel_len[0], g3->indel_len[1],
 110                            g3->indel_len[0]? g3->indel_seq[0] : "*", g3->indel_len[1]? g3->indel_seq[1] : "*");
 111         }
 112         printf("\n");
 113 }
 114
 115 int glf3_write1(glfFile fp, const glf3_t *g3)
 116 {
 117         int r;
 118         uint8_t c;
 119         uint32_t y[2];
 120         c = g3->rtype<<4 | g3->ref_base;
 121         r = bgzf_write(fp, &c, 1);
 122         if (g3->rtype == GLF3_RTYPE_END) return r;
 123         y[0] = g3->offset;
 124         y[1] = g3->min_lk<<24 | g3->depth;
 125         if (glf3_is_BE) {
 126                 y[0] = bam_swap_endian_4(y[0]);
 127                 y[1] = bam_swap_endian_4(y[1]);
 128         }
 129         r += bgzf_write(fp, y, 8);
 130         r += bgzf_write(fp, &g3->rms_mapQ, 1);
 131         if (g3->rtype == GLF3_RTYPE_SUB) r += bgzf_write(fp, g3->lk, 10);
 132         else {
 133                 int16_t x[2];
 134                 r += bgzf_write(fp, g3->lk, 3);
 135                 x[0] = glf3_is_BE? bam_swap_endian_2(g3->indel_len[0]) : g3->indel_len[0];
 136                 x[1] = glf3_is_BE? bam_swap_endian_2(g3->indel_len[1]) : g3->indel_len[1];
 137                 r += bgzf_write(fp, x, 4);
 138                 if (g3->indel_len[0]) r += bgzf_write(fp, g3->indel_seq[0], abs(g3->indel_len[0]));
 139                 if (g3->indel_len[1]) r += bgzf_write(fp, g3->indel_seq[1], abs(g3->indel_len[1]));
 140         }
 141         return r;
 142 }
 143
 144 #ifndef kv_roundup32
 145 #define kv_roundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
 146 #endif
 147
 148 int glf3_read1(glfFile fp, glf3_t *g3)
 149 {
 150         int r;
 151         uint8_t c;
 152         uint32_t y[2];
 153         r = bgzf_read(fp, &c, 1);
 154         if (r == 0) return 0;
 155         g3->ref_base = c & 0xf;
 156         g3->rtype = c>>4;
 157         if (g3->rtype == GLF3_RTYPE_END) return r;
 158         r += bgzf_read(fp, y, 8);
 159         if (glf3_is_BE) {
 160                 y[0] = bam_swap_endian_4(y[0]);
 161                 y[1] = bam_swap_endian_4(y[1]);
 162         }
 163         g3->offset = y[0];
 164         g3->min_lk = y[1]>>24;
 165         g3->depth = y[1]<<8>>8;
 166         r += bgzf_read(fp, &g3->rms_mapQ, 1);
 167         if (g3->rtype == GLF3_RTYPE_SUB) r += bgzf_read(fp, g3->lk, 10);
 168         else {
 169                 int16_t x[2], max;
 170                 r += bgzf_read(fp, g3->lk, 3);
 171                 r += bgzf_read(fp, x, 4);
 172                 if (glf3_is_BE) {
 173                         x[0] = bam_swap_endian_2(x[0]);
 174                         x[1] = bam_swap_endian_2(x[1]);
 175                 }
 176                 g3->indel_len[0] = x[0];
 177                 g3->indel_len[1] = x[1];
 178                 x[0] = abs(x[0]); x[1] = abs(x[1]);
 179                 max = (x[0] > x[1]? x[0] : x[1]) + 1;
 180                 if (g3->max_len < max) {
 181                         g3->max_len = max;
 182                         kv_roundup32(g3->max_len);
 183                         g3->indel_seq[0] = (char*)realloc(g3->indel_seq[0], g3->max_len);
 184                         g3->indel_seq[1] = (char*)realloc(g3->indel_seq[1], g3->max_len);
 185                 }
 186                 r += bgzf_read(fp, g3->indel_seq[0], x[0]);
 187                 r += bgzf_read(fp, g3->indel_seq[1], x[1]);
 188                 g3->indel_seq[0][x[0]] = g3->indel_seq[1][x[1]] = 0;
 189         }
 190         return r;
 191 }
 192
 193 void glf3_view(glfFile fp)
 194 {
 195         glf3_header_t *h;
 196         char *name;
 197         glf3_t *g3;
 198         int len;
 199         h = glf3_header_read(fp);
 200         g3 = glf3_init1();
 201         while ((name = glf3_ref_read(fp, &len)) != 0) {
 202                 int pos = 0;
 203                 while (glf3_read1(fp, g3) && g3->rtype != GLF3_RTYPE_END) {
 204                         pos += g3->offset;
 205                         glf3_view1(name, g3, pos);
 206                 }
 207                 free(name);
 208         }
 209         glf3_header_destroy(h);
 210         glf3_destroy1(g3);
 211 }
 212
 213 int glf3_view_main(int argc, char *argv[])
 214 {
 215         glfFile fp;
 216         if (argc == 1) {
 217                 fprintf(stderr, "Usage: glfview <in.glf>\n");
 218                 return 1;
 219         }
 220         fp = (strcmp(argv[1], "-") == 0)? bgzf_fdopen(fileno(stdin), "r") : bgzf_open(argv[1], "r");
 221         if (fp == 0) {
 222                 fprintf(stderr, "Fail to open file '%s'\n", argv[1]);
 223                 return 1;
 224         }
 225         glf3_view(fp);
 226         bgzf_close(fp);
 227         return 0;
 228 }
 229
 230 #ifdef GLFVIEW_MAIN
 231 int main(int argc, char *argv[])
 232 {
 233         return glf3_view_main(argc, argv);
 234 }
 235 #endif