Download the allele files
548991
def download_threads(self, url):
        
        # Set the name of the allele file - split the gene name from the URL
        output_file = os.path.join(self.output_path, '{}.tfa'.format(os.path.split(url)[-1]))
        # Check to see whether the file already exists, and if it is unusually small
        size = 0
        try:
            stats = os.stat(output_file)
            size = stats.st_size
        except FileNotFoundError:
            pass
        # If the file doesn't exist, or is truncated, proceed with the download
        if not os.path.isfile(output_file) or size <= 100:
            # Create a new session
            session = OAuth1Session(self.consumer_key,
                                    self.consumer_secret,
                                    access_token=self.session_token,
                                    access_token_secret=self.session_secret)
            # The allele file on the server is called alleles_fasta. Update the URL appropriately
            r = session.get(url + '/alleles_fasta')
            if r.status_code == 200 or r.status_code == 201:
                if re.search('json', r.headers['content-type'], flags=0):
                    decoded = r.json()
                else:
                    decoded = r.text
                # Write the allele to disk
                with open(output_file, 'w') as allele:
                    allele.write(decoded)
1359448
def alleleupdater1(self, sample, gene, targetallele):
        
        from io import StringIO
        from Bio.Blast import NCBIXML
        from Bio.SeqRecord import SeqRecord
        from Bio.Seq import Seq
        from Bio.Alphabet import generic_dna
        # Find the name of the allele file based on whether the gene name is found in the file name
        genefile = [x for x in sample[self.analysistype].alleles if gene in x][0]
        # Remove the extension of the gene file for use in makeblastdb and blastn
        genefilenoext = genefile.split(".")[0]
        # Create the blast databases every time the method is called
        subprocess.call(shlex.split('makeblastdb -in {} -parse_seqids -max_file_sz 2GB -dbtype nucl -out {}'
                                    .format(genefile, genefilenoext)), stdout=self.fnull, stderr=self.fnull)
        # Re-perform BLAST analyses, but using an XML output that stores the alignment
        blastn = NcbiblastnCommandline(query=sample.general.bestassemblyfile, db=genefilenoext, evalue=0.1, outfmt=5)
        # Note that there is no output file specified -  the search results are currently stored in stdout
        stdout, stderr = blastn()
        if stdout.find('Hsp') != -1:
            # Map the blast record to memory
            blast_handle = StringIO(stdout)
            # Open record from memory-mapped file
            records = NCBIXML.parse(blast_handle)
            # Iterate through the records in the blast results
            for record in records:  # This process is just to retrieve HSPs from xml files
                for alignment in record.alignments:
                    for hsp in alignment.hsps:
                        # Extract the allele name from the blast result
                        allele = str(alignment.accession.split("_")[-1]) if "_" in alignment.accession else \
                            str(alignment.accession.split("-")[-1])
                        # The point of this blast is to extract the sequence of the allele
                        # If the allele name matches the target allele (the closest match e.g. BACT000063_20 -> 20)
                        if allele == str(targetallele):
                            # 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 if self.pipeline else \
                                self.referencefilepath
                            # The name of the supplemental allele file (without and with the .fa extension)
                            allelefilenoext = '{}/OLC_{}_alleles'.format(alleledir, analysistype) if self.pipeline \
                                else '{}/rMLST_combined'.format(alleledir)
                            allelefile = allelefilenoext + '.fa' if self.pipeline else allelefilenoext + '.fasta'
                            # Create the file if it doesn't exist
                            open(allelefile, 'a').close()
                            # Create a list of all the blast database files in the folder
                            dbfiles = glob('{}.n*'.format(allelefilenoext))
                            # Remove the database files
                            map(lambda y: os.remove(y), dbfiles)
                            # Create the necessary blast database files
                            subprocess.call(
                                shlex.split('makeblastdb -in {} -parse_seqids -max_file_sz 2GB -dbtype nucl -out {}'
                                            .format(allelefile, allelefilenoext)), stdout=self.fnull, stderr=self.fnull)

                            # Perform BLAST analysis using the supplemental allele file as the database
                            nestedblastn = NcbiblastnCommandline(db=allelefilenoext, evalue=0.1, outfmt=5)
                            # There is no output file specified; the search results are currently stored in stdout
                            # Additionally, the sequence from the previous BLAST query is used as stdin in this BLAST
                            import Bio.Application
                            try:
                                nestedstdout, nestedstderr = nestedblastn(stdin=hsp.query)
                            except Bio.Application.ApplicationError:
                                nestedstdout = ''
                                pass
                            # Search stdout for matches - if the term Hsp appears (the .find function will NOT
                            # return -1), a match has been found, and stdout is written to file
                            if nestedstdout.find('Hsp') != -1:
                                nested_blast_handle = StringIO(nestedstdout)
                                # Open record from memory-mapped file
                                nestedrecords = NCBIXML.parse(nested_blast_handle)
                                # Initialise variables
                                previousallele = ''
                                bitscore = ''
                                # Iterate through the records
                                for nestedrecord in nestedrecords:
                                    for nestedalignment in nestedrecord.alignments:
                                        for nestedhsp in nestedalignment.hsps:
                                            # Calculate the percent identity
                                            percentidentity = float("%.2f" % float(float(nestedhsp.identities) /
                                                                                   float(nestedalignment.length) * 100))
                                            # Only set the previous allele and the bitscore if the % identity is 100%
                                            if percentidentity == 100:
                                                # Get the allele number
                                                previousallele = nestedalignment.accession.split("_")[-1]
                                                bitscore = nestedhsp.score

                                # If this allele has already been found in the supplemental file, return these data
                                if previousallele:
                                    return gene, previousallele, 100.0, bitscore
                                # Otherwise, get the allele sequence into the supplemental file
                                else:
                                    # Initialise variables
                                    allelelist = []
                                    allelenumber = 1000000
                                    # Open the allele file to append
                                    with open(allelefile, 'ab+') as supplemental:
                                        for line in supplemental:
                                            # As there are multiple genes (e.g. BACT000001, BACT000063, etc.) check to
                                            # see if the gene of interest is in the line
                                            if gene in line:
                                                # Append the header to the list
                                                allelelist.append(line)
                                        # Find the last allele number from the header
                                        try:
                                            allelenumber = int(allelelist[-1].split("_")[1]) + 1
                                        # If there are no alleles for the gene of interest, then pass
                                        except IndexError:
                                            pass
                                        # Puts the HSP in the correct order -  hits to the negative strand will be
                                        # reversed compared to what we're looking for
                                        allelesequence = Seq(hsp.query, generic_dna)
                                        if hsp.sbjct_start < hsp.sbjct_end:
                                            end = hsp.sbjct_end
                                        else:
                                            end = hsp.sbjct_start
                                            allelesequence = allelesequence.reverse_complement()
                                        # Screen out hits that are shorter than the targets
                                        # Keeping this format though this statement could be re-written more efficiently
                                        if end < alignment.length:
                                            pass
                                        # The header will be >BACT00001_1000000
                                        definitionline = '{}_{}'.format(gene, allelenumber)
                                        # Create a sequence record using BioPython
                                        fasta = SeqRecord(allelesequence,
                                                          # Without this, the header will be improperly formatted
                                                          description='',
                                                          # Use >:definitionline as the header
                                                          id=definitionline)
                                        # Use the SeqIO module to properly format the new sequence record
                                        SeqIO.write(fasta, supplemental, "fasta")
                                        # Return the necessary values
                                        return gene, allelenumber, 100.0, hsp.score
                            # If there are hits in the supplemental allele file
                            else:
                                # Initialise variables
                                allelelist = []
                                allelenumber = 1000000
                                # Open the allele file to find the last allele associated with the gene of interest
                                with open(allelefile, 'a+') as supplemental:
                                    for line in supplemental:
                                        if gene in line:
                                            allelelist.append(line)
                                    try:
                                        allelenumber = int(allelelist[-1].split("_")[1]) + 1
                                    except IndexError:
                                        pass
                                    # Puts the HSP in the correct order -  hits to the negative strand will be
                                    # reversed compared to what we're looking for
                                    allelesequence = Seq(hsp.query, generic_dna)
                                    if hsp.sbjct_start < hsp.sbjct_end:
                                        end = hsp.sbjct_end
                                    else:
                                        end = hsp.sbjct_start
                                        allelesequence = allelesequence.reverse_complement()
                                    # Screen out hits that are shorter than the targets
                                    # Keeping this format, though this if statement could be re-written more efficiently
                                    if end < alignment.length:
                                        pass
                                    definitionline = '{}_{}'.format(gene, allelenumber)
                                    # Create a sequence record using BioPython
                                    fasta = SeqRecord(allelesequence,
                                                      # If this is not added, the header will not be formatted properly
                                                      description='',
                                                      # Use >:definitionline as the header
                                                      id=definitionline)
                                    # Use the SeqIO module to properly format the new sequence record
                                    SeqIO.write(fasta, supplemental, "fasta")
                                    # Return the appropriate information
                                    return gene, allelenumber, 100.0, hsp.score
