This function will parse the vcf to detect sites having multiple alt alleles and pick out on the
    most likely ones.
    INFILE : open file handle for the input vcf
    OUTFILE : open file handle for the output vcf

    The columns in INFILE are
     [0] CHROM
     [1] POS
     [2] ID
     [3] REF
     [4] ALT
     [5] QUAL
     [6] FILTER
     [7] INFO
     [8] FORMAT
     [9] DNA_NORMAL
    [10] DNA_TUMOR
    [11] RNA_TUMOR  -  Not always present

def parse_radia_multi_alt(infile, outfile):
    """
    This function will parse the vcf to detect sites having multiple alt alleles and pick out on the
    most likely ones.
    INFILE : open file handle for the input vcf
    OUTFILE : open file handle for the output vcf

    The columns in INFILE are
     [0] CHROM
     [1] POS
     [2] ID
     [3] REF
     [4] ALT
     [5] QUAL
     [6] FILTER
     [7] INFO
     [8] FORMAT
     [9] DNA_NORMAL
    [10] DNA_TUMOR
    [11] RNA_TUMOR  -  Not always present
    """
    for line in infile:
        # Print header to putfile
        if line.startswith('#'):
            print(line.strip(), file=outfile)
            continue
        line = line.strip().split('\t')
        # If there is just 1 ALT allele, print and continue
        if len(line[4]) == 1:
            print('\t'.join(line), file=outfile)
        # If not, process
        else:
            seq_field_indeces = [9, 10]
            alleles = [line[3]] + line[4].split(',')  # all alleles, incl. REF
            # collect tumor, normal and (if present) rna AD and AFs
            # AD = Depth of reads supporting each allele
            # AF = Fraction of reads supporting each allele
            normal_AD = line[9].split(':')[5].split(',')
            normal_AF = line[9].split(':')[6].split(',')
            tumor_AD = line[10].split(':')[5].split(',')
            tumor_AF = line[10].split(':')[6].split(',')
            if len(line[11]) > 1:
                rna_AD = line[11].split(':')[5].split(',')
                rna_AF = line[11].split(':')[6].split(',')
                seq_field_indeces += [11]  # append rna since it is present
            else:
                # If rna is missing, set RNA_AD and RNA_AF to null sets for easily integrating into
                # the logic in the following code
                rna_AD = rna_AF = [0, 0, 0, 0]
            # Initialise variables to store the probable ALT alleles and the index values of the
            # same wrt AD and AF
            out_alleles = set([])
            out_AF_AD_index = {0}
            # parse AD and AF to get most probable ALT alleles
            for i in range(1, len(normal_AF)):
                # Criteria for selection = AD > 4 and AF >0.1 in either tumor or RNA, given normal
                # AF < 0.1
                if ((float(tumor_AF[i]) >= 0.1 and int(tumor_AD[i]) >= 4) \
                        or (float(rna_AF[i]) >= 0.1 and int(rna_AD[i]) >= 4)) \
                        and (float(normal_AF[i]) < 0.1):
                    out_alleles.add(alleles[i])
                    out_AF_AD_index.add(i)
            # If the number of probable alleles is greater than 0 the print to outfile with the
            # modified allele fraction representing reads corrresponding to all alleles
            if len(out_alleles) > 0:
                line[4] = ','.join(out_alleles)  # set alt alleles
                # Modify the AD and AF values in the TUMOR/NORMAL/RNA fields
                # one at a time.  Seq fields contain
                # [0] GT* - Genotype
                # [1] DP - Read depth at this position in the sample
                # [2] INDEL - Number of indels
                # [3] START - Number of reads starting at this position
                # [4] STOP - Number of reads stopping at this position
                # [5] AD* - Depth of reads supporting alleles
                # [6] AF* - Fraction of reads supporting alleles
                # [7] BQ* - Avg base quality for reads supporting alleles
                # [8] SB* - Strand Bias for reads supporting alleles
                # Fields marked with *s are teh ones that contain info for each seq field and need
                # to be modified
                for seq_field_index in seq_field_indeces:
                    # Get the details for seq_field
                    deets = line[seq_field_index].split(':')
                    # modify fields 5 thu 8 to hold only info for the probable
                    # alleles
                    for field_index in range(5, 9):
                        field = deets[field_index].split(",")
                        deets[field_index] = ",".join([x for i, x in enumerate(field)
                                                       if i in out_AF_AD_index])
                    # Modify DP to hold the new total of reads
                    deets[1] = str(sum([int(x) for x in deets[5].split(",")]))
                    # get the most likely genotypes based on AD and AF
                    GT_by_AD = set([i for i, x in enumerate(deets[5].split(",")) if int(x) >= 4])
                    GT_by_AF = set([i for i, x in enumerate(deets[6].split(",")) \
                                    if float(x) >= 0.1])
                    # Get the consensus genotype
                    GT = GT_by_AD.intersection(GT_by_AF)
                    if len(GT) == 0:
                        deets[0] = "0/0"
                    elif len(GT) == 1:
                        deets[0] = "/".join([str(x) for x in GT] + [str(x) for x in GT])
                    elif len(GT) == 2:
                        deets[0] = "/".join([str(x) for x in GT])
                    else:
                        print("ERROR : triple genotype detected", file=sys.stderr)
                        print(line, file=sys.stdout)
                    # Rejoin the details line
                    line[seq_field_index] = ":".join(deets)
                # Print the modified line to output
                print("\t".join(line), file=outfile)
            # Else do nothing
            else:
                pass