First read the genome
966991
def main(args):

  
  sys.stderr.write("reading reference genome\n")
  ref = FASTAData(open(args.genome).read())
  sys.stderr.write("read in "+str(len(ref.keys()))+" chromosomes\n")

  """ Next make the transcriptome """
  txome = {}
  sys.stderr.write("write the transcriptome\n")
  inf = None
  if args.gpd[-3:] == '.gz':
    inf = gzip.open(args.gpd)
  else:
    inf = open(args.gpd)
  stream = GPDStream(inf)
  tof = open(args.tempdir+'/transcriptome.fa','w')
  z = 0
  for gpd in stream:
    z += 1
    if gpd.transcript_name in txome:
      sys.stderr.write("WARNING already have a transcript "+gpd.transcript_name+" ignoring line "+str(z)+" of the gpd\n")
      continue
    txome[gpd.transcript_name] = gpd.gene_name
    tof.write('>'+gpd.transcript_name+"\n"+str(gpd.get_sequence(ref))+"\n")
  tof.close()
  inf.close()
  sys.stderr.write("wrote "+str(len(txome.keys()))+" transcripts\n")

  """Build the salmon index"""
  sys.stderr.write("building a salmon index\n")
  cmd = 'salmon index -p '+str(args.numThreads)+' -t '+args.tempdir+'/transcriptome.fa -i '+args.tempdir+'/salmon_index'
  p = Popen(cmd.split())
  p.communicate()
  sys.stderr.write("finished building the index\n")

  """Use the index to quanitfy"""
  sys.stderr.write("quanitfy reads\n")
  reads = ''
  if args.rU:
    reads = '-r '+args.rU
  else:
    reads = '-1 '+args.r1+' -2 '+args.r2
  cmd = 'salmon quant -p '+str(args.numThreads)+' -i '+args.tempdir+'/salmon_index -l A '+reads+' -o '+args.tempdir+'/output_quant'
  p = Popen(cmd.split())
  p.communicate()
  sys.stderr.write("finished quanitfying\n")

  """Now parse the salmon output to add gene name"""
  salmon = {}
  with open(args.tempdir+'/output_quant/quant.sf') as inf:
    header = inf.readline()
    for line in inf:
      f = line.rstrip().split("\t")
      # by each transcript name hold a data strcture of the information
      salmon[f[0]] = {'name':f[0],'length':int(f[1]),'EffectiveLength':float(f[2]),'TPM':float(f[3]),'NumReads':float(f[4])}
  genes = {}
  for name in salmon:
    gene = txome[name]
    if gene not in genes: genes[gene] = []
    genes[gene].append(salmon[name])
  genetot = {}
  for gene in genes:
    tot = sum([x['TPM'] for x in genes[gene]])
    genetot[gene] = tot
  ordered_gene_names = sorted(genetot.keys(), key=lambda x: genetot[x],reverse=True)
      
  """Collected enough information to make output"""
  sys.stderr.write("generating output\n")
  of = sys.stdout
  if args.output != '-':
    of = open(args.output,'w')
  of.write("geneName\ttranscriptName\tlength\tEffectiveLength\ttxTPM\tNumReads\tgeneTPM\n")
  for gene in ordered_gene_names:
    txs = sorted(genes[gene],key=lambda x: x['TPM'],reverse=True)
    for tx in txs:
      of.write(gene+"\t"+tx['name']+"\t"+str(tx['length'])+"\t"+str(tx['EffectiveLength'])+"\t"+str(tx['TPM'])+"\t"+str(tx['NumReads'])+"\t"+str(genetot[gene])+"\n")
  of.close()

  # Temporary working directory step 3 of 3 - Cleanup
  if not args.specific_tempdir:
    rmtree(args.tempdir)