819512
def supported_alleles(self):
        
        if 'supported_alleles' not in self._cache:
            result = set(self.allele_to_allele_specific_models)
            if self.allele_to_fixed_length_sequence:
                result = result.union(self.allele_to_fixed_length_sequence)
            self._cache["supported_alleles"] = sorted(result)
        return self._cache["supported_alleles"]
645065
def mitosomatic(args):
    
    import pandas as pd

    p = OptionParser(mitosomatic.__doc__)
    p.add_option("--minaf", default=.005, type="float",
                 help="Minimum allele fraction")
    p.add_option("--maxaf", default=.1, type="float",
                 help="Maximum allele fraction")
    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(not p.print_help())

    df, = args
    af_file = df.rsplit(".", 1)[0] + ".af"
    fw = open(af_file, "w")
    df = pd.read_csv(df, sep="\t")
    for i, row in df.iterrows():
        na = row["num_A"]
        nt = row["num_T"]
        nc = row["num_C"]
        ng = row["num_G"]
        nd = row["num_D"]
        ni = row["num_I"]
        depth = row["depth"]
        #major, minor = sorted([na, nt, nc, ng], reverse=True)[:2]
        #af = minor * 1. / (major + minor)
        af = (nd + ni) * 1. / depth
        if not (opts.minaf <= af <= opts.maxaf):
            continue
        print("{}\t{}\t{:.6f}".format(row["chrom"], row["start"], af), file=fw)
    fw.close()

    logging.debug("Allele freq written to `{}`".format(af_file))
