mvcf2consensus.py

#!/usr/bin/python
# -*- coding: UTF-8 -*-

'''
Author: Krisian Ullrich
date: January 2017
email: ullrich@evolbio.mpg.de
License: MIT

The MIT License (MIT)

Copyright (c) 2017 Kristian Ullrich

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

This software relies on the following python packages:
sys
os
re
numpy https://pypi.python.org/pypi/numpy
random
argparse
tabix https://pypi.python.org/pypi/pytabix
gzip 
itertools
Bio https://github.com/biopython/DIST
collections
'''


import sys
import os
import re
import numpy as np
import random
import argparse
import tabix
import gzip
import itertools
from itertools import repeat
from Bio import SeqIO
import collections


IUPACdict = {'A':'A','C':'C','G':'G','T':'T','-':'-','N':'N','.':'.','AG':'R','GA':'R','CT':'Y','TC':'Y','CG':'S',
'GC':'S','AT':'W', 'TA': 'W', 'GT': 'K', 'TG': 'K', 'AC': 'M', 'CA': 'M','CGT':'B','CTG':'B','GCT':'B','GTC':'B',
'TCG':'B','TGC':'B','AGT':'D','ATG':'D','GAT':'D','GTA':'D','TAG':'D','TGA':'D','ACT':'H','ATC':'H','CAT':'H',
'CTA':'H','TAC':'H','TCA':'H','ACG':'V','AGC':'V','CAG':'V','CGA':'V','GAC':'V','GCA':'V','ACGT':'N','ACTG':'N',
'AGCT':'N','AGTC':'N','ATCG':'N','ATGC':'N','CAGT':'N','CATG':'N','CGAT':'N','CGTA':'N','CTAG':'N','CTGA':'N',
'GACT':'N','GATC':'N','GCAT':'N','GCTA':'N','GTAC':'N','GTCA':'N','TACG':'N','TAGC':'N','TCAG':'N','TCGA':'N',
'TGAC':'N','TGCA':'N','A*':'N','C*':'N','G*':'N','T*':'N'}


def getIUPAC(x):
    return IUPACdict[x]


def POINT2zero(x):
    if x == '.':
        return '0'
    else:
        return x


def POINTNA2zero(x):
    if x == '.':
        return '0'
    if x == 'NA':
        return '0'
    else:
        return x


def NA2zero(x):
    if x == 'NA':
        return int(0)
    else:
        return int(x)


def nCk(n, k):
    return len(list(itertools.combinations(range(n), k)))


def evaluate_ref_alt_max(x,y):
    if x[0]==0 and x[1]==0:
        return './.'
    if x[0]>0 and x[1]==0:
        return str(y[0])+'/'+str(y[0])
    if x[0]==0 and x[1]>0:
        return str(y[1])+'/'+str(y[1])
    if x[0]>0 and x[1]>0:
        return str(y[0])+'/'+str(y[1])


def get_header_from_vcf( argsDICT ):
    HEADER = []
    with open(argsDICT['ivcf'],'rb') as f:
        for flines in f:
           if flines[:2]=='##':
               HEADER.append(flines)
           if flines[:4]=='#CHR':
               HEADER.append('##INFO=<ID=NCFnumber,Number=1,Type=String,Description="NotCalledFraction number">\n')
               HEADER.append('##INFO=<ID=NCFfraction,Number=1,Type=Float,Description="NotCalledFraction fraction">\n')
               HEADER.append('##mvcf2consensus.py '+argsDICT['args'][0]+'\n')
               HEADER.append(changeHEADER(flines, argsDICT['add'], argsDICT['id'], argsDICT['samples']))
           if flines[0]!='#':
               break
    return HEADER


def get_header_from_gz( argsDICT ):
    HEADER = []
    with gzip.open(argsDICT['ivcf'],'rb') as f:
        for flines in f:
           if flines[:2]=='##':
               HEADER.append(flines)
           if flines[:4]=='#CHR':
               HEADER.append('##INFO=<ID=NCFnumber,Number=1,Type=String,Description="NotCalledFraction number">\n')
               HEADER.append('##INFO=<ID=NCFfraction,Number=1,Type=Float,Description="NotCalledFraction fraction">\n')
               HEADER.append('##mvcf2consensus.py '+argsDICT['args'][0]+'\n')
               HEADER.append(changeHEADER(flines, argsDICT['add'], argsDICT['id'], argsDICT['samples']))
           if flines[0]!='#':
               break
    return HEADER


def get_sample_pos_from_header( argsDICT ):
    SAMPLE_POS = []
    if argsDICT['tabix']:
        with gzip.open(argsDICT['ivcf'],'rb') as f:
            for flines in f:
               if flines[:4]=='#CHR':
                   flines_stripped = flines.strip().split('\t')
               if flines[0]!='#':
                   break
    if not argsDICT['tabix']:
        with open(argsDICT['ivcf'],'rb') as f:
            for flines in f:
                if flines[:4]=='#CHR':
                    flines_stripped = flines.strip().split('\t')
                if flines[0]!='#':
                    break
    for s in argsDICT['samples']:
        SAMPLE_POS.append([x for x,y in enumerate(flines_stripped) if y == s][0])
    return SAMPLE_POS


def get_sample_pos_dict_from_header( argsDICT ):
    sampleposDICT = {}
    if argsDICT['tabix']:
        with gzip.open(argsDICT['ivcf'],'rb') as f:
            for flines in f:
               if flines[:4]=='#CHR':
                   flines_stripped = flines.strip().split('\t')
               if flines[0]!='#':
                   break
    if not argsDICT['tabix']:
        with open(argsDICT['ivcf'],'rb') as f:
            for flines in f:
                if flines[:4]=='#CHR':
                    flines_stripped = flines.strip().split('\t')
                if flines[0]!='#':
                    break
    for s in argsDICT['samples']:
        sampleposDICT[s]=[x for x,y in enumerate(flines_stripped) if y == s][0]
    return sampleposDICT