1251120
def reporter(self):
        
        # Create a set of all the gene names without alleles or accessions e.g. sul1_18_AY260546 becomes sul1
        genedict = dict()
        # Load the notes file to a dictionary
        notefile = os.path.join(self.targetpath, 'notes.txt')
        with open(notefile, 'r') as notes:
            for line in notes:
                # Ignore comment lines - they will break the parsing
                if line.startswith('#'):
                    continue
                # Split the line on colons e.g. stx1Aa:  Shiga toxin 1, subunit A, variant a: has three variables after
                # the split: gene(stx1Aa), description(Shiga toxin 1, subunit A, variant a), and _(\n)
                try:
                    gene, description, _ = line.split(':')
                # There are exceptions to the parsing. Some lines only have one :, while others have three. Allow for
                # these possibilities.
                except ValueError:
                    try:
                        gene, description = line.split(':')
                    except ValueError:
                        gene, description, _, _ = line.split(':')
                # Set up the description dictionary
                genedict[gene] = description.replace(', ', '_').strip()
        # Find unique gene names with the highest percent identity
        for sample in self.runmetadata.samples:
            try:
                if sample[self.analysistype].results:
                    # Initialise a dictionary to store the unique genes, and their percent identities
                    sample[self.analysistype].uniquegenes = dict()
                    for name, identity in sample[self.analysistype].results.items():
                        # Split the name of the gene from the string e.g. stx1:11:Z36899:11 yields stx1
                        if ':' in name:
                            sample[self.analysistype].delimiter = ':'
                        else:
                            sample[self.analysistype].delimiter = '_'
                        genename = name.split(sample[self.analysistype].delimiter)[0]
                        # Set the best observed percent identity for each unique gene
                        try:
                            # Pull the previous best identity from the dictionary
                            bestidentity = sample[self.analysistype].uniquegenes[genename]
                            # If the current identity is better than the old identity, save it
                            if float(identity) > float(bestidentity):
                                sample[self.analysistype].uniquegenes[genename] = float(identity)
                        # Initialise the dictionary if necessary
                        except KeyError:
                            sample[self.analysistype].uniquegenes[genename] = float(identity)
            except AttributeError:
                raise
        # Create the path in which the reports are stored
        make_path(self.reportpath)
        # Initialise strings to store the results
        data = 'Strain,Gene,Subtype/Allele,Description,Accession,PercentIdentity,FoldCoverage\n'
        with open(os.path.join(self.reportpath, self.analysistype + '.csv'), 'w') as report:
            for sample in self.runmetadata.samples:
                try:
                    if sample[self.analysistype].results:
                        # If there are many results for a sample, don't write the sample name in each line of the report
                        for name, identity in sorted(sample[self.analysistype].results.items()):
                            # Check to see which delimiter is used to separate the gene name, allele, accession, and
                            # subtype information in the header
                            if len(name.split(sample[self.analysistype].delimiter)) == 4:
                                # Split the name on the delimiter: stx2A:63:AF500190:d; gene: stx2A, allele: 63,
                                # accession: AF500190, subtype: d
                                genename, allele, accession, subtype = name.split(sample[self.analysistype].delimiter)
                            elif len(name.split(sample[self.analysistype].delimiter)) == 3:
                                # Treat samples without a subtype e.g. icaC:intercellular adhesion protein C: differently.
                                # Extract the allele as the 'subtype', and the gene name, and accession as above
                                genename, subtype, accession = name.split(sample[self.analysistype].delimiter)
                            else:
                                genename = name
                                subtype = ''
                                accession = ''
                            # Retrieve the best identity for each gene
                            percentid = sample[self.analysistype].uniquegenes[genename]
                            # If the percent identity of the current gene matches the best percent identity, add it to
                            # the report - there can be multiple occurrences of genes e.g.
                            # sul1,1,AY224185,100.00,840 and sul1,2,CP002151,100.00,927 are both included because they
                            # have the same 100% percent identity
                            if float(identity) == percentid:
                                # Treat the initial vs subsequent results for each sample slightly differently - instead
                                # of including the sample name, use an empty cell instead
                                try:
                                    description = genedict[genename]
                                except KeyError:
                                    description = 'na'
                                # Populate the results
                                data += '{samplename},{gene},{subtype},{description},{accession},{identity},{depth}\n'\
                                    .format(samplename=sample.name,
                                            gene=genename,
                                            subtype=subtype,
                                            description=description,
                                            accession=accession,
                                            identity=identity,
                                            depth=sample[self.analysistype].avgdepth[name])
                    else:
                        data += sample.name + '\n'
                except (KeyError, AttributeError):
                    data += sample.name + '\n'
            # Write the strings to the file
            report.write(data)