126603
def _organize_calls(out_file, hla_base, data):
    
    hla_truth = get_hla_truthset(data)
    sample = dd.get_sample_name(data)
    with file_transaction(data, out_file) as tx_out_file:
        with open(tx_out_file, "w") as out_handle:
            writer = csv.writer(out_handle)
            writer.writerow(["sample", "locus", "mismatches", "options", "alleles", "p-groups", "expected",
                             "validates"])
            for genotype_file in glob.glob("%s.HLA-*.gt" % (hla_base)):
                hla_locus = os.path.basename(genotype_file).replace(
                        "%s.HLA-" % os.path.basename(hla_base), "").replace(".gt", "")
                with open(genotype_file) as in_handle:
                    total_options = set([])
                    for i, line in enumerate(in_handle):
                        _, aone, atwo, m = line.split("\t")[:4]
                        pgroups = (hla_groups.hla_protein(aone, data), hla_groups.hla_protein(atwo, data))
                        if i == 0:
                            call_alleles = [aone, atwo]
                            call_pgroups = pgroups
                            mismatches = m
                        total_options.add(pgroups)
                    if len(total_options) > 0:
                        truth_alleles = tz.get_in([sample, hla_locus], hla_truth, [])
                        writer.writerow([sample, hla_locus, mismatches, len(total_options),
                                         ";".join(call_alleles), ";".join(call_pgroups),
                                         ";".join(truth_alleles), matches_truth(call_alleles, truth_alleles, data)])
    return out_file
