If given an integer, return associated human EnsemblRelease for that
    Ensembl version.

    If given a string, return latest EnsemblRelease which has a reference
    of the same name.

    If given a PyEnsembl Genome, simply return it.
975904
def infer_genome(genome_object_string_or_int):
    
    if isinstance(genome_object_string_or_int, Genome):
        return genome_object_string_or_int
    if is_integer(genome_object_string_or_int):
        return cached_release(genome_object_string_or_int)
    elif is_string(genome_object_string_or_int):
        # first infer the canonical reference name, e.g. mapping hg19 -> GRCh37
        # and then get the associated PyEnsembl Genome object
        reference_name = infer_reference_name(genome_object_string_or_int)
        return genome_for_reference_name(reference_name)
    else:
        raise TypeError(
            ("Expected genome to be an int, string, or pyensembl.Genome "
                "instance, got %s : %s") % (
                str(genome_object_string_or_int),
                type(genome_object_string_or_int)))
1374981
def get_pubtator(pmid):
    
    r = get_url(PUBTATOR_TMPL.replace("PMID", pmid), timeout=10)
    if r and r.status_code == 200:
        pubtator = r.json()[0]
    else:
        log.error(
            f"Cannot access Pubtator, status: {r.status_code} url: {PUBTATOR_TMPL.replace('PMID', pmid)}"
        )
        return None

    known_types = ["CHEBI", "Chemical", "Disease", "Gene", "Species"]

    for idx, anno in enumerate(pubtator["denotations"]):
        s_match = re.match(r"(\w+):(\w+)", anno["obj"])
        c_match = re.match(r"(\w+):(\w+):(\w+)", anno["obj"])
        if c_match:
            (ctype, namespace, cid) = (
                c_match.group(1),
                c_match.group(2),
                c_match.group(3),
            )

            if ctype not in known_types:
                log.info(f"{ctype} not in known_types for Pubtator")
            if namespace not in known_types:
                log.info(f"{namespace} not in known_types for Pubtator")

            pubtator["denotations"][idx][
                "obj"
            ] = f'{pubtator_ns_convert.get(namespace, "UNKNOWN")}:{cid}'
            pubtator["denotations"][idx]["entity_type"] = pubtator_entity_convert.get(
                ctype, None
            )
            pubtator["denotations"][idx][
                "annotation_type"
            ] = pubtator_annotation_convert.get(ctype, None)
        elif s_match:
            (ctype, cid) = (s_match.group(1), s_match.group(2))

            if ctype not in known_types:
                log.info(f"{ctype} not in known_types for Pubtator")

            pubtator["denotations"][idx][
                "obj"
            ] = f'{pubtator_ns_convert.get(ctype, "UNKNOWN")}:{cid}'
            pubtator["denotations"][idx]["entity_type"] = pubtator_entity_convert.get(
                ctype, None
            )
            pubtator["denotations"][idx][
                "annotation_type"
            ] = pubtator_annotation_convert.get(ctype, None)

    annotations = {}
    for anno in pubtator["denotations"]:
        log.info(anno)
        if anno["obj"] not in annotations:
            annotations[anno["obj"]] = {"spans": [anno["span"]]}
            annotations[anno["obj"]]["entity_types"] = [anno.get("entity_type", [])]
            annotations[anno["obj"]]["annotation_types"] = [
                anno.get("annotation_type", [])
            ]

        else:
            annotations[anno["obj"]]["spans"].append(anno["span"])

    del pubtator["denotations"]
    pubtator["annotations"] = copy.deepcopy(annotations)

    return pubtator
