treetime/vcf_utils.py
def write_vcf(tree_dict, file_name):#, compress=False):
    """
    Writes out a VCF-style file (which seems to be minimally handleable
    by vcftools and pyvcf) of the alignment. This is created from a dict
    in a similar format to what's created by :py:meth:`treetime.vcf_utils.read_vcf`

    Positions of variable sites are transformed to start at 1 to match
    VCF convention.

    Parameters
    ----------
     tree_dict: nested dict
        A nested dict with keys 'sequence' 'reference' and 'positions',
        as is created by :py:meth:`treetime.TreeAnc.get_tree_dict`

     file_name: str
        File to which the new VCF should be written out. File names ending with
        '.gz' will result in the VCF automatically being gzipped.

    """

#   Programming Logic Note:
#
#    For a sequence like:
#    Pos     1 2 3 4 5 6
#    Ref     A C T T A C
#    Seq1    A C - - - G
#
#    In a dict it is stored:
#    Seq1:{3:'-', 4:'-', 5:'-', 6:'G'}  (Numbering from 1 for simplicity)
#
#    In a VCF it needs to be:
#    POS REF     ALT     Seq1
#    2   CTTA    C       1/1
#    6   C       G       1/1
#
#    If a position is deleted (pos 3), need to get invariable position preceeding it
#
#    However, in alternative case, the base before a deletion is mutant, so need to check
#        that next position isn't a deletion (as otherwise won't be found until after the
#        current single bp mutation is written out)
#
#    When deleted position found, need to gather up all adjacent mutant positions with deletions,
#        but not include adjacent mutant positions that aren't deletions (pos 6)
#
#    Don't run off the 'end' of the position list if deletion is the last thing to be included
#        in the VCF file

    sequences = tree_dict['sequences']
    ref = tree_dict['reference']
    positions = tree_dict['positions']

    def handleDeletions(i, pi, pos, ref, delete, pattern):
        refb = ref[pi]
        if delete: #Need to get the position before
            i-=1    #As we'll next go to this position again
            pi-=1
            pos = pi+1
            refb = ref[pi]
            #re-get pattern
            pattern = []
            for k,v in sequences.items():
                try:
                    pattern.append(sequences[k][pi])
                except KeyError:
                    pattern.append(ref[pi])
            pattern = np.array(pattern)

        sites = []
        sites.append(pattern)

        #Gather all positions affected by deletion - but don't run off end of position list
        while (i+1) < len(positions) and positions[i+1] == pi+1:
            i+=1
            pi = positions[i]
            pattern = []
            for k,v in sequences.items():
                try:
                    pattern.append(sequences[k][pi])
                except KeyError:
                    pattern.append(ref[pi])
            pattern = np.array(pattern)

            #Stops 'greedy' behaviour from adding mutations adjacent to deletions
            if any(pattern == '-'): #if part of deletion, append
                sites.append(pattern)
                refb = refb+ref[pi]
            else: #this is another mutation next to the deletion!
                i-=1    #don't append, break this loop

        #Rotate them into 'calls'
        sites = np.asarray(sites)
        align = np.rot90(sites)
        align = np.flipud(align)

        #Get rid of '-', and put '.' for calls that match ref
        #Only removes trailing '-'. This breaks VCF convension, but the standard
        #VCF way of handling this* is really complicated, and the situation is rare.
        #(*deletions and mutations at the same locations)
        fullpat = []
        for pt in align:
            gp = len(pt)-1
            while pt[gp] == '-':
                pt[gp] = ''
                gp-=1
            pat = "".join(pt)
            if pat == refb:
                fullpat.append('.')
            else:
                fullpat.append(pat)

        pattern = np.array(fullpat)

        return i, pi, pos, refb, pattern


    #prepare the header of the VCF & write out
    header=["#CHROM","POS","ID","REF","ALT","QUAL","FILTER","INFO","FORMAT"]+list(sequences.keys())
    with open(file_name, 'w') as the_file:
        the_file.write( "##fileformat=VCFv4.2\n"+
                        "##source=NextStrain\n"+
                        "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n")
        the_file.write("\t".join(header)+"\n")

    vcfWrite = []
    errorPositions = []
    explainedErrors = 0

    #Why so basic? Because we sometimes have to back up a position!
    i=0
    while i < len(positions):
        #Get the 'pattern' of all calls at this position.
        #Look out specifically for current (this pos) or upcoming (next pos) deletions
        #But also distinguish these two, as handled differently.

        pi = positions[i]
        pos = pi+1 #change numbering to match VCF, not python, for output
        refb = ref[pi] #reference base at this position
        delete = False #deletion at this position - need to grab previous base (invariable)
        deleteGroup = False #deletion at next position (mutation at this pos) - do not need to get prev base

        #try/except is much more efficient than 'if' statements for constructing patterns,
        #as on average a 'variable' location will not be variable for any given sequence
        pattern = []
        #pattern2 gets the pattern at next position to check for upcoming deletions
        #it's more efficient to get both here rather than loop through sequences twice!
        pattern2 = []
        for k,v in sequences.items():
            try:
                pattern.append(sequences[k][pi])
            except KeyError:
                pattern.append(ref[pi])

            try:
                pattern2.append(sequences[k][pi+1])
            except KeyError:
                pattern2.append(ref[pi+1])

        pattern = np.array(pattern)
        pattern2 = np.array(pattern2)

        #If a deletion here, need to gather up all bases, and position before
        if any(pattern == '-'):
            if pos != 1:
                deleteGroup = True
                delete = True
            else:
                #If theres a deletion in 1st pos, VCF files do not handle this well.
                #Proceed keeping it as '-' for alt (violates VCF), but warn user to check output.
                #(This is rare)
                print ("WARNING: You have a deletion in the first position of your alignment. VCF format does not handle this well. Please check the output to ensure it is correct.")
        else:
            #If a deletion in next pos, need to gather up all bases
            if any(pattern2 == '-'):
                deleteGroup = True

        #If deletion, treat affected bases as 1 'call':
        if delete or deleteGroup:
            i, pi, pos, refb, pattern = handleDeletions(i, pi, pos, ref, delete, pattern)
        #If no deletion, replace ref with '.', as in VCF format
        else:
            pattern[pattern==refb] = '.'

        #Get the list of ALTs - minus any '.'!
        uniques = np.unique(pattern)
        uniques = uniques[np.where(uniques!='.')]

        #Convert bases to the number that matches the ALT
        j=1
        for u in uniques:
            pattern[np.where(pattern==u)[0]] = str(j)
            j+=1
        #Now convert these calls to #/# (VCF format)
        calls = [ j+"/"+j if j!='.' else '.' for j in pattern ]

        #What if there's no variation at a variable site??
        #This can happen when sites are modified by TreeTime - see below.
        printPos = True
        if len(uniques)==0:
            #If we expect it (it was made constant by TreeTime), it's fine.
            if 'inferred_const_sites' in tree_dict and pi in tree_dict['inferred_const_sites']:
                explainedErrors += 1
                printPos = False #and don't output position to the VCF
            else:
                #If we don't expect, raise an error
                errorPositions.append(str(pi))

        #Write it out - Increment positions by 1 so it's in VCF numbering
        #If no longer variable, and explained, don't write it out
        if printPos:
            output = ["MTB_anc", str(pos), ".", refb, ",".join(uniques), ".", "PASS", ".", "GT"] + calls
            vcfWrite.append("\t".join(output))

        i+=1

    #Note: The number of 'inferred_const_sites' passed back by TreeTime will often be longer
    #than the number of 'site that were made constant' that prints below. This is because given the site:
    # Ref   Alt     Seq
    # G     A       AANAA
    #This will be converted to 'AAAAA' and listed as an 'inferred_const_sites'. However, for VCF
    #purposes, because the site is 'variant' against the ref, it is variant, as expected, and so
    #won't be counted in the below list, which is only sites removed from the VCF.

    if 'inferred_const_sites' in tree_dict and explainedErrors != 0:
        print ( "Sites that were constant except for ambiguous bases were made constant by TreeTime. This happened {} times. These sites are now excluded from the VCF.".format(explainedErrors))

    if len(errorPositions) != 0:
        print ("\n***WARNING: vcf_utils.py"
            "\n{} sites were found that had no alternative bases. If this data has been "
            "run through TreeTime and contains ambiguous bases, try calling get_tree_dict with "
            "var_ambigs=True to see if this clears the error."
            "\n\nAlternative causes:"
            "\n- Not all sequences in your alignment are in the tree (if you are running TreeTime via commandline "
            "this is most likely)"
            "\n- In TreeTime, can be caused by overwriting variants in tips with small branch lengths (debug)"
            "\n\nThese are the positions affected (numbering starts at 0):".format(str(len(errorPositions))))
        print (",".join(errorPositions))

    with open(file_name, 'a') as the_file:
        the_file.write("\n".join(vcfWrite))

    if file_name.endswith(('.gz', '.GZ')):
        import os
        #must temporarily remove .gz ending, or gzip won't zip it!
        os.rename(file_name, file_name[:-3])
        call = ["gzip", file_name[:-3]]
        os.system(" ".join(call))
