Acquire a genome from Entrez
190793
def fetch_genome(genome_id):
    '''

    '''
    # TODO: Can strandedness by found in fetched genome attributes?
    # TODO: skip read/write step?
    # Using a dummy email for now - does this violate NCBI guidelines?
    email = 'loremipsum@gmail.com'
    Entrez.email = email

    print 'Downloading Genome...'
    handle = Entrez.efetch(db='nucleotide', id=str(genome_id), rettype='gb',
                           retmode='text')
    print 'Genome Downloaded...'
    tmpfile = os.path.join(mkdtemp(), 'tmp.gb')
    with open(tmpfile, 'w') as f:
        f.write(handle.read())
    genome = coral.seqio.read_dna(tmpfile)

    return genome
518088
def importGenome(name, batchSize = 100) :
	""""""
	path = os.path.join(this_dir, "bootstrap_data", "genomes/" + name)
	PG.importGenome(path, batchSize)
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]
788414
def _incorporate_dup(self):
        """"""
        seq, cds_start, cds_stop, start, end = self._setup_incorporate()

        dup_seq = seq[start:end]
        seq[end:end] = dup_seq

        is_frameshift = len(dup_seq) % 3 != 0
        variant_start_aa = int(math.ceil((self._var_c.posedit.pos.end.base + 1) / 3.0))

        alt_data = AltTranscriptData(
            seq,
            cds_start,
            cds_stop,
            is_frameshift,
            variant_start_aa,
            self._transcript_data.protein_accession,
            is_ambiguous=self._ref_has_multiple_stops)
        return alt_data
1223522
def _upgrade_genome_resources(galaxy_dir, base_url):
    
    import requests
    for dbkey, ref_file in genome.get_builds(galaxy_dir):
        # Check for a remote genome resources file
        remote_url = base_url % dbkey
        requests.packages.urllib3.disable_warnings()
        r = requests.get(remote_url, verify=False)
        if r.status_code == requests.codes.ok:
            local_file = os.path.join(os.path.dirname(ref_file), os.path.basename(remote_url))
            if os.path.exists(local_file):
                with open(local_file) as in_handle:
                    local_config = yaml.safe_load(in_handle)
                remote_config = yaml.safe_load(r.text)
                needs_update = remote_config["version"] > local_config.get("version", 0)
                if needs_update:
                    shutil.move(local_file, local_file + ".old%s" % local_config.get("version", 0))
            else:
                needs_update = True
            if needs_update:
                print("Updating %s genome resources configuration" % dbkey)
                with open(local_file, "w") as out_handle:
                    out_handle.write(r.text)
341093
def from_command_line():
    
    # Parse options
    parser = argparse.ArgumentParser(
        description='Convert Complete Genomics var files to gVCF format.')
    parser.add_argument(
        '-d', '--refseqdir', metavar='REFSEQDIR', required=True,
        dest='refseqdir',
        help='Directory twobit reference genomes files are stored.')
    parser.add_argument(
        '-i', '--input', metavar='INPUTVARFILE',
        dest='cgivarfile',
        help='Path to Complete Genomics var file to convert. If omitted, data '
        ' also be piped in as standard input.')
    parser.add_argument(
        '-o', '--output', metavar='OUTPUTVCFFILE',
        dest='vcfoutfile',
        help='Path to where to save output VCF file.')
    parser.add_argument(
        '-D', '--download', action='store_true', dest='downloadrefseq',
        help='Download the 2bit file from UCSC to REFSEQDIR, if needed.')
    parser.add_argument(
        '-v', '--var-only', action='store_true', dest='varonly',
        help='Only report variant lines (i.e. VCF, but not gVCF)')
    args = parser.parse_args()

    # Get local twobit file from its directory. Download and store if needed.
    twobit_path, twobit_name = get_reference_genome_file(
        args.refseqdir, build='b37')
    # Handle input
    if sys.stdin.isatty():  # false if data is piped in
        var_input = args.cgivarfile
    else:
        var_input = sys.stdin
    # Handle output
    if args.vcfoutfile:
        convert_to_file(var_input,
                        args.vcfoutfile,
                        twobit_path,
                        twobit_name,
                        args.varonly)
    else:
        for line in convert(
                cgi_input=var_input,
                twobit_ref=twobit_path,
                twobit_name=twobit_name,
                var_only=args.varonly):
            print(line)