915286
def alleleupdater(self, sample, gene, targetallele):
        
        from Bio.Seq import Seq
        from Bio.Alphabet import IUPAC
        from Bio.SeqRecord import SeqRecord
        # As there is some discrepancy with the capitalisation of terms, make sure it is consistent
        analysistype = 'rMLST' if self.analysistype.lower() == 'rmlst' else 'MLST'

        # Set the directory containing the profile and alleles
        alleledir = self.referencefilepath + analysistype + '/local' if self.pipeline else self.referencefilepath
        allelefile = glob('{}/{}.fa'.format(alleledir, gene))[0]
        # Create a string to store the last local allele number
        nextallele = str()
        # Check the local allele file to see if this allele has already been recorded
        for record in SeqIO.parse(allelefile, 'fasta'):
            # String match of the sequence to the sequence of the alleles - if they match then set the name of the
            # matching allele to this local allele
            if sample[self.analysistype].queryseq[gene] == str(record.seq):

                allelenumber = record.id.split('_')[-1]
                # Return allelenumber, percent identity (100.0%), hsp.score - 100?
                return '', '', ''
            # Record the allele number + 1; following the last record, this number will represent the next allele number
            nextallele = int(record.id.split('_')[-1]) + 1
        # Translate the nucleotide sequence to determine if there are any internal stop codons
        dnaseq = Seq(sample[self.analysistype].queryseq[gene], IUPAC.unambiguous_dna)
        protseq = str(dnaseq.translate())
        # There should be only one stop codon per sequence. In sequences with more than one stop codon, this additional
        # stop codon must be internal
        internalstop = True if protseq.count('*') > 1 else False
        if not internalstop:

            # If the alignment length is less than the subject length, then the query sequence is truncated.
            # This may not be an issue, but it must still be flagged
            truncated = True if sample[self.analysistype].alignmentlength[gene] < sample[self.analysistype].subjectlength[gene] else False
            if not truncated:
                print('full length', sample.name, gene, nextallele, targetallele, alleledir, allelefile, protseq)
                # The header will be >BACT00001_1000000
                definitionline = '{}_{} {} NT from allele {}### no internal stop codons'\
                    .format(gene, nextallele, sample[self.analysistype].mismatches[gene], sample[self.analysistype].closealleles[gene])
                # Create a sequence record using BioPython
                fasta = SeqRecord(dnaseq,
                                  # Without this, the header will be improperly formatted
                                  description='',
                                  # Use >:definitionline as the header
                                  id=definitionline)
                print(definitionline)
                # with open(sample[self.analysistype].supplementalalleles) as supplemental:
                #     # Use the SeqIO module to properly format the new sequence record
                #     SeqIO.write(fasta, supplemental, 'fasta')
                # with open(allelefile) as localalleles:
                #     # Use the SeqIO module to properly format the new sequence record
                #     SeqIO.write(fasta, localalleles, 'fasta')
            else:
                print('truncated', sample.name, gene, nextallele, targetallele, alleledir, allelefile, protseq)

        else:
            print('internal stop', sample.name, gene, nextallele, targetallele, alleledir, allelefile, protseq)

        return '', '', ''
693187
def find_seq_homologues(self, return_raw=False):
        
        
        
        # Get sequence from structure/chain
        # We could use Bio.PDB.Polypeptide?
        
        from Bio.SCOP.Raf import to_one_letter_code
        
        s = self
        seq_iter = s.get_residues()
        seq_str = ''

        for aa in seq_iter:
            if aa.resname in to_one_letter_code:
                seq_str += to_one_letter_code[aa.resname]
        
        # Use BLAST to find homologous sequences with associated
        # structures in the PDB database.
        # Perhaps include local BLAST?
        
        from Bio.Blast.NCBIWWW import qblast

        # Adapt for short query sequences if needed
        # From http://www.ncbi.nlm.nih.gov/blast/producttable.shtml#shortp
        
        if len(seq_str) < 15:
            word_size = 2
            expect = 20000
            matrix_name = 'PAM30'
            filter = None
        else:
            word_size = 3
            expect = 10.0
            matrix_name = 'BLOSUM62'
            filter = 'SEG'
             
        query_result = qblast(  "blastp", "pdb", seq_str, 
                                word_size, expect, matrix_name, filter)
        
        if return_raw:
           return query_result 
            
        # Parse BLAST result to yield results
        # PDBID : (E-Value, Identity, Positives, Gaps, Alignment)
        
        from Bio.Blast import NCBIXML
        
        blast_records = NCBIXML.read(query_result)
        
        results = []
        
        for alignment in blast_records.alignments:
            for hsps in alignment.hsps:
                id_perc = "%s/%s" %(hsps.identities, alignment.length)
                pos_perc = "%s/%s" %(hsps.positives, alignment.length)
                gaps_perc = "%s/%s" %(hsps.gaps, alignment.length)

                pdb_id = alignment.title.split('|')[3]
                e_value = hsps.expect
                
                results.append( ( pdb_id,
                                "%2.5e" %e_value,
                                id_perc, 
                                pos_perc, 
                                gaps_perc,
                                '\n'.join([hsps.query, hsps.match, hsps.sbjct]) ) )

        return results
