Finds the closest reference genome to the profile of interest

def referencegenomefinder(self):
        """
        Finds the closest reference genome to the profile of interest
        """
        # Initialise dictionaries
        referencematch = defaultdict(make_dict)
        referencehits = defaultdict(make_dict)
        # Set the name of the reference profile file
        referencegenomeprofile = '{}rMLST_referenceprofile.json'.format(self.referenceprofilepath)
        # Open the reference profile and load the profile into memory
        with open(referencegenomeprofile) as referencefile:
            referencetypes = json.load(referencefile)
        # Iterate through the samples
        for sample in self.metadata:
            if sample[self.analysistype].reportdir != 'NA':
                # Iterate through the reference genomes in the profile
                for genome in referencetypes:
                    # Initialise the number of identical alleles between the assembly of interest and the
                    # reference genome to 0
                    referencehits[sample.name][genome] = 0
                    # Iterate through all the genes in the analysis
                    for gene in self.bestdict[sample.name]:
                        # If the alleles match between the assembly of interest and the reference genome, increment
                        # the number of matches by 1
                        if list(self.bestdict[sample.name][gene].keys())[0] == referencetypes[genome][gene]:
                            referencematch[sample.name][genome][gene] = 1
                            referencehits[sample.name][genome] += 1
                        else:
                            referencematch[sample.name][genome][gene] = 0
        #
        for sample in self.metadata:
            if sample[self.analysistype].reportdir != 'NA':
                # Get the best number of matches
                # From: https://stackoverflow.com/questions/613183/sort-a-python-dictionary-by-value
                try:
                    matches = sorted(referencehits[sample.name].items(),
                                     key=operator.itemgetter(1), reverse=True)
                    most_matches = matches[0][1]
                    i = 0
                    match_list = list()
                    while matches[i][1] == most_matches:
                        match_list.append(matches[i])
                        i += 1
                    sorted_list = sorted(match_list)
                    sortedmatches = sorted_list[0]
                except IndexError:
                    sortedmatches = (0, 0)
                # If there are fewer matches than the total number of genes in the typing scheme
                if 0 < int(sortedmatches[1]) < len(sample[self.analysistype].allelenames):
                    mismatches = []
                    # Iterate through the gene in the analysis
                    for gene, allele in referencetypes[sortedmatches[0]].items():
                        # Populate :self.referencegenome with the genome name, best reference match, number of matches,
                        # gene, query allele(s), and percent identity
                        percentidentity = '{:.2f}'.format(list(self.bestdict[sample.name][gene].values())[0])
                        self.referencegenome[sample.name][sortedmatches[0]][sortedmatches[1]][gene][list(self.bestdict[
                            sample.name][gene].keys())[0]] = percentidentity
                        if list(self.bestdict[sample.name][gene].keys())[0] != allele:
                            sample[self.analysistype].referencegenome = sortedmatches[0]
                            sample.general.referencegenus = sortedmatches[0].split('_')[0]
                            sample[self.analysistype].referencegenomepath = '{}{}.fa' \
                                .format(self.referenceprofilepath, sortedmatches[0])
                            sample[self.analysistype].matchestoreferencegenome = sortedmatches[1]
                            mismatches.append(({gene: ('{} ({})'.format(list(self.bestdict[sample.name][gene]
                                                                        .keys())[0], allele))}))
                        # sample[self.analysistype].mismatchestoreferencegenome = sorted(mismatches)
                        sample[self.analysistype].mismatchestoreferencegenome = mismatches
                elif sortedmatches == 0:
                    for gene in sample[self.analysistype].allelenames:
                        # Populate the profile of results with 'negative' values for sequence type and sorted matches
                        self.referencegenome[sample.name][sortedmatches[0]][0][gene]['NA'] = 0
                        sample[self.analysistype].referencegenome = 'NA'
                        sample.general.referencegenus = 'NA'
                        sample[self.analysistype].referencegenomepath = 'NA'
                        sample[self.analysistype].matchestoreferencegenome = 0
                        sample[self.analysistype].mismatchestoreferencegenome = [0]
                # Otherwise, the query profile matches the reference profile
                else:
                    for gene in referencetypes[sortedmatches[0]]:
                        # Populate self.referencegenome as above
                        self.referencegenome[sample.name][sortedmatches[0]][sortedmatches[1]][gene][list(self.bestdict[
                            sample.name][gene].keys())[0]] = '{:.2f}'.format(list(self.bestdict[
                                                                                sample.name][gene].values())[0])
                        sample[self.analysistype].referencegenome = sortedmatches[0]
                        sample[self.analysistype].referencegenomepath = '{}{}.fa' \
                            .format(self.referenceprofilepath, sortedmatches[0])
                        sample.general.referencegenus = sortedmatches[0].split('_')[0]
                        sample[self.analysistype].matchestoreferencegenome = sortedmatches[1]
                        sample[self.analysistype].mismatchestoreferencegenome = [0]
        # Print the results to file
        make_path(self.reportpath)
        with open('{}referencegenomes.csv'.format(self.reportpath), 'w') as referencegenomereport:
            row = 'Strain,referencegenome,numberofmatches\n'
            for sample in self.metadata:
                if sample[self.analysistype].reportdir != 'NA':
                    row += '{},{},{}\n'.format(sample.name, sample[self.analysistype].referencegenome,
                                               sample[self.analysistype].matchestoreferencegenome)
            referencegenomereport.write(row)
            dotter()