samtools/glf.c.pysam.c

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