#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <strings.h>

#include <assert.h>

#include <errno.h>

#include <getopt.h>

#include <unistd.h>

#include <ctype.h>

#include <htslib/vcf.h>

#include <htslib/kstring.h>

#include <htslib/synced_bcf_reader.h>

#include <htslib/kseq.h>

#include <htslib/bgzf.h>

#include "regidx.h"

#include "bcftools.h"

#include "rbuf.h"

#include "filter.h"

#include "smpl_ilist.h"

#define FLT_INCLUDE 1

#define FLT_EXCLUDE 2

#define PICK_REF 1

#define PICK_ALT 2

#define PICK_LONG 4

#define PICK_SHORT 8

#define PICK_IUPAC 16

#define TO_UPPER 1

#define TO_LOWER 2

typedef struct

{

int num; // number of ungapped blocks in this chain

int *block_lengths; // length of the ungapped blocks in this chain

int *ref_gaps; // length of the gaps on the reference sequence between blocks

int *alt_gaps; // length of the gaps on the alternative sequence between blocks

int ori_pos;

int ref_last_block_ori; // start position on the reference sequence of the following ungapped block (0-based)

int alt_last_block_ori; // start position on the alternative sequence of the following ungapped block (0-based)

}

chain_t;

#define MASK_LC 1

#define MASK_UC 2

#define MASK_SKIP(x) (((x)->with!=MASK_LC && (x)->with!=MASK_UC) ? 1 : 0)

typedef struct

{

char *fname, with;

regidx_t *idx;

regitr_t *itr;

}

mask_t;

typedef struct

{

kstring_t fa_buf; // buffered reference sequence

int fa_ori_pos; // start position of the fa_buffer (wrt original sequence)

int fa_frz_pos; // protected position to avoid conflicting variants (last pos for SNPs/ins)

int fa_mod_off; // position difference of fa_frz_pos in the ori and modified sequence (ins positive)

int fa_frz_mod; // the fa_buf offset of the protected fa_frz_pos position, includes the modified sequence

int fa_end_pos; // region's end position in the original sequence

int fa_length; // region's length in the original sequence (in case end_pos not provided in the FASTA header)

int fa_case; // output upper case or lower case: TO_UPPER|TO_LOWER

int fa_src_pos; // last genomic coordinate read from the input fasta (0-based)

char prev_base; // this is only to validate the REF allele in the VCF - the modified fa_buf cannot be used for inserts following deletions, see 600#issuecomment-383186778

int prev_base_pos; // the position of prev_base

int prev_is_insert;

rbuf_t vcf_rbuf;

bcf1_t **vcf_buf;

int nvcf_buf, rid;

char *chr, *chr_prefix;

mask_t *mask;

int nmask;

int chain_id; // chain_id, to provide a unique ID to each chain in the chain output

chain_t *chain; // chain structure to store the sequence of ungapped blocks between the ref and alt sequences

// Note that the chain is re-initialised for each chromosome/seq_region

filter_t *filter;

char *filter_str;

int filter_logic; // include or exclude sites which match the filters? One of FLT_INCLUDE/FLT_EXCLUDE

bcf_srs_t *files;

bcf_hdr_t *hdr;

FILE *fp_out;

FILE *fp_chain;

char **argv;

int argc, output_iupac, iupac_GTs, haplotype, allele, isample, napplied;

uint8_t *iupac_bitmask, *iupac_als;

int miupac_bitmask, miupac_als, regions_overlap;

char *fname, *ref_fname, *sample, *sample_fname, *output_fname, *mask_fname, *chain_fname, missing_allele, absent_allele;

char mark_del, mark_ins, mark_snv;

smpl_ilist_t *smpl;

}

args_t;

static void destroy_chain(chain_t *chain)

{

if ( !chain ) return;

free(chain->ref_gaps);

free(chain->alt_gaps);

free(chain->block_lengths);

free(chain);

}

static chain_t* init_chain(chain_t *chain, int ref_ori_pos)

{

if ( chain ) destroy_chain(chain);

chain = (chain_t*) calloc(1,sizeof(chain_t));

chain->num = 0;

chain->block_lengths = NULL;

chain->ref_gaps = NULL;

chain->alt_gaps = NULL;

chain->ori_pos = ref_ori_pos;

chain->ref_last_block_ori = ref_ori_pos;

chain->alt_last_block_ori = ref_ori_pos;

return chain;

}