1207778
def getrmlsthelper(self):
        
        # Set the path/name of the folder to contain the new alleles and profile
        newfolder = os.path.join(self.path, self.analysistype)
        # Create the path
        make_path(newfolder)
        # Create arguments to feed into the rest_auth_class script
        args = ArgumentParser
        args.secret_file = os.path.join(self.credentials, 'secret.txt')
        args.file_path = os.path.join(self.credentials)
        args.output_path = newfolder
        args.logging = self.logging
        rmlst = rest_auth_class.REST(args)
        rmlst.secret_finder()
        rmlst.get_request_token()
        rmlst.get_access_token()
        # Download the profile and alleles
        rmlst.main()
        # Get the new alleles into a list, and create the combinedAlleles file
        alleles = glob(os.path.join(newfolder, '*.tfa'))
        self.combinealleles(newfolder, alleles)
1249468
def mconsensus(args):
    
    p = OptionParser(mconsensus.__doc__)
    p.add_option("--allele_counts", default="allele_counts",
                 help="Directory to generate allele counts")
    add_consensus_options(p)
    opts, args = p.parse_args(args)

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

    consensusfiles = args
    prefix = opts.prefix
    acdir = opts.allele_counts
    store = ClustStores(consensusfiles)
    pctid = find_pctid(consensusfiles)
    pf = prefix + ".P{0}".format(pctid)

    clustSfile = pf + ".clustS"
    AC = makeloci(clustSfile, store, prefix,
                  minsamp=opts.minsamp, pctid=pctid)

    mkdir(acdir)
    acfile = pf + ".allele_counts"
    fw = open(acfile, "w")
    seen = DefaultOrderedDict(list)        # chr, pos => taxa
    print("# " + "\t".join(ACHEADER), file=fw)
    # Sort allele counts into separate files
    for ac in AC:
        chrpos = ac.chr, ac.pos
        seen[chrpos].append(ac)
        print(ac.tostring(taxon=True), file=fw)
    fw.close()

    logging.debug("Populate all taxa and instantiate empty vector if missing")
    all_taxa = set([op.basename(x).split(".")[0] for x in consensusfiles])
    taxon_to_ac = defaultdict(list)
    for chrpos, aclist in seen.items():
        included_taxa = set([x.taxon for x in aclist])
        missing_taxa = all_taxa - included_taxa
        template = deepcopy(aclist[0])
        template.clear()
        for ac in aclist:
            taxon_to_ac[ac.taxon].append(ac)
        for tx in missing_taxa:
            taxon_to_ac[tx].append(template)

    logging.debug("Write allele counts for all taxa")
    for tx, aclist in sorted(taxon_to_ac.items()):
        tx_acfile = op.join(acdir, tx + ".allele_counts")
        fw = open(tx_acfile, "w")
        print("# " + "\t".join(ACHEADER_NO_TAXON), file=fw)
        for ac in aclist:
            print(ac.tostring(), file=fw)
        fw.close()
        logging.debug("Written {0} sites in `{1}`".\
                format(len(aclist), tx_acfile))
686769
def do_download_datafiles(self, _):
        '''  '''
        contents = {"trialdata": lambda p: p.get_trial_data(), "eventdata": \
                    lambda p: p.get_event_data(), "questiondata": lambda p: \
                    p.get_question_data()}
        query = Participant.query.all()
        for k in contents:
            ret = "".join([contents[k](p) for p in query])
            temp_file = open(k + '.csv', 'w')
            temp_file.write(ret)
            temp_file.close()
743969
def getgenes(args):
    
    p = OptionParser(getgenes.__doc__)
    p.add_option("--prefix", default="gbout",
            help="prefix of output files [default: %default]")
    p.add_option("--nopep", default=False, action="store_true",
            help="Only get cds and bed, no pep [default: %default]")
    filenames, accessions, idfile, opts, args = preparegb(p, args)
    prefix = opts.prefix

    GenBank(filenames=filenames, accessions=accessions, idfile=idfile).\
        write_genes(output=prefix, individual=opts.individual, \
        pep=(not opts.nopep))

    if opts.individual:
        logging.debug("Output written dir {0}".format(prefix))
    elif opts.nopep:
        logging.debug("Output written to {0}.bed, {0}.cds".format(prefix,))
    else:
        logging.debug("Output written to {0}.bed, {0}.cds, {0}.pep".format(prefix,))