1008397
def generate_combined_fasta(self, genome_list, genome_dir):
        '''
        '''

        fasta = []
            
        for genome in genome_list:
            
            full_path = genome_dir + genome
            handle = open(full_path, "rU")
            print 'making combined fasta for', genome
            try:
                seq_record = SeqIO.read(handle, 'genbank')
                org_accession = seq_record.name
                
            except AssertionError,e:
                print str(e), genome
                        
            for feature in seq_record.features:
                if feature.type == 'CDS':
                    try:                                  
                        prot_accession = feature.qualifiers['protein_id'][0]
                        prot_translation = feature.qualifiers['translation'][0]
                        
                        newfast = '>' + org_accession + ',' + prot_accession + \
                         '\n' + prot_translation + '\n'
                        
                        #if newfast not in fasta:
                        fasta.append(newfast)
                            
                    except AttributeError,e:
                        print "organism %s, protein %s did not have \
                        the right attributes" % (org_accession, prot_accession)
                        print str(e)
                    
                    except KeyError,e:
                        print "organism %s, protein %s did not have \
                        the right key" % (org_accession, prot_accession)
                        print str(e)
                          
            handle.close()
            
            print "%s proteins were added" % len(fasta)

            set_fasta = set(fasta)

            print "%s unique proteins were added -- dropping redundant ones" % len(set_fasta)
                                                
        faastring = "".join(set_fasta)

        write_fasta = open('combined_fasta', 'w')
        write_fasta.write(faastring)
        write_fasta.close()

        return set_fasta
1327341
def best_genomes(self, n):
        """"""
        def key(g):
            return g.fitness

        return sorted(self.most_fit_genomes, key=key, reverse=True)[:n]
176247
def load_sequences_to_reference(self, sc=None, force_rerun=False):
        """"""
        log.info('Loading sequences to reference GEM-PRO...')
        from random import shuffle
        g_ids = [g.id for g in self.reference_gempro.functional_genes]
        shuffle(g_ids)

        def _load_sequences_to_reference_gene_sc(g_id, outdir=self.sequences_by_gene_dir,
                                                 g_to_pickle=self.gene_protein_pickles,
                                                 strain_infodict=self.strain_infodict,
                                                 orth_matrix=self.df_orthology_matrix, force_rerun=force_rerun):
            """Load orthologous strain sequences to reference Protein object, save as new pickle"""
            import ssbio.utils
            import ssbio.io
            from Bio import SeqIO
            import os.path as op

            protein_seqs_pickle_path = op.join(outdir, '{}_protein_withseqs.pckl'.format(g_id))

            if ssbio.utils.force_rerun(flag=force_rerun, outfile=protein_seqs_pickle_path):
                protein_pickle_path = g_to_pickle[g_id]
                protein_pickle = ssbio.io.load_pickle(protein_pickle_path)

                for strain, info in strain_infodict.items():
                    strain_sequences = SeqIO.index(info['genome_path'], 'fasta')
                    strain_gene_functional = info['functional_genes'][g_id]
                    if strain_gene_functional:
                        # Pull the gene ID of the strain from the orthology matrix
                        strain_gene_key = orth_matrix.at[g_id, strain]
                        new_id = '{}_{}'.format(g_id, strain)
                        if protein_pickle.sequences.has_id(new_id):
                            continue
                        protein_pickle.load_manual_sequence(seq=strain_sequences[strain_gene_key],
                                                            ident=new_id,
                                                            set_as_representative=False)
                protein_pickle.save_pickle(outfile=protein_seqs_pickle_path)

            return g_id, protein_seqs_pickle_path

        if sc:
            genes_rdd = sc.parallelize(g_ids)
            result = genes_rdd.map(_load_sequences_to_reference_gene_sc).collect()
        else:
            result = []
            for g in tqdm(g_ids):
                result.append(self._load_sequences_to_reference_gene(g, force_rerun))

        log.info('Storing paths to new Protein objects in self.gene_protein_pickles...')
        updated = []
        for g_id, protein_pickle in result:
            self.gene_protein_pickles[g_id] = protein_pickle
            updated.append(g_id)
        not_updated = set(list(self.gene_protein_pickles.keys())).difference(updated)
        log.info('No change to {} genes, removing from gene_protein_pickles'.format(len(not_updated)))
        log.debug(not_updated)
        for rem in not_updated:
            del self.gene_protein_pickles[rem]