1083947
def sequencetyper(self):
        """"""
        logging.info('Performing sequence typing')
        for sample in self.runmetadata.samples:
            if sample.general.bestassemblyfile != 'NA':
                if type(sample[self.analysistype].allelenames) == list:
                    # Initialise variables
                    header = 0
                    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
                        # Calculate the number of allele matches between each sequence type and the results
                        best_seqtype = dict()
                        for sequencetype in sample[self.analysistype].profiledata:
                            # Initialise a counter
                            count = 0
                            # Iterate through each gene for the sequence type
                            for gene, refallele in sample[self.analysistype].profiledata[sequencetype].items():
                                # Use the gene to extract the calculated allele
                                allele = self.matchdict[genome][gene]
                                # Increment the count on a match
                                if refallele == allele:
                                    count += 1
                            # Add the sequence type to the set of sequence types with the number of matches as the key
                            try:
                                best_seqtype[count].add(sequencetype)
                            except KeyError:
                                best_seqtype[count] = set()
                                best_seqtype[count].add(sequencetype)
                        # Find the highest number of matches from the dictionary
                        best = sorted(best_seqtype.items(), key=operator.itemgetter(0), reverse=True)[0][1]
                        # Deal with multiple allele matches
                        for gene in sample[self.analysistype].allelenames:
                            # Clear the appropriate count and lists
                            multiallele = list()
                            multipercent = list()
                            # 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 not in ['N', 'NA']:
                                    # 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)
                            # 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
                            try:
                                self.bestdict[genome][gene][" ".join(str(allele)
                                                                     for allele in sorted(multiallele))] = \
                                    multipercent[0]
                            except IndexError:
                                self.bestdict[genome][gene]['NA'] = 0
                            # Find the profile with the most alleles in common with the query genome
                            for sequencetype in best:
                                # 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 for E. coli and BACT000014
                                    # for Listeria. 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 sortedrefallele == 'N' and allele != 'N':
                                        # Increment the number of matches to each profile
                                        self.bestmatch[genome][sequencetype] += 1
                                    # Consider cases with multiple allele matches
                                    elif len(allele.split(' ')) > 1:
                                        # Also increment the number of matches if one of the alleles matches the
                                        # reference allele e.g. 41 16665 will match either 41 or 16665
                                        if sortedrefallele != 'N' and allele != 'N':
                                            match = False
                                            for sub_allele in allele.split(' '):
                                                if sub_allele == refallele:
                                                    match = True
                                            if match:
                                                # 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.analysistype == 'mlst':
                                            self.resultprofile[genome][sequencetype][sortedmatches][gene][
                                                list(self.bestdict[genome][gene]
                                                     .keys())[0]] = 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].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))}))
                                        sample[self.analysistype].mismatchestosequencetype = mismatches
                                        sample[self.analysistype].sequencetype = sequencetype
                                        sample[self.analysistype].matchestosequencetype = matches
                        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)
                            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'
                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'
            # Clear out the reverse_profiledata attribute - frozen sets can not be .json encoded
            try:
                delattr(sample[self.analysistype], 'reverse_profiledata')
            except AttributeError:
                pass