static void print_chain(args_t *args)

{

/*

chain_t *chain = args->chain;

int n = chain->num;

int ref_end_pos = args->fa_length + chain->ori_pos;

int last_block_size = ref_end_pos - chain->ref_last_block_ori;

int alt_end_pos = chain->alt_last_block_ori + last_block_size;

int score = 0;

for (n=0; n<chain->num; n++) {

score += chain->block_lengths[n];

}

score += last_block_size;

fprintf(args->fp_chain, "chain %d %s %d + %d %d %s %d + %d %d %d\n", score, args->chr, ref_end_pos, chain->ori_pos, ref_end_pos, args->chr, alt_end_pos, chain->ori_pos, alt_end_pos, ++args->chain_id);

for (n=0; n<chain->num; n++) {

fprintf(args->fp_chain, "%d %d %d\n", chain->block_lengths[n], chain->ref_gaps[n], chain->alt_gaps[n]);

}

fprintf(args->fp_chain, "%d\n\n", last_block_size);

}

static void push_chain_gap(chain_t *chain, int ref_start, int ref_len, int alt_start, int alt_len)

{

// fprintf(stderr, "push_chain_gap(chain=%p, ref_start=%d, ref_len=%d, alt_start=%d, alt_len=%d)\n", chain, ref_start, ref_len, alt_start, alt_len);

int num = chain->num;

if (num && ref_start <= chain->ref_last_block_ori) {

// In case this variant is back-to-back with the previous one

chain->ref_last_block_ori = ref_start + ref_len;

chain->alt_last_block_ori = alt_start + alt_len;

chain->ref_gaps[num-1] += ref_len;

chain->alt_gaps[num-1] += alt_len;

} else {

// Extend the ungapped blocks, store the gap length

chain->block_lengths = (int*) realloc(chain->block_lengths, (num + 1) * sizeof(int));

chain->ref_gaps = (int*) realloc(chain->ref_gaps, (num + 1) * sizeof(int));

chain->alt_gaps = (int*) realloc(chain->alt_gaps, (num + 1) * sizeof(int));

chain->block_lengths[num] = ref_start - chain->ref_last_block_ori;

chain->ref_gaps[num] = ref_len;

chain->alt_gaps[num] = alt_len;

// Update the start positions of the next block

chain->ref_last_block_ori = ref_start + ref_len;

chain->alt_last_block_ori = alt_start + alt_len;

// Increment the number of ungapped blocks

chain->num++;

}

}

static void init_data(args_t *args)

{

args->files = bcf_sr_init();

args->files->require_index = 1;

if ( !bcf_sr_add_reader(args->files,args->fname) ) error("Failed to read from %s: %s\n", !strcmp("-",args->fname)?"standard input":args->fname, bcf_sr_strerror(args->files->errnum));

args->hdr = args->files->readers[0].header;

args->isample = -1;

if ( args->sample_fname )

{

args->smpl = smpl_ilist_init(args->hdr,args->sample_fname,1,SMPL_NONE|SMPL_VERBOSE);

if ( args->smpl && !args->smpl->n ) error("No matching sample found\n");

}

else if ( args->sample && strcmp("-",args->sample) )

{

args->smpl = smpl_ilist_init(args->hdr,args->sample,0,SMPL_NONE|SMPL_VERBOSE);

if ( args->smpl && !args->smpl->n ) error("No matching sample found\n");

}

else if ( !args->sample )

{

args->smpl = smpl_ilist_init(args->hdr,NULL,0,SMPL_NONE|SMPL_VERBOSE);

if ( !args->smpl->n )

{

smpl_ilist_destroy(args->smpl);

args->smpl = NULL;

}

}

if ( args->smpl )

{

if ( args->haplotype || args->allele )

{

if ( args->smpl->n > 1 ) error("Too many samples, only one can be used with -H; check the -s,-S options\n");

args->isample = args->smpl->idx[0];

}

else

{

args->iupac_GTs = 1;

if ( args->smpl->n==1 )

fprintf(stderr,"Note: applying IUPAC codes based on FORMAT/GT in sample %s\n",bcf_hdr_int2id(args->hdr,BCF_DT_SAMPLE,args->smpl->idx[0]));

else

fprintf(stderr,"Note: applying IUPAC codes based on FORMAT/GT in %d samples\n",args->smpl->n);

}

}

else if ( args->output_iupac )

fprintf(stderr,"Note: applying IUPAC codes based on REF,ALT%s\n",bcf_hdr_nsamples(args->hdr)?", ignoring samples":"");

else

fprintf(stderr,"Note: applying REF,ALT variants%s\n",bcf_hdr_nsamples(args->hdr)?", ignoring samples":"");

int i;

for (i=0; i<args->nmask; i++)

{

mask_t *mask = &args->mask[i];

mask->idx = regidx_init(mask->fname,NULL,NULL,0,NULL);

if ( !mask->idx ) error("Failed to initialize mask regions\n");

mask->itr = regitr_init(mask->idx);

}

// In case we want to store the chains

if ( args->chain_fname )

{

args->fp_chain = fopen(args->chain_fname,"w");

if ( ! args->fp_chain ) error("Failed to create %s: %s\n", args->chain_fname, strerror(errno));

args->chain_id = 0;

}

rbuf_init(&args->vcf_rbuf, 100);

args->vcf_buf = (bcf1_t**) calloc(args->vcf_rbuf.m, sizeof(bcf1_t*));

if ( args->output_fname ) {

args->fp_out = fopen(args->output_fname,"w");

if ( ! args->fp_out ) error("Failed to create %s: %s\n", args->output_fname, strerror(errno));

}

else args->fp_out = stdout;

if ( args->filter_str )

args->filter = filter_init(args->hdr, args->filter_str);

args->rid = -1;

}

static void add_mask(args_t *args, char *fname)

{

args->nmask++;

args->mask = (mask_t*)realloc(args->mask,args->nmask*sizeof(*args->mask));

mask_t *mask = &args->mask[args->nmask-1];

mask->fname = fname;

mask->with = 'N';

}

static void add_mask_with(args_t *args, char *with)

{

if ( !args->nmask ) error("The --mask-with option must follow --mask\n");

mask_t *mask = &args->mask[args->nmask-1];

if ( !strcasecmp(with,"uc") ) mask->with = MASK_UC;

else if ( !strcasecmp(with,"lc") ) mask->with = MASK_LC;

else if ( strlen(with)!=1 ) error("Expected \"lc\", \"uc\", or a single character with the --mask-with option\n");

else mask->with = *with;

}

static void destroy_data(args_t *args)

{

free(args->iupac_als);

free(args->iupac_bitmask);

if (args->filter) filter_destroy(args->filter);

if ( args->smpl ) smpl_ilist_destroy(args->smpl);

bcf_sr_destroy(args->files);

int i;

for (i=0; i<args->vcf_rbuf.m; i++)

if ( args->vcf_buf[i] ) bcf_destroy1(args->vcf_buf[i]);

free(args->vcf_buf);

free(args->fa_buf.s);

free(args->chr);

for (i=0; i<args->nmask; i++)

{

mask_t *mask = &args->mask[i];

regidx_destroy(mask->idx);

regitr_destroy(mask->itr);

}

free(args->mask);

if ( args->chain_fname )

if ( fclose(args->fp_chain) ) error("Close failed: %s\n", args->chain_fname);

if ( fclose(args->fp_out) ) error("Close failed: %s\n", args->output_fname);

destroy_chain(args->chain);

}

static void init_region(args_t *args, char *line)

{

char *ss, *se = line;

while ( *se && !isspace_c(*se) && *se!=':' ) se++;

hts_pos_t from = 0, to = 0;

char tmp = 0, *tmp_ptr = NULL;

if ( *se )

{

tmp = *se; *se = 0; tmp_ptr = se;

ss = ++se;

from = strtol(ss,&se,10);

if ( ss==se || !*se || *se!='-' ) from = 0;

else

{

from--;

ss = ++se;

to = strtol(ss,&se,10);

if ( ss==se || (*se && !isspace_c(*se)) ) { from = 0; to = 0; }

else to--;

}

}

free(args->chr);

args->chr = strdup(line);

args->rid = bcf_hdr_name2id(args->hdr,line);

if ( args->rid<0 ) fprintf(stderr,"Warning: Sequence \"%s\" not in %s\n", line,args->fname);

args->prev_base_pos = -1;

args->fa_buf.l = 0;

args->fa_length = 0;

args->fa_end_pos = to; // 0-based

args->fa_ori_pos = from; // 0-based

args->fa_src_pos = from; // 0-based

args->fa_mod_off = 0;

args->fa_frz_pos = -1;

args->fa_frz_mod = -1;

args->fa_case = -1;

args->vcf_rbuf.n = 0;

// bcf_sr_set_regions accepts 1-based coordinates

kstring_t str = {0,0,0};

if ( !from ) from = 1;

else from++;

#ifndef MAX_CSI_COOR

#define MAX_CSI_COOR ((1LL << (14 + 30)) - 1)

if ( to==0 ) to = MAX_CSI_COOR - 1;

else to++;

ksprintf(&str,"%s:%"PRIhts_pos"-%"PRIhts_pos,line,from,to);

bcf_sr_set_opt(args->files,BCF_SR_REGIONS_OVERLAP,args->regions_overlap);

bcf_sr_set_regions(args->files,str.s,0);

free(str.s);

if ( tmp_ptr ) *tmp_ptr = tmp;

fprintf(args->fp_out,">%s%s\n",args->chr_prefix?args->chr_prefix:"",line);

if ( args->chain_fname )

args->chain = init_chain(args->chain, args->fa_ori_pos);

}

static bcf1_t **next_vcf_line(args_t *args)

{

if ( args->vcf_rbuf.n )

{

int i = rbuf_shift(&args->vcf_rbuf);

return &args->vcf_buf[i];

}

while ( bcf_sr_next_line(args->files) )

{

if ( args->filter )

{

int is_ok = filter_test(args->filter, bcf_sr_get_line(args->files,0), NULL);

if ( args->filter_logic & FLT_EXCLUDE ) is_ok = is_ok ? 0 : 1;

if ( !is_ok ) continue;

}

return &args->files->readers[0].buffer[0];

}

if ( args->files->errnum ) error("Error: %s\n", bcf_sr_strerror(args->files->errnum));

return NULL;

}

static void unread_vcf_line(args_t *args, bcf1_t **rec_ptr)

{

bcf1_t *rec = *rec_ptr;

if ( args->vcf_rbuf.n >= args->vcf_rbuf.m )

error("FIXME: too many overlapping records near %s:%"PRId64"\n", bcf_seqname(args->hdr,rec),(int64_t) rec->pos+1);

// Insert the new record in the buffer. The line would be overwritten in

// the next bcf_sr_next_line call, therefore we need to swap it with an

// unused one

int i = rbuf_append(&args->vcf_rbuf);

if ( !args->vcf_buf[i] ) args->vcf_buf[i] = bcf_init1();

bcf1_t *tmp = rec; *rec_ptr = args->vcf_buf[i]; args->vcf_buf[i] = tmp;

}

static void flush_fa_buffer(args_t *args, int len)

{

if ( !args->fa_buf.l ) return;

int nwr = 0;

while ( nwr + 60 <= args->fa_buf.l && (!len || args->fa_ori_pos + nwr + 60 < args->fa_frz_pos) )

{

if ( fwrite(args->fa_buf.s+nwr,1,60,args->fp_out) != 60 ) error("Could not write: %s\n", args->output_fname);

if ( fwrite("\n",1,1,args->fp_out) != 1 ) error("Could not write: %s\n", args->output_fname);

nwr += 60;

}

if ( nwr )

args->fa_ori_pos += nwr;

args->fa_frz_mod -= nwr;

if ( len )

{

// not finished on this chr yet and the buffer cannot be emptied completely

if ( nwr && nwr < args->fa_buf.l )

memmove(args->fa_buf.s,args->fa_buf.s+nwr,args->fa_buf.l-nwr);

args->fa_buf.l -= nwr;

return;

}

// empty the whole buffer

if ( nwr == args->fa_buf.l ) { args->fa_buf.l = 0; return; }

if ( fwrite(args->fa_buf.s+nwr,1,args->fa_buf.l - nwr,args->fp_out) != args->fa_buf.l - nwr ) error("Could not write: %s\n", args->output_fname);

if ( fwrite("\n",1,1,args->fp_out) != 1 ) error("Could not write: %s\n", args->output_fname);

args->fa_ori_pos += args->fa_buf.l - nwr - args->fa_mod_off;

args->fa_mod_off = 0;

args->fa_buf.l = 0;

}

static void apply_absent(args_t *args, hts_pos_t pos)

{

if ( !args->fa_buf.l ) return;

if ( pos <= args->fa_frz_pos + 1 ) return; // if pos==frz+1, then there is no gap, ie=ib, nothing to fill

if ( pos <= args->fa_ori_pos ) return;

int ie = pos && pos - args->fa_ori_pos + args->fa_mod_off < args->fa_buf.l ? pos - args->fa_ori_pos + args->fa_mod_off : args->fa_buf.l;

int ib = args->fa_frz_mod < 0 ? 0 : args->fa_frz_mod;

int i;

for (i=ib; i<ie; i++)

args->fa_buf.s[i] = args->absent_allele;

}

static void freeze_ref(args_t *args, bcf1_t *rec)

{

if ( args->fa_frz_pos >= rec->pos + rec->rlen - 1 ) return;

args->fa_frz_pos = rec->pos + rec->rlen - 1;

args->fa_frz_mod = rec->pos - args->fa_ori_pos + args->fa_mod_off + rec->rlen;

}

static char *mark_del(char *ref, int rlen, char *alt, int mark)

{

char *out = malloc(rlen+1);

int i;

if ( alt )

{

int nalt = strlen(alt);

for (i=0; i<nalt; i++) out[i] = alt[i];

}

else // symbolic <DEL>

{

int nref = strlen(ref);

for (i=0; i<nref; i++) out[i] = ref[i];

}

for (; i<rlen; i++) out[i] = mark;

out[rlen] = 0;

return out;

}

static void mark_ins(char *ref, char *alt, char mark)

{

int i, nref = strlen(ref), nalt = strlen(alt);

if ( mark==TO_LOWER )

for (i=nref; i<nalt; i++) alt[i] = tolower_c(alt[i]);

else if ( mark==TO_UPPER )

for (i=nref; i<nalt; i++) alt[i] = toupper_c(alt[i]);

else if ( mark )

for (i=nref; i<nalt; i++) alt[i] = mark;

}

static void mark_snv(char *ref, char *alt, char mark)

{

int i, nref = strlen(ref), nalt = strlen(alt);

int n = nref < nalt ? nref : nalt;

if ( mark==TO_LOWER )

{

for (i=0; i<n; i++)

if ( tolower_c(ref[i])!=tolower_c(alt[i]) ) alt[i] = tolower_c(alt[i]);

}

else if ( mark==TO_UPPER)

{

for (i=0; i<n; i++)

if ( tolower_c(ref[i])!=tolower_c(alt[i]) ) alt[i] = toupper_c(alt[i]);

}

else if ( mark )

{

for (i=0; i<n; i++)

if ( tolower_c(ref[i])!=tolower_c(alt[i]) ) alt[i] = mark;

}

}

static void iupac_init(args_t *args, bcf1_t *rec)

{

int i;

hts_resize(uint8_t, rec->n_allele, &args->miupac_als, &args->iupac_als, 0);

for (i=0; i<args->miupac_als; i++) args->iupac_als[i] = 0;

}

static int iupac_add_gt(args_t *args, bcf1_t *rec, uint8_t *gt, int ngt)

{

int i, is_set = 0;

for (i=0; i<ngt; i++)

{

if ( bcf_gt_is_missing(gt[i]) ) continue;

if ( gt[i]==(uint8_t)bcf_int8_vector_end ) break;

int ial = bcf_gt_allele(gt[i]);

if ( ial >= rec->n_allele ) error("Invalid VCF, too few ALT alleles at %s:%"PRId64"\n", bcf_seqname(args->hdr,rec),(int64_t) rec->pos+1);

args->iupac_als[ial] = 1;

is_set = 1;

}

return is_set;

}

static int iupac_set_allele(args_t *args, bcf1_t *rec)

{

int i,j, max_len = 0, alt_len = 0, ialt = -1, fallback_alt = -1;

for (i=0; i<rec->n_allele; i++)

{

if ( !args->iupac_als[i] ) continue;

if ( fallback_alt <=0 ) fallback_alt = i;

int l = strlen(rec->d.allele[i]);

for (j=0; j<l; j++)

if ( iupac2bitmask(rec->d.allele[i][j]) < 0 ) break;

if ( j<l ) continue; // symbolic allele, breakpoint or invalid character in the allele

if ( l>max_len )

{

hts_resize(uint8_t, l, &args->miupac_bitmask, &args->iupac_bitmask, HTS_RESIZE_CLEAR);

for (j=max_len; j<l; j++) args->iupac_bitmask[j] = 0;

max_len = l;

}

if ( i>0 && l>alt_len )

{

alt_len = l;

ialt = i;

}

for (j=0; j<l; j++)

args->iupac_bitmask[j] |= iupac2bitmask(rec->d.allele[i][j]);

}

if ( alt_len > 0 )

{

for (j=0; j<alt_len; j++) rec->d.allele[ialt][j] = bitmask2iupac(args->iupac_bitmask[j]);

return ialt;

}

if ( fallback_alt >= 0 ) return fallback_alt;

return ialt;

}

static void apply_variant(args_t *args, bcf1_t *rec)

{

static int warned_haplotype = 0;

if ( args->absent_allele ) apply_absent(args, rec->pos);

if ( rec->n_allele==1 && !args->missing_allele && !args->absent_allele ) { return; }

int i,j;

if ( args->mask )

{

char *chr = (char*)bcf_hdr_id2name(args->hdr,args->rid);

int start = rec->pos;

int end = rec->pos + rec->rlen - 1;

for (i=0; i<args->nmask; i++)

{

mask_t *mask = &args->mask[i];

if ( MASK_SKIP(mask) && regidx_overlap(mask->idx, chr,start,end,NULL) ) return;

}

}

int ialt = 1; // the alternate allele

if ( args->iupac_GTs )

{

bcf_unpack(rec, BCF_UN_FMT);

bcf_fmt_t *fmt = bcf_get_fmt(args->hdr, rec, "GT");

if ( !fmt ) return;

if ( fmt->type!=BCF_BT_INT8 )

error("Todo: GT field represented with BCF_BT_INT8, too many alleles at %s:%"PRId64"?\n",bcf_seqname(args->hdr,rec),(int64_t) rec->pos+1);

ialt = -1;

int is_set = 0;

iupac_init(args,rec);

for (i=0; i<args->smpl->n; i++)

{

uint8_t *ptr = fmt->p + fmt->size*args->smpl->idx[i];

is_set += iupac_add_gt(args, rec, ptr, fmt->n);

}

if ( !is_set && !args->missing_allele ) return;

if ( is_set ) ialt = iupac_set_allele(args, rec);

}

else if ( args->isample >= 0 )

{

bcf_unpack(rec, BCF_UN_FMT);

bcf_fmt_t *fmt = bcf_get_fmt(args->hdr, rec, "GT");

if ( !fmt ) return;

if ( fmt->type!=BCF_BT_INT8 )

error("Todo: GT field represented with BCF_BT_INT8, too many alleles at %s:%"PRId64"?\n",bcf_seqname(args->hdr,rec),(int64_t) rec->pos+1);

uint8_t *ptr = fmt->p + fmt->size*args->isample;

enum { use_hap, use_iupac, pick_one } action = use_hap;

if ( args->allele==PICK_IUPAC )

{

if ( !args->haplotype ) action = use_iupac;

if ( !bcf_gt_is_phased(ptr[0]) && !bcf_gt_is_phased(ptr[fmt->n-1]) ) action = use_iupac;

}

else if ( args->output_iupac ) action = use_iupac;

else if ( !args->haplotype ) action = pick_one;

if ( action==use_hap )

{

if ( args->haplotype > fmt->n )

{

if ( bcf_gt_is_missing(ptr[fmt->n-1]) || bcf_gt_is_missing(ptr[0]) )

{

if ( !args->missing_allele ) return;

ialt = -1;

}

else

{

if ( !warned_haplotype )

{

fprintf(stderr, "Can't apply %d-th haplotype at %s:%"PRId64". (This warning is printed only once.)\n", args->haplotype,bcf_seqname(args->hdr,rec),(int64_t) rec->pos+1);

warned_haplotype = 1;

}

return;

}

}

else

{

ialt = (int8_t)ptr[args->haplotype-1];

if ( bcf_gt_is_missing(ialt) || ialt==bcf_int8_vector_end )

{

if ( !args->missing_allele ) return;

ialt = -1;

}

else

ialt = bcf_gt_allele(ialt);

}

}

else if ( action==use_iupac )

{

ialt = -1;

iupac_init(args,rec);

int is_set = iupac_add_gt(args, rec, ptr, fmt->n);

if ( !is_set && !args->missing_allele ) return;

if ( is_set ) ialt = iupac_set_allele(args, rec);

}

else

{

int is_hom = 1;

for (i=0; i<fmt->n; i++)

{

if ( bcf_gt_is_missing(ptr[i]) )

{

if ( !args->missing_allele ) return; // ignore missing or half-missing genotypes

ialt = -1;

break;

}

if ( ptr[i]==(uint8_t)bcf_int8_vector_end ) break;

ialt = bcf_gt_allele(ptr[i]);

if ( i>0 && ialt!=bcf_gt_allele(ptr[i-1]) ) { is_hom = 0; break; }

}

if ( !is_hom )

{

int prev_len = 0, jalt;

for (i=0; i<fmt->n; i++)

{

if ( ptr[i]==(uint8_t)bcf_int8_vector_end ) break;

jalt = bcf_gt_allele(ptr[i]);

if ( rec->n_allele <= jalt ) error("Broken VCF, too few alts at %s:%"PRId64"\n", bcf_seqname(args->hdr,rec),(int64_t) rec->pos+1);

if ( args->allele & (PICK_LONG|PICK_SHORT) )

{

int len = jalt==0 ? rec->rlen : strlen(rec->d.allele[jalt]);

if ( i==0 ) ialt = jalt, prev_len = len;

else if ( len == prev_len )

{

if ( args->allele & PICK_REF && jalt==0 ) ialt = jalt, prev_len = len;

else if ( args->allele & PICK_ALT && ialt==0 ) ialt = jalt, prev_len = len;

}

else if ( args->allele & PICK_LONG && len > prev_len ) ialt = jalt, prev_len = len;

else if ( args->allele & PICK_SHORT && len < prev_len ) ialt = jalt, prev_len = len;

}

else

{

if ( args->allele & PICK_REF && jalt==0 ) ialt = jalt;

else if ( args->allele & PICK_ALT && ialt==0 ) ialt = jalt;

}

}

}

}

if ( !ialt )

{

// ref allele

if ( args->absent_allele ) freeze_ref(args,rec);

return;

}

if ( rec->n_allele <= ialt ) error("Broken VCF, too few alts at %s:%"PRId64"\n", bcf_seqname(args->hdr,rec),(int64_t) rec->pos+1);

}

else if ( args->output_iupac && rec->n_allele>1 )

{

int ialt, alen = 0, mlen = 0;

for (i=0; i<rec->n_allele; i++)

{

int l = strlen(rec->d.allele[i]);

for (j=0; j<l; j++)

if ( iupac2bitmask(rec->d.allele[i][j]) < 0 ) break;

if ( j<l ) continue; // symbolic allele, breakpoint or invalid character in the allele

if ( l > mlen )

{

hts_expand(uint8_t,l,args->miupac_bitmask,args->iupac_bitmask);

for (j=mlen; j<l; j++) args->iupac_bitmask[j] = 0;

mlen = l;

}

if ( i>0 && l>alen )

{

alen = l;

ialt = i;

}

for (j=0; j<l; j++)

args->iupac_bitmask[j] |= iupac2bitmask(rec->d.allele[i][j]);

}

if ( alen > 0 )

for (j=0; j<alen; j++) rec->d.allele[ialt][j] = bitmask2iupac(args->iupac_bitmask[j]);

else

ialt = 1;

}

if ( rec->n_allele==1 && ialt!=-1 )

{

// non-missing reference

if ( args->absent_allele ) freeze_ref(args,rec);

return;

}

if ( ialt==-1 ) // missing allele, it can be a single position or an entire gvcf block

{

if ( !args->missing_allele ) return;

// Missing alleles cannot overlap an already applied variant.

if ( rec->pos <= args->fa_frz_pos )

{

fprintf(stderr,"The site %s:%"PRId64" overlaps with another variant, skipping...\n", bcf_seqname(args->hdr,rec),(int64_t)rec->pos+1);

return;

}

// Mask the full span of reference blocks and possible deletions.

int set_missing = 0;

if ( rec->rlen>1 )

{

if ( bcf_has_variant_types(rec,VCF_REF,bcf_match_exact)>0 ) set_missing = 1;

else if ( bcf_has_variant_types(rec,VCF_DEL,bcf_match_overlap)>0 ) set_missing = 1;

for (j=1; !set_missing && j<rec->n_allele; j++)

if ( !strncasecmp(rec->d.allele[j],"<DEL",4) ) set_missing = 1;

}

hts_pos_t pos = rec->pos;

hts_pos_t len = set_missing ? rec->rlen : 1;

hts_pos_t idx = pos - args->fa_ori_pos + args->fa_mod_off; // position of the variant within the modified fasta sequence

// Clip to the available FASTA interval before accessing its buffer.

if ( idx<0 )

{

pos -= idx;

len += idx;

idx = 0;

}

if ( len<=0 || idx >= (hts_pos_t)args->fa_buf.l ) return;

if ( len > (hts_pos_t)args->fa_buf.l - idx ) len = (hts_pos_t)args->fa_buf.l - idx;

kstring_t str = {0,0,0};

if ( set_missing )

kputsn(args->fa_buf.s+idx,len,&str);

else

kputc(rec->d.allele[0][0],&str);

kputc(',',&str);

hts_pos_t k;

for (k=0; k<len; k++) kputc(args->missing_allele,&str);

rec->pos = pos;

if ( bcf_update_alleles_str(args->hdr,rec,str.s)<0 )

error("Could not update alleles at %s:%"PRId64"\n", bcf_seqname(args->hdr,rec),(int64_t)rec->pos+1);

free(str.s);

// Use exactly the span represented by the replacement alleles.

rec->rlen = len;

ialt = 1;

}

// An explicit sequence ALT replaces the REF

if ( ialt>0 && rec->d.allele[ialt][0]!='<' ) rec->rlen = strlen(rec->d.allele[0]);

// For some variant types POS+REF refer to the base *before* the event; in such case set trim_beg

int trim_beg = 0;

int var_type = bcf_get_variant_type(rec,ialt);

int var_len = rec->d.var[ialt].n;

if ( var_type & VCF_INDEL )

{

// normally indel starts one base after, but not if the first base of the fa reference is deleted

if ( rec->d.allele[0][0] == rec->d.allele[ialt][0] )

trim_beg = 1;

else

trim_beg = 0;

}

else if ( (var_type & VCF_OTHER) && !strcasecmp(rec->d.allele[ialt],"<DEL>") )

{

trim_beg = 1;

var_len = 1 - rec->rlen;

}

else if ( var_type & VCF_OTHER )

{

if ( !strncasecmp(rec->d.allele[ialt],"<INS",4) ) trim_beg = 1;

//else if ( !strncasecmp(rec->d.allele[ialt],"<DEL",4) ) trim_beg = 1;

}

// Overlapping variant?

if ( ialt==0 && rec->pos <= args->fa_frz_pos && rec->pos + rec->rlen - 1 > args->fa_frz_pos )

{

// Applying the reference allele which overlaps a previous deletion. If we are here, it

// means it goes beyond the freezed position, hence the record can be trimmed and moved

// forward

int ntrim = args->fa_frz_pos - rec->pos + 1;

int nref = strlen(rec->d.allele[0]);

if ( ntrim >= nref )

error("Error: failed to trim overlapping variant at %s:%"PRId64", ntrim=%d >= nref=%d. Is the VCF normalized?\n",

bcf_seqname(args->hdr,rec),(int64_t)rec->pos+1,ntrim,nref);

rec->pos += ntrim;

rec->rlen -= ntrim;

memmove(rec->d.allele[0],rec->d.allele[0]+ntrim,nref-ntrim);

rec->d.allele[0][nref-ntrim] = 0;

}

if ( rec->pos <= args->fa_frz_pos )

{

// Can be still OK iff this is an insertion (and which does not follow another insertion, see #888).

// This still may not be enough for more complicated cases with multiple duplicate positions

// and other types in between. In such case let the user normalize the VCF and remove duplicates.

int overlap = 0;

if ( rec->pos < args->fa_frz_pos || !trim_beg || var_len==0 || args->prev_is_insert ) overlap = 1;

if ( overlap )

{

fprintf(stderr,"The site %s:%"PRId64" overlaps with another variant, skipping...\n", bcf_seqname(args->hdr,rec),(int64_t) rec->pos+1);

return;