1159976
def liftover_to_genome(pass_pos, gtf):
    """"""

    fixed_pos = []
    for pos in pass_pos:
        if pos["chrom"] not in gtf:
            continue
        db_pos = gtf[pos["chrom"]][0]
        mut = _parse_mut(pos["sv"])
        print([db_pos, pos])
        if db_pos[3] == "+":
            pos['pre_pos'] = db_pos[1] + pos["pre_pos"] + 1
        else:
            pos['pre_pos'] = db_pos[2] - (pos["pre_pos"] - 1)
        pos['chrom'] = db_pos[0]
        pos['nt'] = list(mut[0])
        fixed_pos.append(pos)

    _print_header(fixed_pos)
    for pos in fixed_pos:
        print_vcf(pos)
1019750
def apply(self, read_tuple_name, read_tuple_id=None, synchronize_widths=True):
        
        parts = read_tuple_name.split("__")
        parts[0] = self._fill_right(parts[0], "-", self.prefix_width)
        if read_tuple_id is not None:
            parts[1] = "{:x}".format(read_tuple_id)
        parts[1] = self._fill_left(parts[1], "0", self.read_tuple_id_width)

        if synchronize_widths:
            new_segments = []
            segments = parts[2][1:-1].split("),(")
            for segment in segments:
                values = segment.split(",")
                values[0] = values[0].zfill(self.genome_id_width)
                values[1] = values[1].zfill(self.chr_id_width)
                values[3] = values[3].zfill(self.coor_width)
                values[4] = values[4].zfill(self.coor_width)
                new_segments.append("(" + ",".join(values) + ")")
            parts[2] = ",".join(new_segments)

        return "__".join(parts)
466433
def load(self, read_tuple_name):
        
        self.prefix_width = 0
        self.read_tuple_id_width = 0
        self.genome_id_width = 0
        self.chr_id_width = 0
        self.coor_width = 0

        parts = read_tuple_name.split("__")
        self.prefix_width = len(parts[0])
        self.read_tuple_id_width = len(parts[1])

        segments = parts[2][1:-1].split("),(")
        for segment in segments:
            int_widths = list(map(len, segment.split(",")))
            self.genome_id_width = max(self.genome_id_width, int_widths[0])
            self.chr_id_width = max(self.chr_id_width, int_widths[1])
            self.coor_width = max(self.coor_width, int_widths[2], int_widths[3])
291260
def _guess_seqtype_from_file(handle):
    ""
    if isinstance(handle, basestring):
        handle = StringIO(handle)

    for line in handle:
        if not line.strip():
            continue
        if line.lstrip().split()[0] in ('LOCUS', 'FEATURES', 'source', 'CDS',
                                        'gene'):
            return 'genbank'
        if len(line) > 2 and line[:3] in ('ID ', 'FT '):
            return 'embl'
        if line.startswith('>'):
            return 'fasta'
    handle.seek(0)
    import string
    from Bio.Data import IUPACData as iupac
    all_input_letters = set(handle.read().lower())
    all_valid = set(string.digits)
    all_valid.update(set(iupac.protein_letters.lower()))
    all_valid.update(set(iupac.unambiguous_dna_letters.lower()))
    all_valid.update(set('- \n'))
    if all_valid.issuperset(all_input_letters):
        return 'fasta'

    raise ValueError("Failed to guess format for input")