def changeHEADER(header, add, name, samples):
    if add:
        headersplit = header.strip().split('\tFORMAT')
        return('\t'.join([headersplit[0]]+['FORMAT']+samples+[name])+'\n')
    if not add:
        headersplit = header.strip().split('\tFORMAT')
        return('\t'.join([headersplit[0]]+['FORMAT']+[name])+'\n')


def get_missing_type(x):
    if x == 'keep':
        return ['keep', 0, 0]
    if x == 'zero':
        return ['zero', 0, 0]
    if 'set' in x:
        return ['set', int(x.split(':')[1]), int(x.split(':')[2])]


def get_reference_len( argsDICT ):
    print 'start obtaining reference length'
    length_dict = {}
    for seq in SeqIO.parse( argsDICT['R'], 'fasta' ):
        if argsDICT['verbose']:
            print seq.id
        if seq.id not in argsDICT['chr']:
            continue
        if seq.id in argsDICT['chr']:
            length_dict[seq.id] = len( seq )
    print 'finished obtaining reference length'    
    return length_dict


class vcfRecord:
    
    
    def __init__(self, record, CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, FORMAT, SAMPLES, tabix):
        if not tabix:
            self.splitted = record.strip().split('\t')
        if tabix:
            self.splitted = record
        self.CHROM = self.splitted[CHROM]
        self.POS = self.splitted[POS]
        self.ID = self.splitted[ID]
        self.REF = [self.splitted[REF]]
        self.ALT = self.splitted[ALT].split(',')
        self.QUAL = self.splitted[QUAL]
        self.FILTER = self.splitted[FILTER]
        self.INFO = self.splitted[INFO]
        self.FORMAT = self.splitted[FORMAT].split(':')
        self.SAMPLES = collections.OrderedDict()
        for s in SAMPLES:
            self.SAMPLES[s] = collections.OrderedDict()
            if len(self.splitted[s].split(':')) != len(self.FORMAT):
                self.splitted[s]+=':'+':'.join(list(np.repeat("0",len(self.FORMAT)-len(self.splitted[s].split(':')))))
            for p,f in enumerate(self.FORMAT):
                self.SAMPLES[s][f] = self.splitted[s].split(':')[p].split(',')
    
    
    def get_var_len(self):
        """" Return length of alleles """
        var_len = [len(self.REF)]
        var_len += [len(x) for x in self.ALT]
        return var_len
    
    
    def get_var_number(self):
        """ Return total number of alleles """
        return len([self.REF])+len(self.ALT)
    
    
    def get_ref_alt(self):
        """ Return alleles """
        REF_ALT = self.REF+self.ALT
        return REF_ALT
    
    
    def is_filtered(self, argsDICT):
        """ Return True if FILTER of record is in filterarg """
        if len(argsDICT['filterarg'])!=1:
            if self.FILTER not in argsDICT['filterarg']:
                return True
            if self.FILTER in argsDICT['filterarg']:
                return False
        if len(argsDICT['filterarg']) == 1 and argsDICT['filterarg'][0]!='':
            if self.FILTER not in argsDICT['filterarg']:
                return True
            if self.FILTER in argsDICT['filterarg']:
                return False
        if len(argsDICT['filterarg']) == 1 and argsDICT['filterarg'][0] == '':
            return False
    
    
    def get_var_type(self):
        """ Return variant type according to GATK """
        var_type = 'NA'
        var_len = self.get_var_len()
        var_number = self.get_var_number()
        if len(self.ALT) == 1 and self.ALT[0] == '.':
            var_type = 'noVAR'
        if len(self.ALT) == 1 and self.ALT[0] == '*':
            var_type = 'complexVAR'
        if var_number == 2:
            if var_len[0] == 1 and var_len[1] == 1:
                #equals GATK SNP
                var_type = 'biSNP'
            if var_len[0] > 1 and var_len[1] > 1 and var_len[0] == var_len[1]:
                #equals GATK MNP
                var_type = 'repSNP'
            if var_len[0] > var_len[1]:
                var_type = 'deletion'
            if var_len[0] < var_len[1]:
                var_type = 'insertion'
        #equals GATK MULTI-ALLELIC
        if var_number != 2:
            if all( [x == var_len[0] for x in var_len[1::]] ) and all( [x == 1 for x in var_len] ):
                #equals GATK SNP
                var_type = 'muSNP'
            #TODO 'murepSNP'
            if all( [x < var_len[0] for x in var_len[1::]] ):
                var_type = 'mudeletion'
            if all( [x > var_len[0] for x in var_len[1::]] ):
                var_type = 'muinsertion'
            if any( [x < var_len[0] for x in var_len[1::]] ) and any( [x == var_len[0] for x in var_len[1::]] ):
                var_type = 'complexDS'
            if any( [x > var_len[0] for x in var_len[1::]] ) and any( [x == var_len[0] for x in var_len[1::]] ):
                var_type = 'complexIS'
            if any( [x < var_len[0] for x in var_len[1::]] ) and any( [x > var_len[0] for x in var_len[1::]] ):
                var_type = 'complexDI'
            if any( [x < var_len[0] for x in var_len[1::]] ) and any( [x > var_len[0] for x in var_len[1::]] ) and any( [x == var_len[0] for x in var_len[1::]] ):
                var_type = 'complexDIS'
        return var_type
    
    
    def keep_samples(self, argsDICT):
        """ Keep only samples in record which are listed """
        for s in self.SAMPLES.copy():
            if s not in argsDICT['samples_pos']:
                self.SAMPLES.pop(s)
    
    
    def get_info(self, value):
        """" Return INFO value """
        return [x[1] for x in [x.split('=') for x in self.INFO.split(';')] if x[0] == value][0]
    
    
    def add_info(self, value):
        """ Add INFO value """
        self.INFO+=str(value)
    
    
    def add_format_on_samples(self, s, f, v):
        """ Add new FORMAT field to sample """
        self.SAMPLES[s][f] = v
    
    
    def replace_format_on_samples(self, s, f, v):
        """ Replace FORMAT field of sample """ 
        self.SAMPLES[s][f] = v
    
    
    def get_samples_ad(self):
        """ Return samples AD as list """
        return [self.SAMPLES[x]['AD'] for x in self.SAMPLES]
    
    
    def get_samples_gt(self):
        """ Return samples GT as list """
        return [self.SAMPLES[x]['GT'] for x in self.SAMPLES]
    
    
    def get_samples_gq(self):
        """ Return samples GQ as list """
        return [self.SAMPLES[x]['GQ'] for x in self.SAMPLES]
    
    
    def get_samples_dp(self):
        """ Return samples GQ as list """
        return [self.SAMPLES[x]['DP'] for x in self.SAMPLES]
    
    
    def get_samples_pl(self):
        """ Return samples GQ as list """
        return [self.SAMPLES[x]['PL'] for x in self.SAMPLES]
    
    
    def get_samples_ad_len(self):
        """ Return samples AD len as list """
        return [len(self.SAMPLES[x]['AD']) for x in self.SAMPLES]
    
    
    def get_samples_ad_ref(self):
        """ Return samples AD REF as list """
        return [self.SAMPLES[x]['AD'][0] for x in self.SAMPLES]
    
    
    def get_samples_ad_alt(self):
        """ Return samples AD ALT as list """
        return [self.SAMPLES[x]['AD'][1:] for x in self.SAMPLES]
    
    
    def gt_type(self, x):
        if len(x)==1:
            if x[0] == '.':
                return 'NA'
            if x[0] != '.':
                return 'hap'
        if x[0] == '.' and x[1] == '.':
            return 'NA'
        if x[0] != '.' and x[1] != '.' and x[0] == x[1]:
            return 'hom'
        if x[0] != '.' and x[1] != '.' and x[0] != x[1]:
            return 'het'
    
    
    def get_gt_type(self):
        """ Return GT type """
        gt_split = re.compile("[|/]")
        GTs = [re.split(gt_split,x[0]) for x in self.get_samples_gt()]
        return [self.gt_type(x) for x in GTs]
    
    
    def get_gt_na_pos(self):
        """ Return index of GT which is NA """
        samplesgttype = self.get_gt_type()
        samplesgtnapos = [x for x,y in enumerate(samplesgttype) if y == 'NA']
        return samplesgtnapos
    
    
    def fill_empty_ad(self):
        """ Replaces "." with 0 and fills to var_number if less """
        x_ = [x+list(np.repeat(0,self.get_var_number()-len(x))) for x in [[int(x.replace(".","0")) for x in y] for y in self.get_samples_ad()]]
        for x,y in zip(x_, self.SAMPLES.keys()):
            self.SAMPLES[y]['AD'] = x
     
    
    def fill_empty_pl(self):
        """ Replaces "." with [0,0,0] """
        x_ = [x+list(np.repeat(0,nCk(self.get_var_number()+1,2)-len(x))) for x in [[int(x.replace(".","0")) for x in y] for y in self.get_samples_pl()]]
        for x,y in zip(x_, self.SAMPLES.keys()):
            self.SAMPLES[y]['PL'] = x
    
    
    def set_samples_ad_by_gt(self, argsDICT):
        if argsDICT['missingGT'][0] == 'keep':
            return
        if argsDICT['missingGT'][0] == 'zero':
            samplesgtnapos = slef.get_gt_na_pos()
            #set AD to zero for all samples with GT equals NA
            for s in [x.sample for y,x in enumerate(self.samples) if y in samplesgtnapos]:
                replace_format_on_samples(self, s, 'AD', list(np.repeat(0, get_var_number(self))))
        if argsDICT['missingGT'][0] == 'set':
            samplesgtnapos = self.get_gt_na_pos()
            #set AD to zero for all samples with GT equals NA
            for s in self.get_gt_na_pos():
                self.replace_format_on_samples(self.SAMPLES.keys()[s], 'AD', [argsDICT['missingGT'][1]]+list(np.repeat(argsDICT['missingGT'][2], self.get_var_number()-1)))
    
    
    def set_dp_by_ad(self):
        samplesad = self.get_samples_ad()
        samplesadsum = [np.sum([int(x) for x in x]) for x in samplesad]
        #set DP to ADSUM
        for i in range(0, len(self.SAMPLES)):
            self.replace_format_on_samples(self.SAMPLES.keys()[i], 'DP', samplesadsum[i])
    
    
    def set_samples_gt_by_ad(self, argsDICT):
        if argsDICT['resetGT'] == 'keep':
            return
        if argsDICT['resetGT'] == 'AD':
            #setGTbyAD
            samplesad_ref = self.get_samples_ad_ref()
            samplesad_alt = self.get_samples_ad_alt()
            #NOTE np.argmax return the first element for equal values which would bias GT setting always to the first ALT entry if equal depth exist
            #NOTE this is only used to determine NotCalledFraction, AD for other alleles will be considered
            samplesad_alt_maxpos = [np.argmax(x) for x in samplesad_alt]
            samplesad_alt_max = [np.max(x) for x in samplesad_alt]
            samplesad_ref_alt_max = map(lambda a, b: [a,b], samplesad_ref, samplesad_alt_max)
            samplesad_ref_alt_maxpos = map(lambda a, b: [a,b], list(np.repeat(0,len(samplesad_ref))), [x+1 for x in samplesad_alt_maxpos])
            samplesgt = []
            for a, b in zip(samplesad_ref_alt_max, samplesad_ref_alt_maxpos):
                samplesgt.append(evaluate_ref_alt_max(a, b))
            #get NotCalled positions
            minDP_NCpos = [x for x,y in enumerate(self.get_samples_dp()) if y < argsDICT['minDP']]
            maxDP_NCpos = [x for x,y in enumerate(self.get_samples_dp()) if y > argsDICT['maxDP']]
            minGQ_NCpos = [x for x,y in enumerate(self.get_samples_gq()) if y < argsDICT['minGQ']]
            #combine and unique NotCalled positions
            NCpos = [x for x in set(minDP_NCpos+maxDP_NCpos+minGQ_NCpos)]
            for x in NCpos:
                samplesgt[x]='./.'
            for i in range(0, len(self.SAMPLES)):
                self.replace_format_on_samples(self.SAMPLES.keys()[i], 'GT', samplesgt[i])
    
    
    def get_consensus_calldata(self, argsDICT):
        """ Return consensus of existing record samples """
        consensusdict = collections.OrderedDict()
        for rf in self.FORMAT:
            consensusdict[rf] = None
        consensusdict['AD'] = self.merge_ad()
        consensusdict['DP'] = np.sum(consensusdict['AD'])
        consensusdict['GQ'] = self.merge_gq()
        #evaluate CONSENSUSPL from SAMPLESPL AND SAMPLESGTNAPOS
        consensusdict['PL'] = self.merge_pl()
        consensusdict['GT'] = self.set_consensus_gt_by_ad(consensusdict)
        return consensusdict
    
    
    def set_consensus_gt_by_ad(self, consensusdict):
        consensusad_ref = consensusdict['AD'][0]
        consensusad_alt = consensusdict['AD'][1:]
        #NOTE np.argmax return the first element for equal values which would bias GT setting always to the first ALT entry if equal depth exist
        consensusad_alt_maxpos = np.argmax(consensusad_alt)
        consensusad_alt_max = np.max(consensusad_alt)
        consensusad_ref_alt_max = [consensusad_ref, consensusad_alt_max]
        consensusad_ref_alt_maxpos = [0, consensusad_alt_maxpos+1]
        consensusgt = evaluate_ref_alt_max(consensusad_ref_alt_max, consensusad_ref_alt_maxpos)
        return consensusgt
    
    
    def merge_ad(self):
        ADOUT = []
        for i in range(0, self.get_var_number()):
            ADOUT.append(np.sum([int(y[i]) for x,y in enumerate(self.get_samples_ad()) if x not in self.get_gt_na_pos()]))
        return ADOUT
    
    
    def merge_gq(self):
        return int(np.median([[0 if z == "." else int(z) for z in y] for x,y in enumerate(self.get_samples_gq()) if x not in self.get_gt_na_pos()]))
    
    
    def merge_dp(self):
        return np.sum([y for x,y in enumerate(self.get_samples_dp()) if x not in self.get_gt_na_pos()])
    
    
    def merge_pl(self):
        PLOUT = []
        for i in range(0,nCk(self.get_var_number()+1,2)):
            PLOUT.append(int(np.median([y[i] for x,y in enumerate(self.get_samples_pl()) if x not in self.get_gt_na_pos()])))
        return PLOUT
    
    
    def get_alleles_muSNP(self):
        MAJOUT = list(np.repeat('N',len(self.SAMPLES)))
        MINOUT = list(np.repeat('N',len(self.SAMPLES)))
        IUPACOUT = list(np.repeat('N',len(self.SAMPLES)))
        for s_pos, sample_ in enumerate(self.SAMPLES.keys()):
            MAJOUT[s_pos], MINOUT[s_pos], IUPACOUT[s_pos] = _muSNP(self, sample_, 'refmajorsample')
        return [MAJOUT,MINOUT,IUPACOUT]
    
    def get_alleles_biSNP(self):
        MAJOUT = list(np.repeat('N',len(self.SAMPLES)))
        MINOUT = list(np.repeat('N',len(self.SAMPLES)))
        IUPACOUT = list(np.repeat('N',len(self.SAMPLES)))
        for s_pos, sample_ in enumerate(self.SAMPLES.keys()):
            MAJOUT[s_pos], MINOUT[s_pos], IUPACOUT[s_pos] = _biSNP(self, sample_, 'refmajorsample')
        return [MAJOUT,MINOUT,IUPACOUT]
    
    
    def get_alleles_repSNP(self):
        MAJOUT = list(np.repeat('N'*len(get_ref_alt(self)[0]),len(self.SAMPLES)))
        MINOUT = list(np.repeat('N'*len(get_ref_alt(self)[0]),len(self.SAMPLES)))
        IUPACOUT = list(np.repeat('N'*len(get_ref_alt(self)[0]),len(self.SAMPLES)))
        for s_pos, sample_ in enumerate(self.SAMPLES.keys()):
            MAJOUT[s_pos], MINOUT[s_pos], IUPACOUT[s_pos] = _repSNP(self, sample_, 'refmajorsample')
        return [MAJOUT,MINOUT,IUPACOUT]