1044526
def load(self,cache_genotype=False,cache_phenotype=True):
        
        self.f = h5py.File(self.file_name,'r')
        self.pheno = self.f['phenotype']
        self.geno  = self.f['genotype']
        #TODO: load all row and column headers for genotype and phenotype

        #parse out thse we alwasy need for convenience
        self.genoM = self.geno['matrix']

        self.phenoM = self.pheno['matrix']
        self.sample_ID = self.geno['row_header']['sample_ID'][:]
        self.genoChrom = self.geno['col_header']['chrom'][:]
        self.genoPos   = self.geno['col_header']['pos'][:]
        if 'pos_cum' in list(self.geno['col_header'].keys()):
            self.genoPos_cum   = self.geno['col_header']['pos_cum'][:]
        else:
            self.genoPos_cum = None
        self.phenotype_ID = self.pheno['col_header']['phenotype_ID'][:]


        #cache?
        if cache_genotype:
            self.genoM = self.genoM[:]
        if cache_phenotype:
            self.phenoM = self.phenoM[:]
    
        # Additional pheno col header
        headers = list(self.pheno['col_header'].keys())
        if 'gene_ID' in headers:
            self.eqtl = True
            self.geneID = self.pheno['col_header']['gene_ID'][:]
            self.gene_pos = SP.array([self.pheno['col_header']['gene_chrom'][:],self.pheno['col_header']['gene_start'][:],self.pheno['col_header']['gene_end']],dtype='int').T
            self.geneIDs= list(set(self.geneID))
        else:
            self.eqtl = False
        if 'environment' in headers:
            self.E  = self.pheno['col_header/environment'][:]
            self.Es = list(set(self.E))
        else:
            self.E = None

        #dimensions
        self.N = self.genoM.shape[0]
        self.S = self.genoM.shape[1]
        self.P = self.phenoM.shape[1]
        assert (self.genoM.shape[0]==self.phenoM.shape[0]), 'dimension missmatch'
695618
def get_all_resources(self):
        
        resources = {}
        resources[
            "genetic_map_HapMapII_GRCh37"
        ] = self.get_genetic_map_HapMapII_GRCh37()
        resources["cytoBand_hg19"] = self.get_cytoBand_hg19()
        resources["knownGene_hg19"] = self.get_knownGene_hg19()
        resources["kgXref_hg19"] = self.get_kgXref_hg19()
        for source, target in itertools.permutations(["NCBI36", "GRCh37", "GRCh38"], 2):
            resources[source + "_" + target] = self.get_assembly_mapping_data(
                source, target
            )
        return resources
1121853
def download_kegg_gene_metadata(gene_id, outdir=None, force_rerun=False):
    
    if not outdir:
        outdir = ''

    # Replace colon with dash in the KEGG gene ID
    outfile = op.join(outdir, '{}.kegg'.format(custom_slugify(gene_id)))

    if ssbio.utils.force_rerun(flag=force_rerun, outfile=outfile):
        raw_text = bs_kegg.get("{}".format(gene_id))
        if raw_text == 404:
            return

        with io.open(outfile, mode='wt', encoding='utf-8') as f:
            f.write(raw_text)

        log.debug('{}: downloaded KEGG metadata file'.format(outfile))
    else:
        log.debug('{}: KEGG metadata file already exists'.format(outfile))

    return outfile
185465
function getFiles () {
  const FUCK_ENV = process.env.FUCK_ENV
    || process.env.npm_package_config_FUCK_ENV
    || FUCK_ENV_DEFAULT
  return FUCK_ENV.split(',').map(file => path.resolve(root, file))
}
669109
def group_reads_by_allele(allele_reads):
    
    allele_to_reads_dict = defaultdict(list)
    for allele_read in allele_reads:
        allele_to_reads_dict[allele_read.allele].append(allele_read)
    return allele_to_reads_dict