258465
def _run_purecn(paired, work_dir):
    
    segfns = {"cnvkit": _segment_normalized_cnvkit, "gatk-cnv": _segment_normalized_gatk}
    out_base, out, all_files = _get_purecn_files(paired, work_dir)
    failed_file = out_base + "-failed.log"
    cnr_file = tz.get_in(["depth", "bins", "normalized"], paired.tumor_data)
    if not utils.file_uptodate(out["rds"], cnr_file) and not utils.file_exists(failed_file):
        cnr_file, seg_file = segfns[cnvkit.bin_approach(paired.tumor_data)](cnr_file, work_dir, paired)
        from bcbio import heterogeneity
        vcf_file = heterogeneity.get_variants(paired.tumor_data, include_germline=False)[0]["vrn_file"]
        vcf_file = germline.filter_to_pass_and_reject(vcf_file, paired, out_dir=work_dir)
        with file_transaction(paired.tumor_data, out_base) as tx_out_base:
            # Use UCSC style naming for human builds to support BSgenome
            genome = ("hg19" if dd.get_genome_build(paired.tumor_data) in ["GRCh37", "hg19"]
                      else dd.get_genome_build(paired.tumor_data))
            cmd = ["PureCN.R", "--seed", "42", "--out", tx_out_base, "--rds", "%s.rds" % tx_out_base,
                   "--sampleid", dd.get_sample_name(paired.tumor_data),
                   "--genome", genome,
                   "--vcf", vcf_file, "--tumor", cnr_file,
                   "--segfile", seg_file, "--funsegmentation", "Hclust", "--maxnonclonal", "0.3"]
            if dd.get_num_cores(paired.tumor_data) > 1:
                cmd += ["--cores", str(dd.get_num_cores(paired.tumor_data))]
            try:
                cmd = "export R_LIBS_USER=%s && %s && %s" % (utils.R_sitelib(), utils.get_R_exports(),
                                                             " ".join([str(x) for x in cmd]))
                do.run(cmd, "PureCN copy number calling")
            except subprocess.CalledProcessError as msg:
                if _allowed_errors(str(msg)):
                    logger.info("PureCN failed to find solution for %s: skipping" %
                                dd.get_sample_name(paired.tumor_data))
                    with open(failed_file, "w") as out_handle:
                        out_handle.write(str(msg))
                else:
                    logger.exception()
                    raise
            for f in all_files:
                if os.path.exists(os.path.join(os.path.dirname(tx_out_base), f)):
                    shutil.move(os.path.join(os.path.dirname(tx_out_base), f),
                                os.path.join(os.path.dirname(out_base), f))
    out = _get_purecn_files(paired, work_dir, require_exist=True)[1]
    return out if (out.get("rds") and os.path.exists(out["rds"])) else None
254664
def get_genomic_seq_for_transcript(self, transcript_id, expand):
        
        
        headers = {"content-type": "application/json"}
        
        self.attempt = 0
        ext = "/sequence/id/{0}?type=genomic;expand_3prime={1};expand_5prime={1}".format(transcript_id, expand)
        r = self.ensembl_request(ext, headers)
        
        gene = json.loads(r)
        
        seq = gene["seq"]
        seq_id = gene["id"]
        
        if seq_id != transcript_id:
            raise ValueError("ensembl gave the wrong transcript")
        
        desc = gene["desc"].split(":")
        chrom = desc[2]
        start = int(desc[3]) + expand
        end = int(desc[4]) - expand
        strand_temp = int(desc[5])
        
        strand = "+"
        if strand_temp == -1:
            strand = "-"
        
        return (chrom, start, end, strand, seq)
218593
def get_column(self, chrom, position,
                 missing_seqs=MissingSequenceHandler.TREAT_AS_ALL_GAPS,
                 species=None):
    """"""
    blocks = self.get_blocks(chrom, position, position + 1)
    if len(blocks) == 0:
      raise NoSuchAlignmentColumnError("Request for column on chrom " +
                                       chrom + " at position " +
                                       str(position) + " not possible; " +
                                       "genome alignment not defined at " +
                                       "that locus.")
    if len(blocks) > 1:
      raise NoUniqueColumnError("Request for column on chrom " + chrom +
                                " at position " + str(position) + "not " +
                                "possible; ambiguous alignment of that locus.")

    return blocks[0].get_column_absolute(position, missing_seqs, species)