def _biSNP(self, sample_pos, type_):
    if sample_pos not in self.SAMPLES.keys():
        return
    samplesad_ = self.SAMPLES[sample_pos]['AD']
    ref_alt_ = self.get_ref_alt()
    ref_ = ref_alt_[0]
    alt_ = ref_alt_[1:]
    ref_out = 'N'
    alt_out = 'N'
    iupac_out = 'N'
    samplesad_ = [int(x) for x in samplesad_]
    ref_alt_nonzero = [ref_alt_[x] for x,y in enumerate(samplesad_) if samplesad_[x]!=0]
    samplesad_nonzero = [samplesad_[x] for x,y in enumerate(samplesad_) if samplesad_[x]!=0]
    if len(samplesad_nonzero) == 0:
        return [ref_out, alt_out, iupac_out]
    admax = np.max([int(x) for x in samplesad_nonzero])
    admin = np.min([int(x) for x in samplesad_nonzero])
    admaxpos = [x for x,y in enumerate(samplesad_nonzero) if samplesad_nonzero[x]==admax]
    adminpos = [x for x,y in enumerate(samplesad_nonzero) if samplesad_nonzero[x]!=admax]
    if type_ == 'refaltN':
        if len(ref_alt_nonzero)==1:
            ref_out = ref_alt_nonzero[0]
            alt_out = ref_alt_nonzero[0]
            iupac_out = ref_alt_nonzero[0]
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2:
            return [ref_out, alt_out, iupac_out]
    if type_ == 'refaltsample':
        if len(ref_alt_nonzero)==1:
            ref_out = ref_alt_nonzero[0]
            alt_out = ref_alt_nonzero[0]
            iupac_out = ref_alt_nonzero[0]
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2:
            ref_out = random.sample(ref_alt_nonzero,1)[0]
            alt_out = [x for x in ref_alt_nonzero if x != ref_out][0]
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]
    if type_ == 'refmajorsample':
        if len(ref_alt_nonzero)==1:
            ref_out = ref_alt_nonzero[0]
            alt_out = ref_alt_nonzero[0]
            iupac_out = ref_alt_nonzero[0]
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2 and len(set(samplesad_nonzero)) != 1:
            ref_out = ref_alt_nonzero[admaxpos[0]]
            alt_out = ref_alt_nonzero[adminpos[0]]
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2 and len(set(samplesad_nonzero)) == 1:
            ref_out = random.sample(ref_alt_nonzero,1)[0]
            alt_out = [x for x in ref_alt_nonzero if x != ref_out][0]
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]
    if type_ == 'refminorsample':
        if len(ref_alt_nonzero)==1:
            ref_out = ref_alt_nonzero[0]
            alt_out = ref_alt_nonzero[0]
            iupac_out = ref_alt_nonzero[0]
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2 and len(set(samplesad_nonzero)) != 1:
            ref_out = ref_alt_nonzero[adminpos[0]]
            alt_out = ref_alt_nonzero[admaxpos[0]]
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2 and len(set(samplesad_nonzero)) == 1:
            ref_out = random.sample(ref_alt_nonzero,1)[0]
            alt_out = [x for x in ref_alt_nonzero if x != ref_out][0]
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]
    if type_ == 'iupac':
        if len(ref_alt_nonzero)==1:
            ref_out = ref_alt_nonzero[0]
            alt_out = ref_alt_nonzero[0]
            iupac_out = ref_alt_nonzero[0]
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2:
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]