507872
def _prepare_data(self):
        
        self._extract_genes()
        self._extract_total_number_of_chars()
        self._extract_number_of_taxa()
        self._extract_reading_frames()

        Data = namedtuple('Data', ['gene_codes', 'number_taxa', 'number_chars',
                                   'seq_records', 'gene_codes_and_lengths',
                                   'reading_frames'])
        self.data = Data(self.gene_codes, self.number_taxa, self.number_chars,
                         self.seq_records, self._gene_codes_and_lengths,
                         self.reading_frames)
632753
def get_start_stops(transcript_sequence, start_codons=None, stop_codons=None):
        
        transcript_sequence = transcript_sequence.upper()  # for comparison with codons below
        if not start_codons:
            start_codons = ['ATG']
        if not stop_codons:
            stop_codons = ['TAA', 'TAG', 'TGA']

        seq_frames = {1: {'starts': [], 'stops': []},
                      2: {'starts': [], 'stops': []},
                      3: {'starts': [], 'stops': []}}

        for codons, positions in ((start_codons, 'starts'),
                                  (stop_codons, 'stops')):
            if len(codons) > 1:
                pat = re.compile('|'.join(codons))
            else:
                pat = re.compile(codons[0])

            for m in re.finditer(pat, transcript_sequence):
                # Increment position by 1, Frame 1 starts at position 1 not 0
                start = m.start() + 1
                rem = start % 3
                if rem == 1:  # frame 1
                    seq_frames[1][positions].append(start)
                elif rem == 2:  # frame 2
                    seq_frames[2][positions].append(start)
                elif rem == 0:  # frame 3
                    seq_frames[3][positions].append(start)
        return seq_frames
510374
def prepare(args):
    
    p = OptionParser(prepare.__doc__ + FastqNamings)
    p.add_option("-K", default=51, type="int", help="K-mer size")
    p.set_cpus(cpus=32)
    opts, args = p.parse_args(args)

    if len(args) < 2:
        sys.exit(not p.print_help())

    genomesize = float(args[0]) / 1000
    fnames = args[1:]
    for x in fnames:
        assert op.exists(x), "File `{0}` not found.".format(x)

    s = comment_banner("Meraculous params file") + "\n"
    s += comment_banner("Basic parameters") + "\n"
    s += "# Describe the libraries ( one line per library )\n"
    s += "# " + " ".join(header.split()) + "\n"

    libs = get_libs(fnames)
    lib_seqs = []
    rank = 0
    for lib, fs in libs:
        size = lib.size
        if size == 0:
            continue
        rank += 1
        library_name = lib.library_name
        name = library_name.replace("-", "")
        wildcard = "{0}*.1.*,{0}*.2.*".format(library_name)
        rl = max(readlen([x]) for x in fs)
        lib_seq = lib.get_lib_seq(wildcard, name, rl, rank)
        lib_seqs.append(lib_seq)

    s += "\n" + "\n".join(load_csv(None, lib_seqs, sep=" ")) + "\n"
    params = [("genome_size", genomesize),
              ("is_diploid", 0),
              ("mer_size", opts.K),
              ("num_prefix_blocks", 1),
              ("no_read_validation", 0),
              ("local_num_procs", opts.cpus)]
    s += "\n" + "\n".join(load_csv(None, params, sep=" ")) + "\n"

    cfgfile = "meraculous.config"
    write_file(cfgfile, s, tee=True)

    s = "~/export/meraculous/bin/run_meraculous.sh -c {0}"\
                .format(cfgfile)
    runsh = "run.sh"
    write_file(runsh, s)