881597
def _create_genome_builds(self):
        

        # TODO add more species

        graph = self.graph
        geno = Genotype(graph)
        model = Model(graph)
        LOG.info("Adding equivalent assembly identifiers")
        for sp in self.species:
            tax_id = self.globaltt[sp]
            txid_num = tax_id.split(':')[1]
            for key in self.files[txid_num]['assembly']:
                ucsc_id = key
                try:
                    ucsc_label = ucsc_id.split(':')[1]
                except IndexError:
                    LOG.error('%s Assembly id:  "%s" is problematic', sp, key)
                    continue
                if key in self.localtt:
                    mapped_id = self.localtt[key]
                else:
                    LOG.error(
                        '%s Assembly id:  "%s" is not in local translation table',
                        sp, key)

                mapped_label = mapped_id.split(':')[1]

                mapped_label = 'NCBI build ' + str(mapped_label)
                geno.addReferenceGenome(ucsc_id, ucsc_label, tax_id)
                geno.addReferenceGenome(mapped_id, mapped_label, tax_id)
                model.addSameIndividual(ucsc_id, mapped_id)

        return
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
1330502
def parse(self, limit=None):
        
        if limit is not None:
            LOG.info("Only parsing first %d rows", limit)

        protein_paths = self._get_file_paths(self.tax_ids, 'protein_links')
        col = ['NCBI taxid', 'entrez', 'STRING']
        for taxon in protein_paths:
            ensembl = Ensembl(self.graph_type, self.are_bnodes_skized)
            string_file_path = '/'.join((
                self.rawdir, protein_paths[taxon]['file']))

            with gzip.open(string_file_path, 'rb') as reader:
                dataframe = pd.read_csv(reader, sep=r'\s+')
            p2gene_map = dict()

            if taxon in self.id_map_files:
                LOG.info("Using string provided id_map files")
                map_file = '/'.join((self.rawdir, self.id_map_files[taxon]['file']))

                with gzip.open(map_file, 'rt') as reader:
                    line = next(reader).strip()
                    if line != '# NCBI taxid / entrez / STRING':
                        LOG.error(
                            'Expected Headers:\t%s\nRecived Headers:\t%s\n', col, line)
                        exit(-1)

                    for line in reader.readlines():
                        row = line.rstrip('\n').split('\t')
                        # tax = row[col.index(''NCBI taxid')].strip()
                        gene = row[col.index('entrez')].strip()
                        prot = row[col.index('STRING')].strip()

                        genes = gene.split('|')
                        p2gene_map[prot.replace(taxon + '.', '')] = [
                            "NCBIGene:" + entrez_id for entrez_id in genes]
            else:
                LOG.info("Fetching ensembl proteins for taxon %s", taxon)
                p2gene_map = ensembl.fetch_protein_gene_map(taxon)
                for key in p2gene_map:
                    for phen, gene in enumerate(p2gene_map[key]):
                        p2gene_map[key][phen] = "ENSEMBL:{}".format(gene)

            LOG.info(
                "Finished fetching ENSP ID mappings, fetched %i proteins",
                len(p2gene_map))

            LOG.info(
                "Fetching protein protein interactions for taxon %s", taxon)

            self._process_protein_links(dataframe, p2gene_map, taxon, limit)
724720
def add_reference_resources(data, remote_retriever=None):
    
    aligner = data["config"]["algorithm"].get("aligner", None)
    if remote_retriever:
        data["reference"] = remote_retriever.get_refs(data["genome_build"],
                                                      alignment.get_aligner_with_aliases(aligner, data),
                                                      data["config"])
    else:
        data["reference"] = genome.get_refs(data["genome_build"], alignment.get_aligner_with_aliases(aligner, data),
                                            data["dirs"]["galaxy"], data)
        _check_ref_files(data["reference"], data)
    # back compatible `sam_ref` target
    data["sam_ref"] = utils.get_in(data, ("reference", "fasta", "base"))
    ref_loc = utils.get_in(data, ("config", "resources", "species", "dir"),
                           utils.get_in(data, ("reference", "fasta", "base")))
    if remote_retriever:
        data = remote_retriever.get_resources(data["genome_build"], ref_loc, data)
    else:
        data["genome_resources"] = genome.get_resources(data["genome_build"], ref_loc, data)
    data["genome_resources"] = genome.add_required_resources(data["genome_resources"])
    if effects.get_type(data) == "snpeff" and "snpeff" not in data["reference"]:
        data["reference"]["snpeff"] = effects.get_snpeff_files(data)
    if "genome_context" not in data["reference"]:
        data["reference"]["genome_context"] = annotation.get_context_files(data)
    if "viral" not in data["reference"]:
        data["reference"]["viral"] = viral.get_files(data)
    if not data["reference"]["viral"]:
        data["reference"]["viral"] = None
    if "versions" not in data["reference"]:
        data["reference"]["versions"] = _get_data_versions(data)

    data = _fill_validation_targets(data)
    data = _fill_prioritization_targets(data)
    data = _fill_capture_regions(data)
    # Re-enable when we have ability to re-define gemini configuration directory
    if False:
        data["reference"]["gemini"] = population.get_gemini_files(data)
    return data