def _muSNP(self, sample_pos, type_):
    if sample_pos not in self.SAMPLES.keys():
        return
    samplesad_ = self.SAMPLES[sample_pos]['AD']
    ref_alt_ = self.get_ref_alt()
    ref_ = ref_alt_[0]
    alt_ = ref_alt_[1:]
    ref_out = 'N'
    alt_out = 'N'
    iupac_out = 'N'
    samplesad_ = [int(x) for x in samplesad_]
    ref_alt_nonzero = [ref_alt_[x] for x,y in enumerate(samplesad_) if samplesad_[x]!=0]
    samplesad_nonzero = [samplesad_[x] for x,y in enumerate(samplesad_) if samplesad_[x]!=0]
    if len(samplesad_nonzero) == 0:
        return [ref_out, alt_out, iupac_out]
    admax = np.max([int(x) for x in samplesad_nonzero])
    admin = np.min([int(x) for x in samplesad_nonzero])
    admaxpos = [x for x,y in enumerate(samplesad_nonzero) if samplesad_nonzero[x]==admax]
    adminpos = [x for x,y in enumerate(samplesad_nonzero) if samplesad_nonzero[x]!=admax]
    if type_ == 'refaltN':
        if len(ref_alt_nonzero)==1:
            ref_out = ref_alt_nonzero[0]
            alt_out = ref_alt_nonzero[0]
            iupac_out = ref_alt_nonzero[0]
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)>1:
            return [ref_out, alt_out, iupac_out]
    if type_ == 'refaltsample':
        if len(ref_alt_nonzero)==1:
            ref_out = ref_alt_nonzero[0]
            alt_out = ref_alt_nonzero[0]
            iupac_out = ref_alt_nonzero[0]
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2:
            ref_out = random.sample(ref_alt_nonzero,1)[0]
            alt_out = [x for x in ref_alt_nonzero if x != ref_out][0]
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)>2:
            ref_out = random.sample(ref_alt_nonzero,1)[0]
            alt_out = random.sample([x for x in ref_alt_nonzero if x != ref_out],1)[0]
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]
    if type_ == 'refmajorsample':
        if len(ref_alt_nonzero)==1:
            ref_out = ref_alt_nonzero[0]
            alt_out = ref_alt_nonzero[0]
            iupac_out = ref_alt_nonzero[0]
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2 and len(set(samplesad_nonzero)) != 1:
            ref_out = ref_alt_nonzero[admaxpos[0]]
            alt_out = ref_alt_nonzero[adminpos[0]]
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2 and len(set(samplesad_nonzero)) == 1:
            ref_out = random.sample(ref_alt_nonzero,1)[0]
            alt_out = [x for x in ref_alt_nonzero if x != ref_out][0]
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)>2 and len(set(samplesad_nonzero)) != 1:
            max_nonzero = [x for y,x in enumerate(ref_alt_nonzero) if y in admaxpos]
            min_nonzero = [x for y,x in enumerate(ref_alt_nonzero) if y in adminpos]
            ref_out = random.sample(max_nonzero,1)[0]
            alt_out = random.sample(min_nonzero,1)[0]
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)>2 and len(set(samplesad_nonzero)) == 1:
            max_nonzero = [x for y,x in enumerate(ref_alt_nonzero) if y in admaxpos]
            min_nonzero = [x for y,x in enumerate(ref_alt_nonzero) if y in adminpos]
            ref_out = random.sample(max_nonzero,1)[0]
            alt_out = random.sample([x for x in max_nonzero if x != ref_out],1)[0]
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]
    if type_ == 'refminorsample':
        if len(ref_alt_nonzero)==1:
            ref_out = ref_alt_nonzero[0]
            alt_out = ref_alt_nonzero[0]
            iupac_out = ref_alt_nonzero[0]
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2 and len(set(samplesad_nonzero)) != 1:
            ref_out = ref_alt_nonzero[adminpos[0]]
            alt_out = ref_alt_nonzero[admaxpos[0]]
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2 and len(set(samplesad_nonzero)) == 1:
            ref_out = random.sample(ref_alt_nonzero,1)[0]
            alt_out = [x for x in ref_alt_nonzero if x != ref_out][0]
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)>2 and len(set(samplesad_nonzero)) != 1:
            max_nonzero = [x for y,x in enumerate(ref_alt_nonzero) if y in adminpos]
            min_nonzero = [x for y,x in enumerate(ref_alt_nonzero) if y in admaxpos]
            ref_out = random.sample(max_nonzero,1)[0]
            alt_out = random.sample(min_nonzero,1)[0]
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)>2 and len(set(samplesad_nonzero)) == 1:
            max_nonzero = [x for y,x in enumerate(ref_alt_nonzero) if y in adminpos]
            min_nonzero = [x for y,x in enumerate(ref_alt_nonzero) if y in admaxpos]
            ref_out = random.sample(min_nonzero,1)[0]
            alt_out = random.sample([x for x in min_nonzero if x != ref_out],1)[0]
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]
    if type_ == 'iupac':
        if len(ref_alt_nonzero)==1:
            ref_out = ref_alt_nonzero[0]
            alt_out = ref_alt_nonzero[0]
            iupac_out = ref_alt_nonzero[0]
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)>1:
            iupac_out = getIUPAC(''.join(ref_alt_nonzero))
            return [ref_out, alt_out, iupac_out]