1076541
def allele_support_df(loci, sources):
    
    return pandas.DataFrame(
        allele_support_rows(loci, sources),
        columns=EXPECTED_COLUMNS)
1268326
def main():
    """"""
    from argparse import ArgumentParser
    description = "Download a bunch of thumbnails from Wikimedia Commons"
    parser = ArgumentParser(description=description)
    source_group = parser.add_mutually_exclusive_group()
    source_group.add_argument("-l", "--list", metavar="LIST",
                              dest="file_list",
                              type=argparse.FileType('r'),
                              help='A list of files <filename,width>')
    source_group.add_argument("-c", "--category", metavar="CATEGORY",
                              dest="category_name",
                              type=str,
                              help='A category name (without prefix)')
    parser.add_argument("files", nargs='*',
                        metavar="FILES",
                        help='A list of filenames')
    parser.add_argument("-o", "--output", metavar="FOLDER",
                        dest="output_path",
                        action=Folder,
                        default=os.getcwd(),
                        help='The directory to download the files to')
    parser.add_argument("-w", "--width",
                        dest="width",
                        type=int,
                        default=100,
                        help='The width of the thumbnail (default: 100)')
    verbosity_group = parser.add_mutually_exclusive_group()
    verbosity_group.add_argument("-v",
                                 action="count",
                                 dest="verbose",
                                 default=1,
                                 help="Verbosity level. -v for DEBUG")
    verbosity_group.add_argument("-q", "--quiet",
                                 action="store_const",
                                 dest="verbose",
                                 const=0,
                                 help="To silence the INFO messages")
    args = parser.parse_args()
    logging_map = {0: logging.WARNING,
                   1: logging.INFO,
                   2: logging.DEBUG}
    logging_level = logging_map.get(args.verbose, logging.DEBUG)
    logging.basicConfig(level=logging_level)
    logging.info("Starting")

    if args.file_list:
        download_from_file_list(args.file_list, args.output_path)
    elif args.category_name:
        download_from_category(args.category_name, args.output_path, args.width)
    elif args.files:
        download_from_files(args.files, args.output_path, args.width)
    else:
        parser.print_help()
1063316
def process_mgi_note_allele_view(self, limit=None):
        

        line_counter = 0
        if self.test_mode:
            graph = self.testgraph
        else:
            graph = self.graph
        model = Model(graph)
        LOG.info("Assembling notes on alleles")
        raw = '/'.join((self.rawdir, 'mgi_note_allele_view'))

        notehash = {}
        with open(raw, 'r', encoding="utf8") as csvfile:
            filereader = csv.reader(csvfile, delimiter='\t', quotechar='\"')
            for line in filereader:
                line_counter += 1
                if line_counter == 1:
                    continue

                (object_key, notetype, note, sequencenum) = line

                # read all the notes into a hash to concatenate
                if object_key not in notehash:
                    notehash[object_key] = {}
                if notetype not in notehash[object_key]:
                    notehash[object_key][notetype] = []
                if len(notehash[object_key][notetype]) < int(sequencenum):
                    for i in range(
                            len(notehash[object_key][notetype]),
                            int(sequencenum)
                    ):
                        notehash[object_key][notetype].append('')  # ??? I don't get it

                notehash[object_key][notetype][int(sequencenum)-1] = note.strip()

            # finish iteration over notes

        line_counter = 0
        for allele_key in notehash:
            if self.test_mode is True:
                if int(allele_key) not in self.test_keys.get('allele'):
                    continue
            line_counter += 1
            allele_id = self.idhash['allele'].get(allele_key)
            if allele_id is None:
                continue
            for n in notehash[allele_key]:
                LOG.info(
                    "found %d %s notes for %s",
                    len(notehash[allele_key]), n, allele_id)
                notes = ''.join(notehash[allele_key][n])
                notes += ' ['+n+']'
                model.addDescription(allele_id, notes)

            if not self.test_mode and limit is not None and line_counter > limit:
                break

        return