NVBIO
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
vcfutils.c
Go to the documentation of this file.
1 #include "htslib/vcfutils.h"
2 
3 int bcf_calc_ac(const bcf_hdr_t *header, bcf1_t *line, int *ac, int which)
4 {
5  int i;
6  for (i=0; i<line->n_allele; i++) ac[i]=0;
7 
8  // Use INFO/AC,AN field only when asked
9  if ( which&BCF_UN_INFO )
10  {
11  bcf_unpack(line, BCF_UN_INFO);
12  int an_id = bcf_hdr_id2int(header, BCF_DT_ID, "AN");
13  int ac_id = bcf_hdr_id2int(header, BCF_DT_ID, "AC");
14  int i, an=-1, ac_len=0, ac_type=0;
15  uint8_t *ac_ptr=NULL;
16  if ( an_id>=0 && ac_id>=0 )
17  {
18  for (i=0; i<line->n_info; i++)
19  {
20  bcf_info_t *z = &line->d.info[i];
21  if ( z->key == an_id ) an = z->v1.i;
22  else if ( z->key == ac_id ) { ac_ptr = z->vptr; ac_len = z->len; ac_type = z->type; }
23  }
24  }
25  if ( an>=0 && ac_ptr )
26  {
27  int nac = 0;
28  #define BRANCH_INT(type_t) { \
29  type_t *p = (type_t *) ac_ptr; \
30  for (i=0; i<ac_len; i++) \
31  { \
32  ac[i+1] = p[i]; \
33  nac += p[i]; \
34  } \
35  }
36  switch (ac_type) {
37  case BCF_BT_INT8: BRANCH_INT(int8_t); break;
38  case BCF_BT_INT16: BRANCH_INT(int16_t); break;
39  case BCF_BT_INT32: BRANCH_INT(int32_t); break;
40  default: fprintf(stderr, "[E::%s] todo: %d at %s:%d\n", __func__, ac_type, header->id[BCF_DT_CTG][line->rid].key, line->pos+1); exit(1); break;
41  }
42  #undef BRANCH_INT
43  assert( an>=nac ); // sanity check for missing values
44  ac[0] = an - nac;
45  return 1;
46  }
47  }
48 
49  // Split genotype fields only when asked
50  if ( which&BCF_UN_FMT )
51  {
52  int i, gt_id = bcf_hdr_id2int(header,BCF_DT_ID,"GT");
53  if ( gt_id<0 ) return 0;
54  bcf_unpack(line, BCF_UN_FMT);
55  bcf_fmt_t *fmt_gt = NULL;
56  for (i=0; i<(int)line->n_fmt; i++)
57  if ( line->d.fmt[i].id==gt_id ) { fmt_gt = &line->d.fmt[i]; break; }
58  if ( !fmt_gt ) return 0;
59  #define BRANCH_INT(type_t,missing,vector_end) { \
60  for (i=0; i<line->n_sample; i++) \
61  { \
62  type_t *p = (type_t*) (fmt_gt->p + i*fmt_gt->size); \
63  int ial; \
64  for (ial=0; ial<fmt_gt->n; ial++) \
65  { \
66  if ( p[ial]==vector_end ) break; /* smaller ploidy */ \
67  if ( !(p[ial]>>1) || p[ial]==missing ) continue; /* missing allele */ \
68  ac[(p[ial]>>1)-1]++; \
69  } \
70  } \
71  }
72  switch (fmt_gt->type) {
73  case BCF_BT_INT8: BRANCH_INT(int8_t, bcf_int8_missing, bcf_int8_vector_end); break;
74  case BCF_BT_INT16: BRANCH_INT(int16_t, bcf_int16_missing, bcf_int16_vector_end); break;
75  case BCF_BT_INT32: BRANCH_INT(int32_t, bcf_int32_missing, bcf_int32_vector_end); break;
76  default: fprintf(stderr, "[E::%s] todo: %d at %s:%d\n", __func__, fmt_gt->type, header->id[BCF_DT_CTG][line->rid].key, line->pos+1); exit(1); break;
77  }
78  #undef BRANCH_INT
79  return 1;
80  }
81  return 0;
82 }
83 
84 int bcf_gt_type(bcf_fmt_t *fmt_ptr, int isample, int *_ial, int *_jal)
85 {
86  int i, nals = 0, has_ref = 0, has_alt = 0, ial = 0, jal = 0;
87  #define BRANCH_INT(type_t,missing,vector_end) { \
88  type_t *p = (type_t*) (fmt_ptr->p + isample*fmt_ptr->size); \
89  for (i=0; i<fmt_ptr->n; i++) \
90  { \
91  if ( p[i] == vector_end ) break; /* smaller ploidy */ \
92  if ( !p[i] || p[i] == missing ) continue; /* missing allele */ \
93  int tmp = p[i]>>1; \
94  if ( tmp>1 ) \
95  { \
96  if ( !ial ) { ial = tmp; has_alt = 1; } \
97  else if ( tmp!=ial ) \
98  { \
99  if ( tmp<ial ) \
100  { \
101  jal = ial; \
102  ial = tmp; \
103  } \
104  else \
105  { \
106  jal = tmp; \
107  } \
108  has_alt = 2; \
109  } \
110  } \
111  else has_ref = 1; \
112  nals++; \
113  } \
114  }
115  switch (fmt_ptr->type) {
116  case BCF_BT_INT8: BRANCH_INT(int8_t, bcf_int8_missing, bcf_int8_vector_end); break;
117  case BCF_BT_INT16: BRANCH_INT(int16_t, bcf_int16_missing, bcf_int16_vector_end); break;
118  case BCF_BT_INT32: BRANCH_INT(int32_t, bcf_int32_missing, bcf_int32_vector_end); break;
119  default: fprintf(stderr, "[E::%s] todo: fmt_type %d\n", __func__, fmt_ptr->type); exit(1); break;
120  }
121  #undef BRANCH_INT
122 
123  if ( _ial ) *_ial = ial>0 ? ial-1 : ial;
124  if ( _jal ) *_jal = jal>0 ? jal-1 : jal;
125  if ( !nals ) return GT_UNKN;
126  if ( nals==1 )
127  return has_ref ? GT_HAPL_R : GT_HAPL_A;
128  if ( !has_ref )
129  return has_alt==1 ? GT_HOM_AA : GT_HET_AA;
130  if ( !has_alt )
131  return GT_HOM_RR;
132  return GT_HET_RA;
133 }
134 
135 int bcf_trim_alleles(const bcf_hdr_t *header, bcf1_t *line)
136 {
137  int i;
138  bcf_fmt_t *gt = bcf_get_fmt(header, line, "GT");
139  if ( !gt ) return 0;
140 
141  int *ac = (int*) calloc(line->n_allele,sizeof(int));
142 
143  // check if all alleles are populated
144  #define BRANCH(type_t,missing,vector_end) { \
145  for (i=0; i<line->n_sample; i++) \
146  { \
147  type_t *p = (type_t*) (gt->p + i*gt->size); \
148  int ial; \
149  for (ial=0; ial<gt->size; ial++) \
150  { \
151  if ( p[ial]==vector_end ) break; /* smaller ploidy */ \
152  if ( !(p[ial]>>1) || p[ial]==missing ) continue; /* missing allele */ \
153  if ( (p[ial]>>1)-1 >= line->n_allele ) return -1; \
154  ac[(p[ial]>>1)-1]++; \
155  } \
156  } \
157  }
158  switch (gt->type) {
159  case BCF_BT_INT8: BRANCH(int8_t, bcf_int8_missing, bcf_int8_vector_end); break;
160  case BCF_BT_INT16: BRANCH(int16_t, bcf_int16_missing, bcf_int16_vector_end); break;
161  case BCF_BT_INT32: BRANCH(int32_t, bcf_int32_missing, bcf_int32_vector_end); break;
162  default: fprintf(stderr, "[E::%s] todo: %d at %s:%d\n", __func__, gt->type, header->id[BCF_DT_CTG][line->rid].key, line->pos+1); exit(1); break;
163  }
164  #undef BRANCH
165 
166  int rm_als = 0, nrm = 0;
167  for (i=1; i<line->n_allele; i++)
168  {
169  if ( !ac[i] ) { rm_als |= 1<<i; nrm++; }
170  }
171  free(ac);
172 
173  if ( nrm ) bcf_remove_alleles(header, line, rm_als);
174  return nrm;
175 }
176 
177 void bcf_remove_alleles(const bcf_hdr_t *header, bcf1_t *line, int rm_mask)
178 {
179  int *map = (int*) calloc(line->n_allele, sizeof(int));
180 
181  // create map of indexes from old to new ALT numbering and modify ALT
182  kstring_t str = {0,0,0};
183  kputs(line->d.allele[0], &str);
184 
185  int nrm = 0, i,j; // i: ori alleles, j: new alleles
186  for (i=1, j=1; i<line->n_allele; i++)
187  {
188  if ( rm_mask & 1<<i )
189  {
190  // remove this allele
191  line->d.allele[i] = NULL;
192  nrm++;
193  continue;
194  }
195  kputc(',', &str);
196  kputs(line->d.allele[i], &str);
197  map[i] = j;
198  j++;
199  }
200  if ( !nrm ) { free(map); free(str.s); return; }
201 
202  int nR_ori = line->n_allele;
203  int nR_new = line->n_allele-nrm;
204  assert(nR_new > 0); // should not be able to remove reference allele
205  int nA_ori = nR_ori-1;
206  int nA_new = nR_new-1;
207 
208  int nG_ori = nR_ori*(nR_ori + 1)/2;
209  int nG_new = nR_new*(nR_new + 1)/2;
210 
211  bcf_update_alleles_str(header, line, str.s);
212 
213  // remove from Number=G, Number=R and Number=A INFO fields.
214  uint8_t *dat = NULL;
215  int mdat = 0, ndat = 0, mdat_bytes = 0, nret;
216  for (i=0; i<line->n_info; i++)
217  {
218  bcf_info_t *info = &line->d.info[i];
219  int vlen = bcf_hdr_id2length(header,BCF_HL_INFO,info->key);
220 
221  if ( vlen!=BCF_VL_A && vlen!=BCF_VL_G && vlen!=BCF_VL_R ) continue; // no need to change
222 
223  int type = bcf_hdr_id2type(header,BCF_HL_INFO,info->key);
224  if ( type==BCF_HT_FLAG ) continue;
225  int size = 1;
226  if ( type==BCF_HT_REAL || type==BCF_HT_INT ) size = 4;
227 
228  mdat = mdat_bytes / size;
229  nret = bcf_get_info_values(header, line, bcf_hdr_int2id(header,BCF_DT_ID,info->key), (void**)&dat, &mdat, type);
230  mdat_bytes = mdat * size;
231  if ( nret<0 )
232  {
233  fprintf(stderr,"[%s:%d %s] Could not access INFO/%s at %s:%d [%d]\n", __FILE__,__LINE__,__FUNCTION__,
234  bcf_hdr_int2id(header,BCF_DT_ID,info->key), bcf_seqname(header,line), line->pos+1, nret);
235  exit(1);
236  }
237  if ( type==BCF_HT_STR )
238  {
239  str.l = 0;
240  char *ss = (char*) dat, *se = (char*) dat;
241  if ( vlen==BCF_VL_A || vlen==BCF_VL_R )
242  {
243  int nexp, inc = 0;
244  if ( vlen==BCF_VL_A )
245  {
246  nexp = nA_ori;
247  inc = 1;
248  }
249  else
250  nexp = nR_ori;
251  for (j=0; j<nexp; j++)
252  {
253  if ( !*se ) break;
254  while ( *se && *se!=',' ) se++;
255  if ( rm_mask & 1<<(j+inc) )
256  {
257  if ( *se ) se++;
258  ss = se;
259  continue;
260  }
261  if ( str.l ) kputc(',',&str);
262  kputsn(ss,se-ss,&str);
263  if ( *se ) se++;
264  ss = se;
265  }
266  assert( j==nexp );
267  }
268  else // Number=G, assuming diploid genotype
269  {
270  int k = 0, n = 0;
271  for (j=0; j<nR_ori; j++)
272  {
273  for (k=0; k<=j; k++)
274  {
275  if ( !*se ) break;
276  while ( *se && *se!=',' ) se++;
277  n++;
278  if ( rm_mask & 1<<j || rm_mask & 1<<k )
279  {
280  if ( *se ) se++;
281  ss = se;
282  continue;
283  }
284  if ( str.l ) kputc(',',&str);
285  kputsn(ss,se-ss,&str);
286  if ( *se ) se++;
287  ss = se;
288  }
289  if ( !*se ) break;
290  }
291  assert( n=nG_ori );
292  }
293 
294  nret = bcf_update_info(header, line, bcf_hdr_int2id(header,BCF_DT_ID,info->key), (void*)str.s, str.l, type);
295  if ( nret<0 )
296  {
297  fprintf(stderr,"[%s:%d %s] Could not update INFO/%s at %s:%d [%d]\n", __FILE__,__LINE__,__FUNCTION__,
298  bcf_hdr_int2id(header,BCF_DT_ID,info->key), bcf_seqname(header,line), line->pos+1, nret);
299  exit(1);
300  }
301  continue;
302  }
303 
304  if ( vlen==BCF_VL_A || vlen==BCF_VL_R )
305  {
306  int inc = 0, ntop;
307  if ( vlen==BCF_VL_A )
308  {
309  assert( nret==nA_ori );
310  ntop = nA_ori;
311  ndat = nA_new;
312  inc = 1;
313  }
314  else
315  {
316  assert( nret==nR_ori );
317  ntop = nR_ori;
318  ndat = nR_new;
319  }
320  int k = 0;
321 
322  #define BRANCH(type_t,is_vector_end) \
323  { \
324  type_t *ptr = (type_t*) dat; \
325  int size = sizeof(type_t); \
326  for (j=0; j<ntop; j++) /* j:ori, k:new */ \
327  { \
328  if ( is_vector_end ) { memcpy(dat+k*size, dat+j*size, size); break; } \
329  if ( rm_mask & 1<<(j+inc) ) continue; \
330  if ( j!=k ) memcpy(dat+k*size, dat+j*size, size); \
331  k++; \
332  } \
333  }
334  switch (type)
335  {
336  case BCF_HT_INT: BRANCH(int32_t,ptr[j]==bcf_int32_vector_end); break;
337  case BCF_HT_REAL: BRANCH(float,bcf_float_is_vector_end(ptr[j])); break;
338  }
339  #undef BRANCH
340  }
341  else // Number=G
342  {
343  assert( nret==nG_ori );
344  int k, l_ori = -1, l_new = 0;
345  ndat = nG_new;
346 
347  #define BRANCH(type_t,is_vector_end) \
348  { \
349  type_t *ptr = (type_t*) dat; \
350  int size = sizeof(type_t); \
351  for (j=0; j<nR_ori; j++) \
352  { \
353  for (k=0; k<=j; k++) \
354  { \
355  l_ori++; \
356  if ( is_vector_end ) { memcpy(dat+l_new*size, dat+l_ori*size, size); break; } \
357  if ( rm_mask & 1<<j || rm_mask & 1<<k ) continue; \
358  if ( l_ori!=l_new ) memcpy(dat+l_new*size, dat+l_ori*size, size); \
359  l_new++; \
360  } \
361  } \
362  }
363  switch (type)
364  {
365  case BCF_HT_INT: BRANCH(int32_t,ptr[l_ori]==bcf_int32_vector_end); break;
366  case BCF_HT_REAL: BRANCH(float,bcf_float_is_vector_end(ptr[l_ori])); break;
367  }
368  #undef BRANCH
369  }
370 
371  nret = bcf_update_info(header, line, bcf_hdr_int2id(header,BCF_DT_ID,info->key), (void*)dat, ndat, type);
372  if ( nret<0 )
373  {
374  fprintf(stderr,"[%s:%d %s] Could not update INFO/%s at %s:%d [%d]\n", __FILE__,__LINE__,__FUNCTION__,
375  bcf_hdr_int2id(header,BCF_DT_ID,info->key), bcf_seqname(header,line), line->pos+1, nret);
376  exit(1);
377  }
378  }
379 
380  // Update GT fields, the allele indexes might have changed
381  for (i=1; i<line->n_allele; i++) if ( map[i]!=i ) break;
382  if ( i<line->n_allele )
383  {
384  mdat = mdat_bytes / 4; // sizeof(int32_t)
385  nret = bcf_get_genotypes(header,line,(void**)&dat,&mdat);
386  mdat_bytes = mdat * 4;
387  if ( nret>0 )
388  {
389  nret /= line->n_sample;
390  int32_t *ptr = (int32_t*) dat;
391  for (i=0; i<line->n_sample; i++)
392  {
393  for (j=0; j<nret; j++)
394  {
395  if ( ptr[j]==bcf_gt_missing ) continue;
396  if ( ptr[j]==bcf_int32_vector_end ) break;
397  int al = bcf_gt_allele(ptr[j]);
398  assert( al<nR_ori && map[al]>=0 );
399  ptr[j] = (map[al]+1)<<1 | (ptr[j]&1);
400  }
401  ptr += nret;
402  }
403  bcf_update_genotypes(header, line, (void*)dat, nret*line->n_sample);
404  }
405  }
406 
407  // Remove from Number=G, Number=R and Number=A FORMAT fields.
408  // Assuming haploid or diploid GTs
409  for (i=0; i<line->n_fmt; i++)
410  {
411  bcf_fmt_t *fmt = &line->d.fmt[i];
412  int vlen = bcf_hdr_id2length(header,BCF_HL_FMT,fmt->id);
413 
414  if ( vlen!=BCF_VL_A && vlen!=BCF_VL_G && vlen!=BCF_VL_R ) continue; // no need to change
415 
416  int type = bcf_hdr_id2type(header,BCF_HL_FMT,fmt->id);
417  if ( type==BCF_HT_FLAG ) continue;
418 
419  int size = 1;
420  if ( type==BCF_HT_REAL || type==BCF_HT_INT ) size = 4;
421 
422  mdat = mdat_bytes / size;
423  nret = bcf_get_format_values(header, line, bcf_hdr_int2id(header,BCF_DT_ID,fmt->id), (void**)&dat, &mdat, type);
424  mdat_bytes = mdat * size;
425  if ( nret<0 )
426  {
427  fprintf(stderr,"[%s:%d %s] Could not access FORMAT/%s at %s:%d [%d]\n", __FILE__,__LINE__,__FUNCTION__,
428  bcf_hdr_int2id(header,BCF_DT_ID,fmt->id), bcf_seqname(header,line), line->pos+1, nret);
429  exit(1);
430  }
431 
432  if ( type==BCF_HT_STR )
433  {
434  int size = nret/line->n_sample; // number of bytes per sample
435  str.l = 0;
436  if ( vlen==BCF_VL_A || vlen==BCF_VL_R )
437  {
438  int nexp, inc = 0;
439  if ( vlen==BCF_VL_A )
440  {
441  nexp = nA_ori;
442  inc = 1;
443  }
444  else
445  nexp = nR_ori;
446  for (j=0; j<line->n_sample; j++)
447  {
448  char *ss = ((char*)dat) + j*size, *se = ss + size, *ptr = ss;
449  int k_src = 0, k_dst = 0, l = str.l;
450  for (k_src=0; k_src<nexp; k_src++)
451  {
452  if ( ptr>=se || !*ptr) break;
453  while ( ptr<se && *ptr && *ptr!=',' ) ptr++;
454  if ( rm_mask & 1<<(k_src+inc) )
455  {
456  ss = ++ptr;
457  continue;
458  }
459  if ( k_dst ) kputc(',',&str);
460  kputsn(ss,ptr-ss,&str);
461  ss = ++ptr;
462  k_dst++;
463  }
464  assert( k_src==nexp );
465  l = str.l - l;
466  for (; l<size; l++) kputc(0, &str);
467  }
468  }
469  else // Number=G, diploid or haploid
470  {
471  for (j=0; j<line->n_sample; j++)
472  {
473  char *ss = ((char*)dat) + j*size, *se = ss + size, *ptr = ss;
474  int k_src = 0, k_dst = 0, l = str.l;
475  int nexp = 0; // diploid or haploid?
476  while ( ptr<se )
477  {
478  if ( !*ptr ) break;
479  if ( *ptr==',' ) nexp++;
480  ptr++;
481  }
482  if ( ptr!=ss ) nexp++;
483  assert( nexp==nG_ori || nexp==nR_ori );
484  ptr = ss;
485  if ( nexp==nG_ori ) // diploid
486  {
487  int ia, ib;
488  for (ia=0; ia<nR_ori; ia++)
489  {
490  for (ib=0; ib<=ia; ib++)
491  {
492  if ( ptr>=se || !*ptr ) break;
493  while ( ptr<se && *ptr && *ptr!=',' ) ptr++;
494  if ( rm_mask & 1<<ia || rm_mask & 1<<ib )
495  {
496  ss = ++ptr;
497  continue;
498  }
499  if ( k_dst ) kputc(',',&str);
500  kputsn(ss,ptr-ss,&str);
501  ss = ++ptr;
502  k_dst++;
503  }
504  if ( ptr>=se || !*ptr ) break;
505  }
506  }
507  else // haploid
508  {
509  for (k_src=0; k_src<nR_ori; k_src++)
510  {
511  if ( ptr>=se || !*ptr ) break;
512  while ( ptr<se && *ptr && *ptr!=',' ) ptr++;
513  if ( rm_mask & 1<<k_src )
514  {
515  ss = ++ptr;
516  continue;
517  }
518  if ( k_dst ) kputc(',',&str);
519  kputsn(ss,ptr-ss,&str);
520  ss = ++ptr;
521  k_dst++;
522  }
523  assert( k_src==nR_ori );
524  l = str.l - l;
525  for (; l<size; l++) kputc(0, &str);
526  }
527  }
528  }
529  nret = bcf_update_format(header, line, bcf_hdr_int2id(header,BCF_DT_ID,fmt->id), (void*)str.s, str.l, type);
530  if ( nret<0 )
531  {
532  fprintf(stderr,"[%s:%d %s] Could not update FORMAT/%s at %s:%d [%d]\n", __FILE__,__LINE__,__FUNCTION__,
533  bcf_hdr_int2id(header,BCF_DT_ID,fmt->id), bcf_seqname(header,line), line->pos+1, nret);
534  exit(1);
535  }
536  continue;
537  }
538 
539  int nori = nret / line->n_sample;
540  if ( vlen==BCF_VL_A || vlen==BCF_VL_R || (vlen==BCF_VL_G && nori==nR_ori) ) // Number=A, R or haploid Number=G
541  {
542  int ntop, inc = 0;
543  if ( vlen==BCF_VL_A )
544  {
545  assert( nori==nA_ori ); // todo: will fail if all values are missing
546  ntop = nA_ori;
547  ndat = nA_new*line->n_sample;
548  inc = 1;
549  }
550  else
551  {
552  assert( nori==nR_ori ); // todo: will fail if all values are missing
553  ntop = nR_ori;
554  ndat = nR_new*line->n_sample;
555  }
556 
557  #define BRANCH(type_t,is_vector_end) \
558  { \
559  for (j=0; j<line->n_sample; j++) \
560  { \
561  type_t *ptr_src = ((type_t*)dat) + j*nori; \
562  type_t *ptr_dst = ((type_t*)dat) + j*nA_new; \
563  int size = sizeof(type_t); \
564  int k_src, k_dst = 0; \
565  for (k_src=0; k_src<ntop; k_src++) \
566  { \
567  if ( is_vector_end ) { memcpy(ptr_dst+k_dst, ptr_src+k_src, size); break; } \
568  if ( rm_mask & 1<<(k_src+inc) ) continue; \
569  if ( k_src!=k_dst ) memcpy(ptr_dst+k_dst, ptr_src+k_src, size); \
570  k_dst++; \
571  } \
572  } \
573  }
574  switch (type)
575  {
576  case BCF_HT_INT: BRANCH(int32_t,ptr_src[k_src]==bcf_int32_vector_end); break;
577  case BCF_HT_REAL: BRANCH(float,bcf_float_is_vector_end(ptr_src[k_src])); break;
578  }
579  #undef BRANCH
580  }
581  else // Number=G, diploid or mixture of haploid+diploid
582  {
583  assert( nori==nG_ori );
584  ndat = nG_new*line->n_sample;
585 
586  #define BRANCH(type_t,is_vector_end) \
587  { \
588  for (j=0; j<line->n_sample; j++) \
589  { \
590  type_t *ptr_src = ((type_t*)dat) + j*nori; \
591  type_t *ptr_dst = ((type_t*)dat) + j*nG_new; \
592  int size = sizeof(type_t); \
593  int ia, ib, k_dst = 0, k_src; \
594  int nset = 0; /* haploid or diploid? */ \
595  for (k_src=0; k_src<nG_ori; k_src++) { if ( is_vector_end ) break; nset++; } \
596  if ( nset==nR_ori ) /* haploid */ \
597  { \
598  for (k_src=0; k_src<nR_ori; k_src++) \
599  { \
600  if ( rm_mask & 1<<k_src ) continue; \
601  if ( k_src!=k_dst ) memcpy(ptr_dst+k_dst, ptr_src+k_src, size); \
602  k_dst++; \
603  } \
604  memcpy(ptr_dst+k_dst, ptr_src+k_src, size); \
605  } \
606  else /* diploid */ \
607  { \
608  k_src = -1; \
609  for (ia=0; ia<nR_ori; ia++) \
610  { \
611  for (ib=0; ib<=ia; ib++) \
612  { \
613  k_src++; \
614  if ( is_vector_end ) { memcpy(ptr_dst+k_dst, ptr_src+k_src, size); ia = nR_ori; break; } \
615  if ( rm_mask & 1<<ia || rm_mask & 1<<ib ) continue; \
616  if ( k_src!=k_dst ) memcpy(ptr_dst+k_dst, ptr_src+k_src, size); \
617  k_dst++; \
618  } \
619  } \
620  } \
621  } \
622  }
623  switch (type)
624  {
625  case BCF_HT_INT: BRANCH(int32_t,ptr_src[k_src]==bcf_int32_vector_end); break;
626  case BCF_HT_REAL: BRANCH(float,bcf_float_is_vector_end(ptr_src[k_src])); break;
627  }
628  #undef BRANCH
629  }
630  nret = bcf_update_format(header, line, bcf_hdr_int2id(header,BCF_DT_ID,fmt->id), (void*)dat, ndat, type);
631  if ( nret<0 )
632  {
633  fprintf(stderr,"[%s:%d %s] Could not update FORMAT/%s at %s:%d [%d]\n", __FILE__,__LINE__,__FUNCTION__,
634  bcf_hdr_int2id(header,BCF_DT_ID,fmt->id), bcf_seqname(header,line), line->pos+1, nret);
635  exit(1);
636  }
637  }
638  free(dat);
639  free(str.s);
640  free(map);
641 }
642 
Generated on Wed Feb 25 2015 08:32:49 for NVBIO by   doxygen 1.8.4