def _repSNP(self, sample_pos, type_):
    if sample_pos not in self.SAMPLES.keys():
        return
    samplesad_ = self.SAMPLES[sample_pos]['AD']
    ref_alt_ = self.get_ref_alt()
    ref_ = ref_alt_[0]
    alt_ = ref_alt_[1:]
    ref_out = ''.join(list(repeat('N',len(ref_))))
    alt_out = ''.join(list(repeat('N',len(ref_))))
    iupac_out = ''.join(list(repeat('N',len(ref_))))
    samplesad_ = [int(x) for x in samplesad_]
    ref_alt_nonzero = [ref_alt_[x] for x,y in enumerate(samplesad_) if samplesad_[x]!=0]
    samplesad_nonzero = [samplesad_[x] for x,y in enumerate(samplesad_) if samplesad_[x]!=0]
    if len(samplesad_nonzero) == 0:
        return [ref_out, alt_out, iupac_out]
    admax = np.max([int(x) for x in samplesad_nonzero])
    admin = np.min([int(x) for x in samplesad_nonzero])
    admaxpos = [x for x,y in enumerate(samplesad_nonzero) if samplesad_nonzero[x]==admax]
    adminpos = [x for x,y in enumerate(samplesad_nonzero) if samplesad_nonzero[x]==admin]
    if type_ == 'refaltN':
        if len(ref_alt_nonzero)==1:
            ref_out = ref_alt_nonzero[0]
            alt_out = ref_alt_nonzero[0]
            iupac_out = ref_alt_nonzero[0]
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2:
            return [ref_out, alt_out, iupac_out]
    if type_ == 'refaltsample':
        if len(ref_alt_nonzero)==1:
            ref_out = ref_alt_nonzero[0]
            alt_out = ref_alt_nonzero[0]
            iupac_out = ref_alt_nonzero[0]
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2:
            ref_out = random.sample(ref_alt_nonzero,1)[0]
            alt_out = [x for x in ref_alt_nonzero if x != ref_out][0]
            iupac_out = ''.join([getIUPAC(ref_alt_nonzero[0][x]+ref_alt_nonzero[1][x]) for x,y in enumerate(ref_alt_nonzero[0])])
            return [ref_out, alt_out, iupac_out]
    if type_ == 'refmajorsample':
        if len(ref_alt_nonzero)==1:
            ref_out = ref_alt_nonzero[0]
            alt_out = ref_alt_nonzero[0]
            iupac_out = ref_alt_nonzero[0]
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2 and len(set(samplesad_nonzero)) != 1:
            ref_out = ref_alt_nonzero[admaxpos[0]]
            alt_out = ref_alt_nonzero[adminpos[0]]
            iupac_out = ''.join([getIUPAC(ref_alt_nonzero[0][x]+ref_alt_nonzero[1][x]) for x,y in enumerate(ref_alt_nonzero[0])])
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2 and len(set(samplesad_nonzero)) == 1:
            ref_out = random.sample(ref_alt_nonzero,1)[0]
            alt_out = [x for x in ref_alt_nonzero if x != ref_out][0]
            iupac_out = ''.join([getIUPAC(ref_alt_nonzero[0][x]+ref_alt_nonzero[1][x]) for x,y in enumerate(ref_alt_nonzero[0])])
            return [ref_out, alt_out, iupac_out]
    if type_ == 'refminorsample':
        if len(ref_alt_nonzero)==1:
            ref_out = ref_alt_nonzero[0]
            alt_out = ref_alt_nonzero[0]
            iupac_out = ref_alt_nonzero[0]
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2 and len(set(samplesad_nonzero)) != 1:
            ref_out = ref_alt_nonzero[adminpos[0]]
            alt_out = ref_alt_nonzero[admaxpos[0]]
            iupac_out = ''.join([getIUPAC(ref_alt_nonzero[0][x]+ref_alt_nonzero[1][x]) for x,y in enumerate(ref_alt_nonzero[0])])
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2 and len(set(samplesad_nonzero)) == 1:
            ref_out = random.sample(ref_alt_nonzero,1)[0]
            alt_out = [x for x in ref_alt_nonzero if x != ref_out][0]
            iupac_out = ''.join([getIUPAC(ref_alt_nonzero[0][x]+ref_alt_nonzero[1][x]) for x,y in enumerate(ref_alt_nonzero[0])])
            return [ref_out, alt_out, iupac_out]
    if type_ == 'iupac':
        if len(ref_alt_nonzero)==1:
            ref_out = ref_alt_nonzero[0]
            alt_out = ref_alt_nonzero[0]
            iupac_out = ref_alt_nonzero[0]
            return [ref_out, alt_out, iupac_out]
        if len(ref_alt_nonzero)==2:
            iupac_out = ''.join([getIUPAC(ref_alt_nonzero[0][x]+ref_alt_nonzero[1][x]) for x,y in enumerate(ref_alt_nonzero[0])])
            return [ref_out, alt_out, iupac_out]