1298476
def deleteGenome(species, name) :
    """"""

    printf('deleting genome (%s, %s)...' % (species, name))

    conf.db.beginTransaction()
    objs = []
    allGood = True
    try :
        genome = Genome_Raba(name = name, species = species.lower())
        objs.append(genome)
        pBar = ProgressBar(label = 'preparing')
        for typ in (Chromosome_Raba, Gene_Raba, Transcript_Raba, Exon_Raba, Protein_Raba) :
            pBar.update()
            f = RabaQuery(typ, namespace = genome._raba_namespace)
            f.addFilter({'genome' : genome})
            for e in f.iterRun() :
                objs.append(e)
        pBar.close()
        
        pBar = ProgressBar(nbEpochs = len(objs), label = 'deleting objects')
        for e in objs :
            pBar.update()
            e.delete()
        pBar.close()
        
    except KeyError as e :
        #~ printf("\tWARNING, couldn't remove genome form db, maybe it's not there: ", e)
        raise KeyError("\tWARNING, couldn't remove genome form db, maybe it's not there: ", e)
        allGood = False
    printf('\tdeleting folder')
    try :
        shutil.rmtree(conf.getGenomeSequencePath(species, name))
    except OSError as e:
        #~ printf('\tWARNING, Unable to delete folder: ', e)
        OSError('\tWARNING, Unable to delete folder: ', e)
        allGood = False
        
    conf.db.endTransaction()
    return allGood
332779
def read_sequencing(directory):
    '''

    '''
    dirfiles = os.listdir(directory)
    seq_exts = ['.seq', '.abi', '.ab1']
    # Exclude files that aren't sequencing results
    seq_paths = [x for x in dirfiles if os.path.splitext(x)[1] in seq_exts]
    paths = [os.path.join(directory, x) for x in seq_paths]
    sequences = [read_dna(x) for x in paths]

    return sequences
