def METHOD_NAME(args):
    ## read tree and data, if reading data fails, return with error code
    tree = Phylo.read(args.tree, 'newick')
    #check whether VCF - by existance of VCF reference file
    is_vcf = args.vcf_aa_reference is not None
    ## check file format and read in sequences
    node_data = read_node_data(args.mutations, args.tree)
    if node_data is None:
        print("ERROR: could not read mutation data "+ ("(incl sequences)" if not is_vcf else ""))
        return 1
    root_node = tree.root.name
    #if VCF, read in the reference seq for each gene, put on root
    if(is_vcf):
        node_data["nodes"][root_node]['aa_sequences'] = {}
        with open(args.vcf_aa_reference, encoding='utf-8') as handle:
            for record in SeqIO.parse(handle, "fasta"):
                if record.id==args.gene:
                    #'root' may not be same as 'reference', so apply any mutations at root here!
                    node_data["nodes"][root_node]['aa_sequences'][record.id] = get_sequence(str(record.seq), node_data["nodes"][root_node]["aa_muts"][record.id])
    # check that root node has sequences for each requested gene
    if "aa_sequences" not in node_data["nodes"][root_node]:
        print("ERROR: ancestral sequences are not provided")
        return 1
    if not args.gene in node_data["nodes"][root_node]['aa_sequences']:
        print("ERROR: ancestral sequences missing for gene",args.gene)
        return 1
    # gather all reconstructed sequences
    sequences = dict()
    is_terminal = {}
    sequences[root_node] = node_data["nodes"][root_node]['aa_sequences'][args.gene]
    is_terminal[root_node] = False
    for node in tree.get_nonterminals(order='preorder'):
        parent_sequence = sequences[node.name]
        for child in node:
            sequences[child.name] = get_sequence(parent_sequence, node_data["nodes"][child.name]["aa_muts"][args.gene])
            is_terminal[child.name] = child.is_terminal()
    # write alignments to file
    seqs = [SeqRecord.SeqRecord(seq=Seq.Seq(sequences[strain]), id=strain, name=strain, description='')
            for strain in sequences if is_terminal[strain] or args.internal_nodes]
    SeqIO.write(seqs, args.output, 'fasta')