def parse_record(self, argsDICT, polytypeDICT, filterDICT, tsvouthandle, vcfouthandle):
    """ Parse each record and writes into outhandle if variant is kept after filtering """
    REF_ALT = self.get_ref_alt()
    POLYTYPE = self.get_var_type()
    POLYNUMBER = self.get_var_number()
    POLYLEN = self.get_var_len()
    if POLYTYPE in polytypeDICT:
        polytypeDICT[POLYTYPE]+=1
    if POLYTYPE not in polytypeDICT:
        polytypeDICT[POLYTYPE]=1
    #0. REDUCE record TO SAMPLES and ADD mandatory FORMAT TO SAMPLES [GT:AD:DP]:GQ:PL
    self.keep_samples(argsDICT)
    if 'GQ' not in self.FORMAT:
        for s in argsDICT['samples_pos']:
            self.add_format_on_samples(s, 'GQ', 0)
        self.FORMAT.append('GQ')
    if 'PL' not in self.FORMAT:
        for s in argsDICT['samples_pos']:
            self.add_format_on_samples(s, 'PL', list(np.repeat(0,nCk(self.get_var_number()+1,2))))
        self.FORMAT.append('PL')
    #1. FILTER BY filterarg
    if self.is_filtered(argsDICT):
        if argsDICT['verbose']:
            print self
        filterDICT['filterarg']+=1
        return
    #2. FILTER BY minMQ
    MQ = self.get_info('MQ')
    if MQ == ".":
        MQ = 0.0
    if float(MQ) < argsDICT['minMQ']:
        if argsDICT['verbose']:
            print [self.CHROM, self.POS, self.REF, self.ALT, self.INFO, self.FORMAT, self.SAMPLES]
        filterDICT['MQ']+=1
        return
    #3. GET AD
    #fill empty AD
    self.fill_empty_ad()
    SAMPLESAD = self.get_samples_ad()
    SAMPLESADLEN = self.get_samples_ad_len()
    if not all([x==POLYNUMBER for x in SAMPLESADLEN]):
        if argsDICT['verbose']:
            print [self.CHROM, self.POS, self.REF, self.ALT, self.INFO, self.FORMAT, self.SAMPLES]
        sys.exit('\nSAMPLES AD length differs from total POLYNUMBER')
    #4. SET AD BY MISSING GT OPTION
    SAMPLESGT = self.get_samples_gt()
    SAMPLESGQ = self.get_samples_gq()
    self.set_samples_ad_by_gt(argsDICT)
    SAMPLESAD = self.get_samples_ad()
    #set dp based on altered ad
    self.set_dp_by_ad()
    #5. KEEP OR RESET GT USING THE DP OPTION AND AD INFORMATION
    SAMPLESDP = self.get_samples_dp()
    self.set_samples_gt_by_ad(argsDICT)
    #6. FILTER BY NotCalledFraction
    SAMPLESGTTYPE = self.get_gt_type()
    SAMPLESGTNAPOS = self.get_gt_na_pos()
    NOTCALLED = len(SAMPLESGTNAPOS)
    NOTCALLEDOUT = str(NOTCALLED)+'/'+str(len(self.SAMPLES))
    self.add_info(';NCFnumber='+str(NOTCALLEDOUT))
    NOTCALLEDFRACTION = float(NOTCALLED)/float(len(self.SAMPLES))
    self.add_info(';NCFfraction='+str(NOTCALLEDFRACTION))
    if NOTCALLEDFRACTION >= argsDICT['ncf']:
        if argsDICT['verbose']:
            print [self.CHROM, self.POS, self.REF, self.ALT, self.INFO, self.FORMAT, self.SAMPLES]
        filterDICT['NCF']+=1
        return
    #7. GET CONSENSUS AND EVALUATE CDP OPTION
    #fill empty PL
    self.fill_empty_pl()    
    CONSENSUS_CALLDATA = self.get_consensus_calldata(argsDICT)
    if CONSENSUS_CALLDATA['DP'] < argsDICT['cdp']:
        if argsDICT['verbose']:
            print [self.CHROM, self.POS, self.REF, self.ALT, self.INFO, self.FORMAT, self.SAMPLES]
        filterDICT['CDP']+=1
        return
    if argsDICT['add']:
        self.SAMPLES[argsDICT['id']] = CONSENSUS_CALLDATA
    if not argsDICT['add']:
        self.SAMPLES = collections.OrderedDict()
        self.SAMPLES[argsDICT['id']] = CONSENSUS_CALLDATA
    VCFLINE = '\t'.join([self.CHROM, self.POS, self.ID, self.REF[0], ','.join(self.ALT), self.QUAL, self.FILTER, self.INFO, ':'.join(self.FORMAT)])
    SAMPLES_LINE = []
    for s in self.SAMPLES.keys():
        SAMPLE_LINE = []
        for k in self.SAMPLES[s].keys():
            if isinstance(self.SAMPLES[s][k], list):
                SAMPLE_LINE.append(','.join([str(x) for x in self.SAMPLES[s][k]]))
            elif self.SAMPLES[s][k] is None:
                SAMPLE_LINE.append(".")
            else:
                SAMPLE_LINE.append(str(self.SAMPLES[s][k]))
        SAMPLES_LINE.append(':'.join(SAMPLE_LINE))
    vcfouthandle.write(VCFLINE+'\t'+'\t'.join(SAMPLES_LINE)+'\n')
    CONSENSUSOUT = SAMPLES_LINE[-1]
    if POLYTYPE == 'repSNP':
        if argsDICT['add']:
            CONSENSUSMAJ, CONSENSUSMIN, CONSENSUSIUPAC = self.get_alleles_repSNP()
            CONSENSUSMAJ = CONSENSUSMAJ[[x for x,y in enumerate(self.SAMPLES.keys()) if y == argsDICT['id']][0]]
            CONSENSUSMIN = CONSENSUSMIN[[x for x,y in enumerate(self.SAMPLES.keys()) if y == argsDICT['id']][0]]
            CONSENSUSIUPAC = CONSENSUSIUPAC[[x for x,y in enumerate(self.SAMPLES.keys()) if y == argsDICT['id']][0]]
        if not argsDICT['add']:
            CONSENSUSMAJ, CONSENSUSMIN, CONSENSUSIUPAC = self.get_alleles_repSNP()
            CONSENSUSMAJ = CONSENSUSMAJ[0]
            CONSENSUSMIN = CONSENSUSMIN[0]
            CONSENSUSIUPAC = CONSENSUSIUPAC[0]
    if POLYTYPE == 'biSNP':
        if argsDICT['add']:
            CONSENSUSMAJ, CONSENSUSMIN, CONSENSUSIUPAC = self.get_alleles_biSNP()
            CONSENSUSMAJ = CONSENSUSMAJ[[x for x,y in enumerate(self.SAMPLES.keys()) if y == argsDICT['id']][0]]
            CONSENSUSMIN = CONSENSUSMIN[[x for x,y in enumerate(self.SAMPLES.keys()) if y == argsDICT['id']][0]]
            CONSENSUSIUPAC = CONSENSUSIUPAC[[x for x,y in enumerate(self.SAMPLES.keys()) if y == argsDICT['id']][0]]
        if not argsDICT['add']:
            CONSENSUSMAJ, CONSENSUSMIN, CONSENSUSIUPAC = self.get_alleles_biSNP()
            CONSENSUSMAJ = CONSENSUSMAJ[0]
            CONSENSUSMIN = CONSENSUSMIN[0]
            CONSENSUSIUPAC = CONSENSUSIUPAC[0]
    if POLYTYPE == 'muSNP':
        if argsDICT['add']:
            CONSENSUSMAJ, CONSENSUSMIN, CONSENSUSIUPAC = self.get_alleles_muSNP()
            CONSENSUSMAJ = CONSENSUSMAJ[[x for x,y in enumerate(self.SAMPLES.keys()) if y == argsDICT['id']][0]]
            CONSENSUSMIN = CONSENSUSMIN[[x for x,y in enumerate(self.SAMPLES.keys()) if y == argsDICT['id']][0]]
            CONSENSUSIUPAC = CONSENSUSIUPAC[[x for x,y in enumerate(self.SAMPLES.keys()) if y == argsDICT['id']][0]]
        if not argsDICT['add']:
            CONSENSUSMAJ, CONSENSUSMIN, CONSENSUSIUPAC = self.get_alleles_muSNP()
            CONSENSUSMAJ = CONSENSUSMAJ[0]
            CONSENSUSMIN = CONSENSUSMIN[0]
            CONSENSUSIUPAC = CONSENSUSIUPAC[0]
    if POLYTYPE != 'repSNP' and POLYTYPE != 'biSNP' and POLYTYPE != 'muSNP':
        print POLYTYPE
        sys.exit('\nPOLYTYPE not considered for CONSENSUS\nPlease reduce to SNP and MNP')
    #write to tsv
    TSVLINE = '\t'.join([self.CHROM,str(self.POS),NOTCALLEDOUT,CONSENSUSOUT,CONSENSUSMAJ,CONSENSUSMIN,CONSENSUSIUPAC])+'\n'
    tsvouthandle.write(TSVLINE)
    filterDICT['passed']+=1