def write_vcf(tree_dict, file_name):#, compress=False):
    """
    Writes out a VCF-style file (which seems to be minimally handleable
    by vcftools and pyvcf) of the alignment. This is created from a dict
    in a similar format to what's created by :py:meth:`treetime.vcf_utils.read_vcf`

    Positions of variable sites are transformed to start at 1 to match
    VCF convention.

    Parameters
    ----------
     tree_dict: nested dict
        A nested dict with keys 'sequence' 'reference' and 'positions',
        as is created by :py:meth:`treetime.TreeAnc.get_tree_dict`

     file_name: str
        File to which the new VCF should be written out. File names ending with
        '.gz' will result in the VCF automatically being gzipped.

    """

#   Programming Logic Note:
#
#    For a sequence like:
#    Pos     1 2 3 4 5 6
#    Ref     A C T T A C
#    Seq1    A C - - - G
#
#    In a dict it is stored:
#    Seq1:{3:'-', 4:'-', 5:'-', 6:'G'}  (Numbering from 1 for simplicity)
#
#    In a VCF it needs to be:
#    POS REF     ALT     Seq1
#    2   CTTA    C       1/1
#    6   C       G       1/1
#
#    If a position is deleted (pos 3), need to get invariable position preceeding it
#
#    However, in alternative case, the base before a deletion is mutant, so need to check
#        that next position isn't a deletion (as otherwise won't be found until after the
#        current single bp mutation is written out)
#
#    When deleted position found, need to gather up all adjacent mutant positions with deletions,
#        but not include adjacent mutant positions that aren't deletions (pos 6)
#
#    Don't run off the 'end' of the position list if deletion is the last thing to be included
#        in the VCF file

    sequences = tree_dict['sequences']
    ref = tree_dict['reference']
    positions = tree_dict['positions']

    def handleDeletions(i, pi, pos, ref, delete, pattern):
        refb = ref[pi]
        if delete: #Need to get the position before
            i-=1    #As we'll next go to this position again
            pi-=1
            pos = pi+1
            refb = ref[pi]
            #re-get pattern
            pattern = []
            for k,v in sequences.items():
                try:
                    pattern.append(sequences[k][pi])
                except KeyError:
                    pattern.append(ref[pi])
            pattern = np.array(pattern)

        sites = []
        sites.append(pattern)

        #Gather all positions affected by deletion - but don't run off end of position list
        while (i+1) < len(positions) and positions[i+1] == pi+1:
            i+=1
            pi = positions[i]
            pattern = []
            for k,v in sequences.items():
                try:
                    pattern.append(sequences[k][pi])
                except KeyError:
                    pattern.append(ref[pi])
            pattern = np.array(pattern)

            #Stops 'greedy' behaviour from adding mutations adjacent to deletions
            if any(pattern == '-'): #if part of deletion, append
                sites.append(pattern)
                refb = refb+ref[pi]
            else: #this is another mutation next to the deletion!
                i-=1    #don't append, break this loop

        #Rotate them into 'calls'
        sites = np.asarray(sites)
        align = np.rot90(sites)
        align = np.flipud(align)

        #Get rid of '-', and put '.' for calls that match ref
        #Only removes trailing '-'. This breaks VCF convension, but the standard
        #VCF way of handling this* is really complicated, and the situation is rare.
        #(*deletions and mutations at the same locations)
        fullpat = []
        for pt in align:
            gp = len(pt)-1
            while pt[gp] == '-':
                pt[gp] = ''
                gp-=1
            pat = "".join(pt)
            if pat == refb:
                fullpat.append('.')
            else:
                fullpat.append(pat)

        pattern = np.array(fullpat)

        return i, pi, pos, refb, pattern


    #prepare the header of the VCF & write out
    header=["#CHROM","POS","ID","REF","ALT","QUAL","FILTER","INFO","FORMAT"]+list(sequences.keys())
    with open(file_name, 'w') as the_file:
        the_file.write( "##fileformat=VCFv4.2\n"+
                        "##source=NextStrain\n"+
                        "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n")
        the_file.write("\t".join(header)+"\n")

    vcfWrite = []
    errorPositions = []
    explainedErrors = 0

    #Why so basic? Because we sometimes have to back up a position!
    i=0
    while i < len(positions):
        #Get the 'pattern' of all calls at this position.
        #Look out specifically for current (this pos) or upcoming (next pos) deletions
        #But also distinguish these two, as handled differently.

        pi = positions[i]
        pos = pi+1 #change numbering to match VCF, not python, for output
        refb = ref[pi] #reference base at this position
        delete = False #deletion at this position - need to grab previous base (invariable)
        deleteGroup = False #deletion at next position (mutation at this pos) - do not need to get prev base

        #try/except is much more efficient than 'if' statements for constructing patterns,
        #as on average a 'variable' location will not be variable for any given sequence
        pattern = []
        #pattern2 gets the pattern at next position to check for upcoming deletions
        #it's more efficient to get both here rather than loop through sequences twice!
        pattern2 = []
        for k,v in sequences.items():
            try:
                pattern.append(sequences[k][pi])
            except KeyError:
                pattern.append(ref[pi])

            try:
                pattern2.append(sequences[k][pi+1])
            except KeyError:
                pattern2.append(ref[pi+1])

        pattern = np.array(pattern)
        pattern2 = np.array(pattern2)

        #If a deletion here, need to gather up all bases, and position before
        if any(pattern == '-'):
            if pos != 1:
                deleteGroup = True
                delete = True
            else:
                #If theres a deletion in 1st pos, VCF files do not handle this well.
                #Proceed keeping it as '-' for alt (violates VCF), but warn user to check output.
                #(This is rare)
                print ("WARNING: You have a deletion in the first position of your alignment. VCF format does not handle this well. Please check the output to ensure it is correct.")
        else:
            #If a deletion in next pos, need to gather up all bases
            if any(pattern2 == '-'):
                deleteGroup = True

        #If deletion, treat affected bases as 1 'call':
        if delete or deleteGroup:
            i, pi, pos, refb, pattern = handleDeletions(i, pi, pos, ref, delete, pattern)
        #If no deletion, replace ref with '.', as in VCF format
        else:
            pattern[pattern==refb] = '.'

        #Get the list of ALTs - minus any '.'!
        uniques = np.unique(pattern)
        uniques = uniques[np.where(uniques!='.')]

        #Convert bases to the number that matches the ALT
        j=1
        for u in uniques:
            pattern[np.where(pattern==u)[0]] = str(j)
            j+=1
        #Now convert these calls to #/# (VCF format)
        calls = [ j+"/"+j if j!='.' else '.' for j in pattern ]

        #What if there's no variation at a variable site??
        #This can happen when sites are modified by TreeTime - see below.
        printPos = True
        if len(uniques)==0:
            #If we expect it (it was made constant by TreeTime), it's fine.
            if 'inferred_const_sites' in tree_dict and pi in tree_dict['inferred_const_sites']:
                explainedErrors += 1
                printPos = False #and don't output position to the VCF
            else:
                #If we don't expect, raise an error
                errorPositions.append(str(pi))

        #Write it out - Increment positions by 1 so it's in VCF numbering
        #If no longer variable, and explained, don't write it out
        if printPos:
            output = ["MTB_anc", str(pos), ".", refb, ",".join(uniques), ".", "PASS", ".", "GT"] + calls
            vcfWrite.append("\t".join(output))

        i+=1

    #Note: The number of 'inferred_const_sites' passed back by TreeTime will often be longer
    #than the number of 'site that were made constant' that prints below. This is because given the site:
    # Ref   Alt     Seq
    # G     A       AANAA
    #This will be converted to 'AAAAA' and listed as an 'inferred_const_sites'. However, for VCF
    #purposes, because the site is 'variant' against the ref, it is variant, as expected, and so
    #won't be counted in the below list, which is only sites removed from the VCF.

    if 'inferred_const_sites' in tree_dict and explainedErrors != 0:
        print ( "Sites that were constant except for ambiguous bases were made constant by TreeTime. This happened {} times. These sites are now excluded from the VCF.".format(explainedErrors))

    if len(errorPositions) != 0:
        print ("\n***WARNING: vcf_utils.py"
            "\n{} sites were found that had no alternative bases. If this data has been "
            "run through TreeTime and contains ambiguous bases, try calling get_tree_dict with "
            "var_ambigs=True to see if this clears the error."
            "\n\nAlternative causes:"
            "\n- Not all sequences in your alignment are in the tree (if you are running TreeTime via commandline "
            "this is most likely)"
            "\n- In TreeTime, can be caused by overwriting variants in tips with small branch lengths (debug)"
            "\n\nThese are the positions affected (numbering starts at 0):".format(str(len(errorPositions))))
        print (",".join(errorPositions))

    with open(file_name, 'a') as the_file:
        the_file.write("\n".join(vcfWrite))

    if file_name.endswith(('.gz', '.GZ')):
        import os
        #must temporarily remove .gz ending, or gzip won't zip it!
        os.rename(file_name, file_name[:-3])
        call = ["gzip", file_name[:-3]]
        os.system(" ".join(call))