667708
def sequencetyper(self):
        """"""
        for sample in self.metadata:
            if sample.general.bestassemblyfile != 'NA':
                if type(sample[self.analysistype].allelenames) == list:
                    # Initialise variables
                    header = 0
                    # Iterate through the genomes
                    # for sample in self.metadata:
                    genome = sample.name
                    # Initialise self.bestmatch[genome] with an int that will eventually be replaced by the # of matches
                    self.bestmatch[genome] = defaultdict(int)
                    if sample[self.analysistype].profile != 'NA':
                        # Create the profiledata variable to avoid writing self.profiledata[self.analysistype]
                        profiledata = sample[self.analysistype].profiledata
                        # For each gene in plusdict[genome]
                        for gene in sample[self.analysistype].allelenames:
                            # Clear the appropriate count and lists
                            multiallele = []
                            multipercent = []
                            # Go through the alleles in plusdict
                            for allele in self.plusdict[genome][gene]:
                                percentid = list(self.plusdict[genome][gene][allele].keys())[0]
                                # "N" alleles screw up the allele splitter function
                                if allele != "N":
                                    # Use the alleleSplitter function to get the allele number
                                    # allelenumber, alleleprenumber = allelesplitter(allele)
                                    # Append as appropriate - alleleNumber is treated as an integer for proper sorting
                                    multiallele.append(int(allele))
                                    multipercent.append(percentid)
                                # If the allele is "N"
                                else:
                                    # Append "N" and a percent identity of 0
                                    multiallele.append("N")
                                    multipercent.append(0)
                                if not multiallele:
                                    multiallele.append("N")
                                    multipercent.append(0)
                            # if self.analysistype == 'rmlst':
                            #     # For whatever reason, the rMLST profile scheme treat multiple allele hits as 'N's.
                            #     multiallele = multiallele if len(multiallele) == 1 else ['N']
                            #     if multipercent:
                            #         multipercent = multipercent if len(multiallele) == 1 else [0, 0]
                            #     else:
                            #         multipercent = [0]
                            # Populate self.bestdict with genome, gene, alleles joined with a space (this was made like
                            # this because allele is a list generated by the .iteritems() above
                            self.bestdict[genome][gene][" ".join(str(allele)
                                                                 for allele in sorted(multiallele))] = multipercent[0]
                            # Find the profile with the most alleles in common with the query genome
                            for sequencetype in profiledata:
                                # The number of genes in the analysis
                                header = len(profiledata[sequencetype])
                                # refallele is the allele number of the sequence type
                                refallele = profiledata[sequencetype][gene]
                                # If there are multiple allele matches for a gene in the reference profile e.g. 10 692
                                if len(refallele.split(" ")) > 1:
                                    # Map the split (on a space) alleles as integers - if they are treated as integers,
                                    # the alleles will sort properly
                                    intrefallele = map(int, refallele.split(" "))
                                    # Create a string of the joined, sorted alleles
                                    sortedrefallele = " ".join(str(allele) for allele in sorted(intrefallele))
                                else:
                                    # Use the reference allele as the sortedRefAllele
                                    sortedrefallele = refallele
                                for allele, percentid in self.bestdict[genome][gene].items():
                                    # If the allele in the query genome matches the allele in the reference profile, add
                                    # the result to the bestmatch dictionary. Genes with multiple alleles were sorted
                                    # the same, strings with multiple alleles will match: 10 692 will never be 692 10
                                    if allele == sortedrefallele and float(percentid) == 100.00:
                                        # Increment the number of matches to each profile
                                        self.bestmatch[genome][sequencetype] += 1
                                    # Special handling of BACT000060 and BACT000065 genes. When the reference profile
                                    # has an allele of 'N', and the query allele doesn't, set the allele to 'N', and
                                    # count it as a match
                                    elif gene == 'BACT000060' or gene == 'BACT000065':
                                        if sortedrefallele == 'N' and allele != 'N':
                                            # Increment the number of matches to each profile
                                            self.bestmatch[genome][sequencetype] += 1
                                    elif allele == sortedrefallele and sortedrefallele == 'N':
                                        # Increment the number of matches to each profile
                                        self.bestmatch[genome][sequencetype] += 1
                        # Get the best number of matches
                        # From: https://stackoverflow.com/questions/613183/sort-a-python-dictionary-by-value
                        try:
                            sortedmatches = sorted(self.bestmatch[genome].items(), key=operator.itemgetter(1),
                                                   reverse=True)[0][1]
                        # If there are no matches, set :sortedmatches to zero
                        except IndexError:
                            sortedmatches = 0
                        # Otherwise, the query profile matches the reference profile
                        if int(sortedmatches) == header:
                            # Iterate through best match
                            for sequencetype, matches in self.bestmatch[genome].items():
                                if matches == sortedmatches:
                                    for gene in profiledata[sequencetype]:
                                        # Populate resultProfile with the genome, best match to profile, # of matches
                                        # to the profile, gene, query allele(s), reference allele(s), and % identity
                                        self.resultprofile[genome][sequencetype][sortedmatches][gene][
                                            list(self.bestdict[genome][gene]
                                                .keys())[0]] = str(list(self.bestdict[genome][gene].values())[0])
                                    sample[self.analysistype].sequencetype = sequencetype
                                    sample[self.analysistype].matchestosequencetype = matches
                        # If there are fewer matches than the total number of genes in the typing scheme
                        elif 0 < int(sortedmatches) < header:
                            mismatches = []
                            # Iterate through the sequence types and the number of matches in bestDict for each genome
                            for sequencetype, matches in self.bestmatch[genome].items():
                                # If the number of matches for a profile matches the best number of matches
                                if matches == sortedmatches:
                                    # Iterate through the gene in the analysis
                                    for gene in profiledata[sequencetype]:
                                        # Get the reference allele as above
                                        refallele = profiledata[sequencetype][gene]
                                        # As above get the reference allele split and ordered as necessary
                                        if len(refallele.split(" ")) > 1:
                                            intrefallele = map(int, refallele.split(" "))
                                            sortedrefallele = " ".join(str(allele) for allele in sorted(intrefallele))
                                        else:
                                            sortedrefallele = refallele
                                        # Populate self.mlstseqtype with the genome, best match to profile, # of matches
                                        # to the profile, gene, query allele(s), reference allele(s), and % identity
                                        if self.updateprofile:
                                            self.mlstseqtype[genome][sequencetype][sortedmatches][gene][
                                                str(list(self.bestdict[genome][gene]
                                                    .keys())[0])][sortedrefallele] = str(list(self.bestdict[genome][gene]
                                                                                        .values())[0])
                                        else:
                                            self.resultprofile[genome][sequencetype][sortedmatches][gene][
                                                list(self.bestdict[genome][gene].keys())[0]] \
                                                = str(list(self.bestdict[genome][gene])[0])
                                                #= str(list(self.bestdict[genome][gene].values())[0])
                                            if sortedrefallele != list(self.bestdict[sample.name][gene].keys())[0]:
                                                mismatches.append(
                                                    ({gene: ('{} ({})'.format(list(self.bestdict[sample.name][gene]
                                                                              .keys())[0], sortedrefallele))}))
                                        if not self.updateprofile or self.analysistype == 'mlst':
                                            self.resultprofile[genome][sequencetype][sortedmatches][gene][
                                                list(self.bestdict[genome][gene]
                                                    .keys())[0]] = str(list(self.bestdict[genome][gene].values())[0])
                                            sample[self.analysistype].mismatchestosequencetype = mismatches
                                            sample[self.analysistype].sequencetype = sequencetype
                                            sample[self.analysistype].matchestosequencetype = matches
                            # Add the new profile to the profile file (if the option is enabled)
                            if self.updateprofile and self.analysistype != 'mlst':
                                self.reprofiler(int(header), genome, sample)
                        elif sortedmatches == 0:
                            for gene in sample[self.analysistype].allelenames:
                                # Populate the results profile with negative values for sequence type and sorted matches
                                self.resultprofile[genome]['NA'][sortedmatches][gene]['NA'] = 0
                            # Add the new profile to the profile file (if the option is enabled)
                            if self.updateprofile:
                                self.reprofiler(int(header), genome, sample)
                            sample[self.analysistype].sequencetype = 'NA'
                            sample[self.analysistype].matchestosequencetype = 'NA'
                            sample[self.analysistype].mismatchestosequencetype = 'NA'
                        else:
                            sample[self.analysistype].matchestosequencetype = 'NA'
                            sample[self.analysistype].mismatchestosequencetype = 'NA'
                            sample[self.analysistype].sequencetype = 'NA'
                        dotter()
                else:
                    sample[self.analysistype].matchestosequencetype = 'NA'
                    sample[self.analysistype].mismatchestosequencetype = 'NA'
                    sample[self.analysistype].sequencetype = 'NA'

            else:
                sample[self.analysistype].matchestosequencetype = 'NA'
                sample[self.analysistype].mismatchestosequencetype = 'NA'
                sample[self.analysistype].sequencetype = 'NA'