78281
private static ProteinSequence getUniprot(String uniProtID) throws Exception {

		AminoAcidCompoundSet set = AminoAcidCompoundSet.getAminoAcidCompoundSet();
		UniprotProxySequenceReader<AminoAcidCompound> uniprotSequence = new UniprotProxySequenceReader<AminoAcidCompound>(uniProtID,set);

		ProteinSequence seq = new ProteinSequence(uniprotSequence);

		return seq;
	}
886323
def get_bp_stats(sam, genomes, s2b, qual):
    
    for read in parse_sam(sam, qual):
        ref = read[2] # scaffold that read mapped to
        if ref not in s2b:
            continue
        genome = s2b[ref] # genome that scaffold belongs to
        refs = genomes[genome][ref][sam.name]['bp_stats']
        ref_start = int(read[3]) - 1 # position of start of alignment on reference
        ref_pos = int(read[3]) - 1 # for keeping track of reference region
        sequence = list(read[9]) # read sequence
        cigar = parse_cigar(read[5]) # parsed cigar string
        cigar_start = 0 # for keeping track of start of cigar regions
        bases = [] # bases to compare with reference
        for cigar_pos, status in cigar:
            if status == 'D': # deletion compared to reference
                refs[ref_pos - 1]['Del'].append(cigar_pos)
                for b in range(0, cigar_pos):
                    bases.append(False)
                ref_pos += cigar_pos
            else:
                cigar_stop = cigar_start + cigar_pos
                if status == 'M': # aligned to reference
                    for b in sequence[cigar_start:cigar_stop]: # bases for cigar region
                        bases.append(b)
                    ref_pos += cigar_pos
                elif status == 'I': # insertion compared to reference
                    refs[ref_pos - 1]['In'].append(cigar_pos)
                else:
                    print('# unrecognized cigar character: %s' % (status), file=sys.stderr)
                    exit()
                cigar_start += cigar_pos
        # add base to frequency at each position
        for base, position in zip(bases, list(range(ref_start, ref_start + len(bases)))):
            if base is False:
                continue
            try:
                refs[position][base.upper()] += 1
            except IndexError:
                continue
    return genomes
475557
def get_gene_info(ensembl_ids=None, hgnc_symbols=None):
    
    uniq_ensembl_ids = set(ensembl_id for ensembl_id in (ensembl_ids or []))
    uniq_hgnc_symbols = set(hgnc_symbol for hgnc_symbol in (hgnc_symbols or []))
    genes = []
    gene_data = []
    
    if uniq_ensembl_ids:
        for ensembl_id in uniq_ensembl_ids:
            for res in query_gene(ensembl_id=ensembl_id):
                gene_data.append(res) 

    elif uniq_hgnc_symbols:
        for hgnc_symbol in uniq_hgnc_symbols:
            query_res = query_gene(hgnc_symbol=hgnc_symbol)
            if query_res:
                for res in query_res:
                    gene_data.append(res)
            else:
                # If no result we add just the symbol
                gene_data.append({
                    'hgnc_symbol': hgnc_symbol,
                    'hgnc_id': None,
                    'ensembl_id': None,
                    'description': None,
                    'chrom': 'unknown',
                    'start': 0,
                    'stop': 0,
                    'hi_score': None,
                    'constraint_score': None,
                })
    for gene in gene_data:
        genes.append(Gene(
            symbol=gene ['hgnc_symbol'],
            hgnc_id=gene['hgnc_id'],
            ensembl_id=gene['ensembl_id'],
            description=gene['description'],
            chrom=gene['chrom'],
            start=gene['start'],
            stop=gene['stop'],
            location=get_cytoband_coord(gene['chrom'], gene['start']),
            hi_score=gene['hi_score'],
            constraint_score=gene['constraint_score'],
            omim_number=get_omim_number(gene['hgnc_symbol'])
            ))

    return genes