def main():
    """ main """
    parser = argparse.ArgumentParser(usage='%(prog)s [options]',
                                     description='Given a multiple sample vcf file reduced to SNP and MNP and the reference fasta file, ' \
                                                 'this script merges the samples into a consensus. ' \
                                                 'Major, minor and IUPAC code of the consensus will be written to a tab separated file.'
                                                 '\nTo extract a subpopulations from e.g. GATK multiple vcf file use the following command from GATK: ' \
                                                 'java -jar GenomeAnalysisTK.jar -T SelectVariants -sn $SAMPLE1 -sn $SAMPLE2 -V $MULTISAMPLEVCF -R $REFERENCE ' \
                                                 '[-trimAlternates : removeUnusedAlternates] ' \
                                                 '[-env : excludeNonVariants] ' \
                                                 '[-ef : excludeFiltered] ' \
                                                 '-selectType SNP -selectType MNP -o $OUTPUT ' \
                                                 '\nor vcftools: ' \
                                                 'vcftools --vcf $MULTISAMPLEVCF [--gvcf] --out $OUTPUT --recode --indv $SAMPLE1 --indc $SAMPLE2 ' \
                                                 '[--minGQ : filter for genotype quality] ' \
                                                 '[--minDP : filter for depth] ' \
                                                 '--remove-indels ')
    parser.add_argument('-v', '--verbose', help='increase output verbosity', action='store_true')
    #parser.add_argument('-tsv', help='specify if tsv should be written [default: False]', action='store_true')
    parser.add_argument('-tabix', help='specify if input is tabix file. ' \
                                       'NOTE: if set to true reference file needs to be specified to get length of each entry [default: False]'
                                       , action='store_true')
    parser.add_argument('-R', help='specify reference fasta file. Only needed if [tabix] option')
    parser.add_argument('-ivcf', help='specify input vcf file')
    parser.add_argument('-o', help='specify output prefix')
    parser.add_argument('-H', help='specify if HEADER of vcf file should be ommitted [default: False]', action='store_true')
    parser.add_argument('-chr', help='secify chromosomes to be parsed as a comma-seperated list')
    parser.add_argument('-filterarg', help='specify FILTER arguments to keep variants, all variants without these will be discarded. ' \
                                           'use a comma-seperated list e.g. PASS,VQSRTTranchesSNP90. ' \
                                           'By default all variants are retained [default: ""]', default='')
    parser.add_argument('-samples', help='specify sample names to be merged as a comma-seperated list, they need to exactly match sample ids from vcf input')
    parser.add_argument('-add', help='specify if CONSENSUS should be added as a new column to existing vcf file [default: False]', action='store_true')
    parser.add_argument('-id', help='specify CONSENSUS name to be added to HEADER')
    parser.add_argument('-minMQ', help='specify mapping quality threshold necessary to exceeded to retain variant [default: "0.0"]', type=float, default=0.0)
    parser.add_argument('-missingGT', help='specify how missing GT data (based on e.g. calling SNPs with GATK on multiple samples simultaneously; GT: ./.) ' \
                                         'should be handled, either keep AD coverage [default: "keep"], set AD to 0 for all samples which have GT: ./. for all alleles ["zero"], ' \