456724
def seq(self):
        """"""

        if self.sequence_file:
            log.debug('{}: reading sequence from sequence file {}'.format(self.id, self.sequence_path))
            tmp_sr = SeqIO.read(self.sequence_path, 'fasta')
            return tmp_sr.seq

        elif self.metadata_file:
            log.debug('{}: reading sequence from metadata file {}'.format(self.id, self.metadata_path))
            tmp_sr = SeqIO.read(self.metadata_path, 'uniprot-xml')
            return tmp_sr.seq

        else:
            if not self._seq:
                log.debug('{}: no sequence stored in memory'.format(self.id))
            else:
                log.debug('{}: reading sequence from memory'.format(self.id))

            return self._seq
1164728
def create(genome, config):
        
        genome_config = config.genome_config
        required = required_for_output(genome_config.input_keys, genome_config.output_keys, genome.connections)

        # Gather inputs and expressed connections.
        node_inputs = {}
        for cg in itervalues(genome.connections):
            if not cg.enabled:
                continue

            i, o = cg.key
            if o not in required and i not in required:
                continue

            if o not in node_inputs:
                node_inputs[o] = [(i, cg.weight)]
            else:
                node_inputs[o].append((i, cg.weight))

        neurons = {}
        for node_key in required:
            ng = genome.nodes[node_key]
            inputs = node_inputs.get(node_key, [])
            neurons[node_key] = IZNeuron(ng.bias, ng.a, ng.b, ng.c, ng.d, inputs)

        genome_config = config.genome_config
        return IZNN(neurons, genome_config.input_keys, genome_config.output_keys)