639438
def _get_ref_alt(self, var, boundary):
        

        # Get reference allele
        if var.posedit.edit.type == "ins" or var.posedit.edit.type == "dup":
            ref = ""
        else:
            # For NARefAlt and Inv
            if var.posedit.edit.ref_s is None or var.posedit.edit.ref == "":
                ref = self._fetch_bounded_seq(var, var.posedit.pos.start.base - 1,
                                              var.posedit.pos.end.base, 0, boundary)
            else:
                ref = var.posedit.edit.ref

        # Get alternative allele
        if var.posedit.edit.type == "sub" or var.posedit.edit.type == "delins" or var.posedit.edit.type == "ins":
            alt = var.posedit.edit.alt
        elif var.posedit.edit.type == "del":
            alt = ""
        elif var.posedit.edit.type == "dup":
            alt = var.posedit.edit.ref or self._fetch_bounded_seq(
                var, var.posedit.pos.start.base - 1, var.posedit.pos.end.base, 0, boundary)
        elif var.posedit.edit.type == "inv":
            alt = reverse_complement(ref)
        elif var.posedit.edit.type == "identity":
            alt = ref

        return ref, alt
704925
def getmlsthelper(referencefilepath, start, organism, update):
    """"""
    from accessoryFunctions.accessoryFunctions import GenObject
    # Initialise a set to for the organism(s) for which new alleles and profiles are desired
    organismset = set()
    # Allow for Shigella to use the Escherichia MLST profile/alleles
    organism = organism if organism != 'Shigella' else 'Escherichia'
    # As there are multiple profiles for certain organisms, this dictionary has the schemes I use as values
    organismdictionary = {'Escherichia': 'Escherichia coli#1',
                          'Shigella': 'Escherichia coli#1',
                          'Vibrio': 'Vibrio parahaemolyticus',
                          'Campylobacter': 'Campylobacter jejuni',
                          'Listeria': 'Listeria monocytogenes',
                          'Bacillus': 'Bacillus cereus',
                          'Klebsiella': 'Klebsiella pneumoniae'}
    # Allow for a genus not in the dictionary being specified
    try:
        organismset.add(organismdictionary[organism])
    except KeyError:
        # Add the organism to the set
        organismset.add(organism)
    for scheme in organismset:
        organismpath = os.path.join(referencefilepath, 'MLST', organism)
        # Find all folders (with the trailing / in the glob search) and remove the trailing /
        try:
            lastfolder = sorted(glob('{}/*/'.format(organismpath)))[-1].rstrip('/')
        except IndexError:
            lastfolder = []
        # Run the method to determine the most recent folder, and how recently it was updated
        delta, foldersize, d1 = schemedate(lastfolder)
        # Set the path/name of the folder to contain the new alleles and profile
        newfolder = '{}/{}'.format(organismpath, d1)
        if update:
            if delta.days > 7 or foldersize < 100:
                printtime('Downloading {} MLST scheme from pubmlst.org'.format(organism), start)
                # Create the object to store the argument attributes to feed to getmlst
                getmlstargs = GenObject()
                getmlstargs.species = scheme
                getmlstargs.repository_url = 'http://pubmlst.org/data/dbases.xml'
                getmlstargs.force_scheme_name = False
                getmlstargs.path = newfolder
                # Create the path to store the downloaded
                make_path(getmlstargs.path)
                getmlst.main(getmlstargs)
                # Even if there is an issue contacting the database, files are created, however, they are populated
                # with XML strings indicating that the download failed
                # Read the first character in the file
                try:
                    profilestart = open(glob('{}/*.txt'.format(newfolder))[0]).readline()
                except IndexError:
                    profilestart = []
                # If it is a <, then the download failed
                if not profilestart or profilestart[0] == '<':
                    # Delete the folder, and use the previous definitions instead
                    shutil.rmtree(newfolder)
                    newfolder = lastfolder
            # If the profile and alleles are up-to-date, set :newfolder to :lastfolder
            else:
                newfolder = lastfolder
        # If update isn't specified, don't update
        else:
            newfolder = lastfolder
            # Ensure that the profile/alleles updated successfully
            # Calculate the size of the folder by adding the sizes of all the files within the folder together
        try:
            newfoldersize = sum(os.path.getsize('{}/{}'.format(newfolder, f)) for f in os.listdir(newfolder)
                                if os.path.isfile('{}/{}'.format(newfolder, f)))
        except (OSError, TypeError):
            newfoldersize = 100
        # If the profile/allele failed, remove the folder, and use the most recent update
        if newfoldersize < 100:
            shutil.rmtree(newfolder)
            try:
                newfolder = sorted(glob('{}/*/'.format(organismpath)))[-1].rstrip('/')
            except IndexError:
                newfolder = organismpath
        # Return the name/path of the allele-containing folder
        return newfolder