def write_vcf ( tree_dict , file_name ) : #, compress=False):
 #   Programming Logic Note:
 #
 #    For a sequence like:
 #    Pos     1 2 3 4 5 6
 #    Ref     A C T T A C
 #    Seq1    A C - - - G
 #
 #    In a dict it is stored:
 #    Seq1:{3:'-', 4:'-', 5:'-', 6:'G'}  (Numbering from 1 for simplicity)
 #
 #    In a VCF it needs to be:
 #    POS REF     ALT     Seq1
 #    2   CTTA    C       1/1
 #    6   C       G       1/1
 #
 #    If a position is deleted (pos 3), need to get invariable position preceeding it
 #
 #    However, in alternative case, the base before a deletion is mutant, so need to check
 #        that next position isn't a deletion (as otherwise won't be found until after the
 #        current single bp mutation is written out)
 #
 #    When deleted position found, need to gather up all adjacent mutant positions with deletions,
 #        but not include adjacent mutant positions that aren't deletions (pos 6)
 #
 #    Don't run off the 'end' of the position list if deletion is the last thing to be included
 #        in the VCF file
 sequences = tree_dict [ 'sequences' ] ref = tree_dict [ 'reference' ] positions = tree_dict [ 'positions' ] def handleDeletions ( i , pi , pos , ref , delete , pattern ) : refb = ref [ pi ] if delete : #Need to get the position before
 i -= 1 #As we'll next go to this position again
 pi -= 1 pos = pi + 1 refb = ref [ pi ] #re-get pattern
 pattern = [ ] for k , v in sequences . items ( ) : try : pattern . append ( sequences [ k ] [ pi ] ) except KeyError : pattern . append ( ref [ pi ] ) pattern = np . array ( pattern ) sites = [ ] sites . append ( pattern ) #Gather all positions affected by deletion - but don't run off end of position list
 while ( i + 1 ) < len ( positions ) and positions [ i + 1 ] == pi + 1 : i += 1 pi = positions [ i ] pattern = [ ] for k , v in sequences . items ( ) : try : pattern . append ( sequences [ k ] [ pi ] ) except KeyError : pattern . append ( ref [ pi ] ) pattern = np . array ( pattern ) #Stops 'greedy' behaviour from adding mutations adjacent to deletions
 if any ( pattern == '-' ) : #if part of deletion, append
 sites . append ( pattern ) refb = refb + ref [ pi ] else : #this is another mutation next to the deletion!
 i -= 1 #don't append, break this loop
 #Rotate them into 'calls'
 sites = np . asarray ( sites ) align = np . rot90 ( sites ) align = np . flipud ( align ) #Get rid of '-', and put '.' for calls that match ref
 #Only removes trailing '-'. This breaks VCF convension, but the standard
 #VCF way of handling this* is really complicated, and the situation is rare.
 #(*deletions and mutations at the same locations)
 fullpat = [ ] for pt in align : gp = len ( pt ) - 1 while pt [ gp ] == '-' : pt [ gp ] = '' gp -= 1 pat = "" . join ( pt ) if pat == refb : fullpat . append ( '.' ) else : fullpat . append ( pat ) pattern = np . array ( fullpat ) return i , pi , pos , refb , pattern #prepare the header of the VCF & write out
 header = [ "#CHROM" , "POS" , "ID" , "REF" , "ALT" , "QUAL" , "FILTER" , "INFO" , "FORMAT" ] + list ( sequences . keys ( ) ) with open ( file_name , 'w' ) as the_file : the_file . write ( "##fileformat=VCFv4.2\n" + "##source=NextStrain\n" + "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n" ) the_file . write ( "\t" . join ( header ) + "\n" ) vcfWrite = [ ] errorPositions = [ ] explainedErrors = 0 #Why so basic? Because we sometimes have to back up a position!
 i = 0 while i < len ( positions ) : #Get the 'pattern' of all calls at this position.
 #Look out specifically for current (this pos) or upcoming (next pos) deletions
 #But also distinguish these two, as handled differently.
 pi = positions [ i ] pos = pi + 1 #change numbering to match VCF, not python, for output
 refb = ref [ pi ] #reference base at this position
 delete = False #deletion at this position - need to grab previous base (invariable)
 deleteGroup = False #deletion at next position (mutation at this pos) - do not need to get prev base
 #try/except is much more efficient than 'if' statements for constructing patterns,
 #as on average a 'variable' location will not be variable for any given sequence
 pattern = [ ] #pattern2 gets the pattern at next position to check for upcoming deletions
 #it's more efficient to get both here rather than loop through sequences twice!
 pattern2 = [ ] for k , v in sequences . items ( ) : try : pattern . append ( sequences [ k ] [ pi ] ) except KeyError : pattern . append ( ref [ pi ] ) try : pattern2 . append ( sequences [ k ] [ pi + 1 ] ) except KeyError : pattern2 . append ( ref [ pi + 1 ] ) pattern = np . array ( pattern ) pattern2 = np . array ( pattern2 ) #If a deletion here, need to gather up all bases, and position before
 if any ( pattern == '-' ) : if pos != 1 : deleteGroup = True delete = True else : #If theres a deletion in 1st pos, VCF files do not handle this well.
 #Proceed keeping it as '-' for alt (violates VCF), but warn user to check output.
 #(This is rare)
 print ( "WARNING: You have a deletion in the first position of your alignment. VCF format does not handle this well. Please check the output to ensure it is correct." ) else : #If a deletion in next pos, need to gather up all bases
 if any ( pattern2 == '-' ) : deleteGroup = True #If deletion, treat affected bases as 1 'call':
 if delete or deleteGroup : i , pi , pos , refb , pattern = handleDeletions ( i , pi , pos , ref , delete , pattern ) #If no deletion, replace ref with '.', as in VCF format
 else : pattern [ pattern == refb ] = '.' #Get the list of ALTs - minus any '.'!
 uniques = np . unique ( pattern ) uniques = uniques [ np . where ( uniques != '.' ) ] #Convert bases to the number that matches the ALT
 j = 1 for u in uniques : pattern [ np . where ( pattern == u ) [ 0 ] ] = str ( j ) j += 1 #Now convert these calls to #/# (VCF format)
 calls = [ j + "/" + j if j != '.' else '.' for j in pattern ] #What if there's no variation at a variable site??
 #This can happen when sites are modified by TreeTime - see below.
 printPos = True if len ( uniques ) == 0 : #If we expect it (it was made constant by TreeTime), it's fine.
 if 'inferred_const_sites' in tree_dict and pi in tree_dict [ 'inferred_const_sites' ] : explainedErrors += 1 printPos = False #and don't output position to the VCF
 else : #If we don't expect, raise an error
 errorPositions . append ( str ( pi ) ) #Write it out - Increment positions by 1 so it's in VCF numbering
 #If no longer variable, and explained, don't write it out
 if printPos : output = [ "MTB_anc" , str ( pos ) , "." , refb , "," . join ( uniques ) , "." , "PASS" , "." , "GT" ] + calls vcfWrite . append ( "\t" . join ( output ) ) i += 1 #Note: The number of 'inferred_const_sites' passed back by TreeTime will often be longer
 #than the number of 'site that were made constant' that prints below. This is because given the site:
 # Ref   Alt     Seq
 # G     A       AANAA
 #This will be converted to 'AAAAA' and listed as an 'inferred_const_sites'. However, for VCF
 #purposes, because the site is 'variant' against the ref, it is variant, as expected, and so
 #won't be counted in the below list, which is only sites removed from the VCF.
 if 'inferred_const_sites' in tree_dict and explainedErrors != 0 : print ( "Sites that were constant except for ambiguous bases were made constant by TreeTime. This happened {} times. These sites are now excluded from the VCF." . format ( explainedErrors ) ) if len ( errorPositions ) != 0 : print ( "\n***WARNING: vcf_utils.py" "\n{} sites were found that had no alternative bases. If this data has been " "run through TreeTime and contains ambiguous bases, try calling get_tree_dict with " "var_ambigs=True to see if this clears the error." "\n\nAlternative causes:" "\n- Not all sequences in your alignment are in the tree (if you are running TreeTime via commandline " "this is most likely)" "\n- In TreeTime, can be caused by overwriting variants in tips with small branch lengths (debug)" "\n\nThese are the positions affected (numbering starts at 0):" . format ( str ( len ( errorPositions ) ) ) ) print ( "," . join ( errorPositions ) ) with open ( file_name , 'a' ) as the_file : the_file . write ( "\n" . join ( vcfWrite ) ) if file_name . endswith ( ( '.gz' , '.GZ' ) ) : import os #must temporarily remove .gz ending, or gzip won't zip it!
 os . rename ( file_name , file_name [ : - 3 ] ) call = [ "gzip" , file_name [ : - 3 ] ] os . system ( " " . join ( call ) )
Writes out a VCF-style file (which seems to be minimally handleable
    by vcftools and pyvcf) of the alignment. This is created from a dict
    in a similar format to what's created by :py:meth:`treetime.vcf_utils.read_vcf`

    Positions of variable sites are transformed to start at 1 to match
    VCF convention.

    Parameters
    ----------
     tree_dict: nested dict
        A nested dict with keys 'sequence' 'reference' and 'positions',
        as is created by :py:meth:`treetime.TreeAnc.get_tree_dict`

     file_name: str
        File to which the new VCF should be written out. File names ending with
        '.gz' will result in the VCF automatically being gzipped.
Writes out a VCF - style file ( which seems to be minimally handleable by vcftools and pyvcf ) of the alignment . This is created from a dict in a similar format to what s created by : py : meth : treetime . vcf_utils . read_vcf Positions of variable sites are transformed to start at 1 to match VCF convention . Parameters ---------- tree_dict : nested dict A nested dict with keys sequence reference and positions as is created by : py : meth : treetime . TreeAnc . get_tree_dict file_name : str File to which the new VCF should be written out . File names ending with . gz will result in the VCF automatically being gzipped .