1271376
def convert_uniprot_to_entrez(self, uniprot):
        """"""
        # Submit request to NCBI eutils/Gene Database
        server = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?" + self.options + "&db=gene&term={0}".format(
            uniprot)
        r = requests.get(server, headers={"Content-Type": "text/xml"})
        if not r.ok:
            r.raise_for_status()
            sys.exit()
        # Process Request
        response = r.text
        info = xmltodict.parse(response)
        geneId = info['eSearchResult']['IdList']['Id']
        # check to see if more than one result is returned
        # if you have more than more result then check which Entrez Id returns the same uniprot Id entered.
        if len(geneId) > 1:
            for x in geneId:
                c = self.convert_entrez_to_uniprot(x)
                c = c.lower()
                u = uniprot.lower()
                if c == u:
                    return x
        else:
            return geneId
1151992
def _prep_genome(out_dir, data):
    
    genome_name = utils.splitext_plus(os.path.basename(dd.get_ref_file(data)))[0]
    out_dir = utils.safe_makedir(os.path.join(out_dir, genome_name))
    ref_file = dd.get_ref_file(data)
    utils.symlink_plus(ref_file, os.path.join(out_dir, os.path.basename(ref_file)))
    with open(os.path.join(out_dir, "GenomeSize.xml"), "w") as out_handle:
        out_handle.write('<sequenceSizes genomeName="%s">' % genome_name)
        for c in pysam.AlignmentFile("%s.dict" % utils.splitext_plus(ref_file)[0]).header["SQ"]:
            cur_ploidy = ploidy.get_ploidy([data], region=[c["SN"]])
            out_handle.write('<chromosome fileName="%s" contigName="%s" totalBases="%s" knownBases="%s" '
                             'isCircular="false" ploidy="%s" md5="%s"/>' %
                             (os.path.basename(ref_file), c["SN"], c["LN"], c["LN"], cur_ploidy, c["M5"]))
        out_handle.write('</sequenceSizes>')
    return out_dir
356759
def get_genomic_seq_for_region(self, chrom, start_pos, end_pos):
        
        
        headers = {"content-type": "text/plain"}
        
        self.attempt = 0
        ext = "/sequence/region/human/{}:{}..{}:1".format(chrom, start_pos, end_pos)
        
        return self.ensembl_request(ext, headers)
1232060
def createRepo(self):
        
        repo = datarepo.SqlDataRepository(self.repoPath)
        repo.open("w")
        repo.initialise()

        referenceSet = references.HtslibReferenceSet("GRCh37-subset")
        referenceSet.populateFromFile(self.fastaFilePath)
        referenceSet.setDescription("Subset of GRCh37 used for demonstration")
        referenceSet.setSpeciesFromJson(
                '{"id": "9606",'
                + '"term": "Homo sapiens", "source_name": "NCBI"}')
        for reference in referenceSet.getReferences():
            reference.setSpeciesFromJson(
                '{"id": "9606",'
                + '"term": "Homo sapiens", "source_name": "NCBI"}')
            reference.setSourceAccessions(
                self.accessions[reference.getName()] + ".subset")
        repo.insertReferenceSet(referenceSet)

        dataset = datasets.Dataset("1kg-p3-subset")
        dataset.setDescription("Sample data from 1000 Genomes phase 3")
        repo.insertDataset(dataset)

        variantSet = variants.HtslibVariantSet(dataset, "mvncall")
        variantSet.setReferenceSet(referenceSet)
        dataUrls = [vcfFile for vcfFile, _ in self.vcfFilePaths]
        indexFiles = [indexFile for _, indexFile in self.vcfFilePaths]
        variantSet.populateFromFile(dataUrls, indexFiles)
        variantSet.checkConsistency()
        repo.insertVariantSet(variantSet)

        for sample, (bamFile, indexFile) in zip(
                self.samples, self.bamFilePaths):
            readGroupSet = reads.HtslibReadGroupSet(dataset, sample)
            readGroupSet.populateFromFile(bamFile, indexFile)
            readGroupSet.setReferenceSet(referenceSet)
            repo.insertReadGroupSet(readGroupSet)

        repo.commit()
        repo.close()
        self.log("Finished creating the repository; summary:\n")
        repo.open("r")
        repo.printSummary()