184968
def _get_update(self, variant):
        
        update = {
                '$inc': {
                    'homozygote': variant.get('homozygote', 0),
                    'hemizygote': variant.get('hemizygote', 0),
                    'observations': 1
                },
                '$set': {
                    'chrom': variant.get('chrom'),
                    'start': variant.get('pos'),
                    'end': variant.get('end'),
                    'ref': variant.get('ref'),
                    'alt': variant.get('alt'),
                }
             }
        if variant.get('case_id'):
            update['$push'] = {
                                'families': {
                                '$each': [variant.get('case_id')],
                                '$slice': -50
                                }
                            }
        return update
523815
def is_insertion(self):
        
        # An insertion would appear in a VCF like C>CT, so that the
        # alternate allele starts with the reference nucleotides.
        # Since the nucleotide strings may be normalized in the constructor,
        # it's worth noting that the normalized form of this variant would be
        # ''>'T', so that 'T'.startswith('') still holds.
        return (len(self.ref) < len(self.alt)) and self.alt.startswith(self.ref)
1245469
def is_snp(self):
        ''''''
        nucleotides = {'A', 'C', 'G', 'T'}
        return len(self.REF) == 1 and self.REF in nucleotides and set(self.ALT).issubset(nucleotides)
396978
def protein_sequence(self, protein_id):
        
        if self.protein_sequences is None:
            raise ValueError(
                "No protein FASTA supplied to this Genome: %s" % self)
        return self.protein_sequences.get(protein_id)
1252562
def get_version():
    
    with open('sj.py') as source:
        for node in ast.walk(ast.parse(source.read(), 'sj.py')):
            if node.__class__.__name__ == 'Assign' and \
               node.targets[0].__class__.__name__ == 'Name' and \
               node.targets[0].id == '__version__':
                return node.value.s
524291
def get_protein_modification_language(identifier_qualified: ParserElement) -> ParserElement:
    """"""
    pmod_identifier = MatchFirst([
        identifier_qualified,
        pmod_default_ns,
        pmod_legacy_ns
    ])

    return pmod_tag + nest(
        Group(pmod_identifier)(IDENTIFIER) +
        Optional(
            WCW +
            amino_acid(PMOD_CODE) +
            Optional(
                WCW +
                ppc.integer(PMOD_POSITION)
            )
        )
    )
1298880
public DNASequence getSequence5PrimeTo3Prime() {
		String sequence = getSequenceAsString(this.getBioBegin(), this.getBioEnd(), this.getStrand());
		if (getStrand() == Strand.NEGATIVE) {
			//need to take complement of sequence because it is negative and we are returning the gene sequence from the opposite strand
			StringBuilder b = new StringBuilder(getLength());
			CompoundSet<NucleotideCompound> compoundSet = this.getCompoundSet();
			for (int i = 0; i < sequence.length(); i++) {
				String nucleotide = String.valueOf(sequence.charAt(i));
				NucleotideCompound nucleotideCompound = compoundSet.getCompoundForString(nucleotide);
				b.append(nucleotideCompound.getComplement().getShortName());
			}
			sequence = b.toString();
		}
		DNASequence dnaSequence = null;
		try {
			dnaSequence = new DNASequence(sequence.toUpperCase());
		} catch (CompoundNotFoundException e) {
			// this should not happen, the sequence is DNA originally, if it does, there's a bug somewhere
			logger.error("Could not create new DNA sequence in getSequence5PrimeTo3Prime(). Error: {}",e.getMessage());
		}
		dnaSequence.setAccession(new AccessionID(this.getAccession().getID()));
		return dnaSequence;
	}