Parse ensembl gene to reactome pathway file
        :param file: file path (not handle)
        :param limit: limit (int, optional) limit the number of rows processed
        :return: None
502417
def _parse_reactome_association_file(
            self, file, limit=None, subject_prefix=None, object_prefix=None):
        
        eco_map = Reactome.get_eco_map(Reactome.map_files['eco_map'])
        count = 0
        with open(file, 'r') as tsvfile:
            reader = csv.reader(tsvfile, delimiter="\t")
            for row in reader:
                (component, pathway_id, pathway_iri, pathway_label, go_ecode,
                species_name) = row
                count += 1
                self._add_component_pathway_association(
                    eco_map, component, subject_prefix, pathway_id, object_prefix,
                    pathway_label, go_ecode)

                if limit is not None and count >= limit:
                    break

        return
43344
def get_resources(genome, ref_file, data):
    
    base_dir = os.path.normpath(os.path.dirname(ref_file))
    resource_file = os.path.join(base_dir, "%s-resources.yaml" % genome.replace("-test", ""))
    if not os.path.exists(resource_file):
        raise IOError("Did not find resource file for %s: %s\n"
                      "To update bcbio_nextgen.py with genome resources for standard builds, run:\n"
                      "bcbio_nextgen.py upgrade -u skip"
                      % (genome, resource_file))
    with open(resource_file) as in_handle:
        resources = yaml.safe_load(in_handle)

    def resource_file_path(x):
        if isinstance(x, six.string_types) and os.path.exists(os.path.join(base_dir, x)):
            return os.path.normpath(os.path.join(base_dir, x))
        return x
    cleaned = utils.dictapply(resources, resource_file_path)
    return ensure_annotations(cleaned, data)
330209
def _parse_gene_coordinates(file):
        
        id_map = {}
        col = ['gene_curie', 'start', 'end', 'strand', 'build']
        if os.path.exists(os.path.join(os.path.dirname(__file__), file)):
            with open(os.path.join(os.path.dirname(__file__), file)) as tsvfile:
                reader = csv.reader(tsvfile, delimiter="\t")
                for row in reader:
                    id_map[row[col.index('gene_curie')]] = {
                        'start': row[col.index('start')],
                        'end': row[col.index('end')],
                        'strand': row[col.index('strand')],
                        'build': row[col.index('build')]
                    }
        return id_map
100763
def libseq_parse(self, line):
        """"""
        count = line.count(',')
        if count > 1:
            raise ValueError("""ERROR: There are too many commas in your
            lib_seq line. There should only be one, between the two filepath
            globs.""")
        if count == 0:
            raise ValueError("""ERROR: There are no commas in your
            lib_seq line. There should be one, between the two filepath
            globs.""")
        line = [x.strip() for x in line.split() if x]
        #print(line)

        #This block tries to fix user input if there is an inappropriate
        # number of commas or spaces associated with commas
        if len(line) > 13:
            commas = [x for x in line if "," in x]
            if commas:
                #this handles the case if there was an extra space added on
                # both sides of the comma
                if (len(commas) == 1) and (commas[0] != ','):
                    if (commas[0][-1] == ',') or (commas[0][0] == ','):
                        line = [line[0]] + ["".join(line[1:3])] + line[3:]
                #This handles the case where there is an extra space on either
                #side of the comma
                elif line[2] == ',':
                    line = [line[0]] + ["".join(line[1:4])] + line[4::]

        #make sure there is an appropriate numer of elements in list
        if (len(line) ) != 13:
            raise ValueError("""ERROR: There are {0} values in the lib_seq line when
                             there should be 12. Check your input meraculous
                             config file for errors and make sure that it matches
                             the specification for the Meraculous manual.""".format(
                                 len(line)))

        for index in self.indices:
            if index == 1:
                globs = [os.path.abspath(x.strip()) for x in line[1].split(',') if x]
                #make sure there are only two filepaths to search for globs
                if len(globs) != 2:
                    raise ValueError("""ERROR: remember to split your glob
                    filepaths in lib_seq with a comma. This error occured because
                    the program was looking for two glob filepaths, but instead
                    found {}: {}""".format(len(globs), globs))
                self.data["globs"] = globs
                #Note that glob.glob returns an empty list if the file or glob
                # filepath does not exist.
                forwards = glob.glob(os.path.abspath(globs[0]))
                reverses = glob.glob(os.path.abspath(globs[1]))
                # print(os.path.abspath(globs[0]))
                # print("globs: ", globs)
                # print("forwards: ", forwards)
                # print("reverses: ", reverses)
                #At this point the
                if (len(forwards) < 1) or (len(reverses) < 1):
                    # print("forwards: ", forwards)
                    # print("reverses: ", reverses)

                    report = []
                    if len(forwards) < 1:
                        report.append(globs[0])
                    if len(reverses) < 1:
                        report.append(globs[1])
                    raise ValueError("""ERROR: no read files were found for the glob strings,
                    {}.

                    Please check the glob string and/or use absolute
                    file paths. A common source of this error is using
                    relative filepaths from the incorrect directory, or
                    referring to files that do not exist.""".format(report))
                if len(forwards) != len(reverses):
                    raise ValueError("""ERROR: the number of reverse read files does not
                    match the number of forward read files""")
                #zip the pairs together into tuples, not lists
                self.data["pairs"] = [x for x in zip(sorted(forwards), sorted(reverses))]

            self.data[self.indices.get(index)] = line[index]
300029
def _process_pathway_pathway(self, limit):
        
        LOG.info("Processing KEGG pathways to other ids")
        if self.test_mode:
            graph = self.testgraph
        else:
            graph = self.graph
        line_counter = 0

        model = Model(graph)
        raw = '/'.join((self.rawdir, self.files['pathway_pathway']['file']))
        with open(raw, 'r', encoding="iso-8859-1") as csvfile:
            filereader = csv.reader(csvfile, delimiter='\t', quotechar='\"')
            for row in filereader:
                line_counter += 1
                (pathway_id_1, pathway_id_2) = row

                if self.test_mode and pathway_id_1 not in self.test_ids['pathway']:
                    continue

                pathway_id_1 = 'KEGG-' + pathway_id_1
                # will look like KEGG-path:map04130 or KEGG-path:ko04130
                pathway_id_2 = 'KEGG-' + pathway_id_2

                if pathway_id_1 != pathway_id_2:
                    model.addEquivalentClass(pathway_id_1, pathway_id_2)

                if not self.test_mode and limit is not None and line_counter > limit:
                    break
        return
1031391
def _read(
    filename,
    schema,
    seq_label='sequence',
    alphabet=None,
    use_uids=True,
    **kwargs):
    
    # Check Alphabet if given
    if alphabet is None:
        alphabet = Bio.Alphabet.Alphabet()

    elif alphabet in ['dna', 'rna', 'protein', 'nucleotide']:
        alphabet = getattr(Bio.Alphabet, 'generic_{}'.format(alphabet))

    else:
        raise Exception(
            "The alphabet is not recognized. Must be 'dna', 'rna', "
            "'nucleotide', or 'protein'.")

    kwargs.update(alphabet=alphabet)

    # Prepare DataFrame fields.
    data = {
        'id': [],
        seq_label: [],
        'description': [],
        'label': []
    }
    if use_uids:
        data['uid'] = []

    # Parse Fasta file.
    for i, s in enumerate(SeqIO.parse(filename, format=schema, **kwargs)):
        data['id'].append(s.id)
        data[seq_label].append(str(s.seq))
        data['description'].append(s.description)
        data['label'].append(s.name)

        if use_uids:
            data['uid'].append(get_random_id(10))

    # Port to DataFrame.
    return pd.DataFrame(data)
1064170
def load_bioassembly_info_from_file(self, biomol_num):
        

        # current functionality is to take in a pre-assembled bioassembly file, parse it's MODELs and get info from that
        pass
        from Bio.PDB import PDBParser
        p = PDBParser(PERMISSIVE=True, QUIET=True)

        structure = StructureIO(self.structure_path, self.file_type)

        structure = p.get_structure('tmp', infile)
        bioass_to_chain_stoich = defaultdict(int)
        for model in structure:
            for chain in model:
                bioass_to_chain_stoich[chain.id] += 1
        return dict(bioass_to_chain_stoich)
960489
def prep_seq2c_bed(data):
    
    if dd.get_background_cnv_reference(data, "seq2c"):
        bed_file = _background_to_bed(dd.get_background_cnv_reference(data, "seq2c"), data)
    else:
        bed_file = regions.get_sv_bed(data)
    if bed_file:
        bed_file = bedutils.clean_file(bed_file, data, prefix="svregions-")
    else:
        bed_file = bedutils.clean_file(dd.get_variant_regions(data), data)
    if not bed_file:
        return None

    col_num = bt.BedTool(bed_file).field_count()
    if col_num < 4:
        annotated_file = annotate.add_genes(bed_file, data, max_distance=0)
        if annotated_file == bed_file:
            raise ValueError("BED file for Seq2C must be annotated with gene names, "
                             "however the input BED is 3-columns and we have no transcript "
                             "data to annotate with " + bed_file)
        annotated_file = annotate.gene_one_per_line(annotated_file, data)
    else:
        annotated_file = bed_file

    ready_file = "%s-seq2cclean.bed" % (utils.splitext_plus(annotated_file)[0])
    if not utils.file_uptodate(ready_file, annotated_file):
        bed = bt.BedTool(annotated_file)
        if col_num > 4 and col_num != 8:
            bed = bed.cut(range(4))
        bed = bed.filter(lambda x: x.name not in ["", ".", "-"])
        with file_transaction(data, ready_file) as tx_out_file:
            bed.saveas(tx_out_file)
        logger.debug("Saved Seq2C clean annotated ready input BED into " + ready_file)

    return ready_file
1170579
def process_catalog(self, limit=None):
        
        raw = '/'.join((self.rawdir, self.files['catalog']['file']))
        LOG.info("Processing Data from %s", raw)
        efo_ontology = RDFGraph(False, "EFO")
        LOG.info("Loading EFO ontology in separate rdf graph")
        efo_ontology.parse(self.files['efo']['url'], format='xml')
        efo_ontology.bind_all_namespaces()
        LOG.info("Finished loading EFO ontology")

        so_ontology = RDFGraph(False, "SO")
        LOG.info("Loading SO ontology in separate rdf graph")
        so_ontology.parse(self.files['so']['url'], format='xml')
        so_ontology.bind_all_namespaces()
        LOG.info("Finished loading SO ontology")

        with open(raw, 'r', encoding="iso-8859-1") as csvfile:
            filereader = csv.reader(csvfile, delimiter='\t')
            header = next(filereader, None)  # the header row
            header_len = len(header)
            LOG.info('header length:\t %i', header_len)

            for row in filereader:
                if not row:
                    pass
                else:
                    if header_len != len(row):
                        LOG.error('BadRow: %i has %i columns', filereader.line_num, row)

                    (date_added_to_catalog,
                     pubmed_num,
                     first_author,
                     pub_date,
                     journal,
                     link,
                     study_name,
                     disease_or_trait,
                     initial_sample_description,
                     replicate_sample_description,
                     region,
                     chrom_num,
                     chrom_pos,
                     reported_gene_nums,
                     mapped_gene,
                     upstream_gene_num,
                     downstream_gene_num,
                     snp_gene_nums,
                     upstream_gene_distance,
                     downstream_gene_distance,
                     strongest_snp_risk_allele,
                     snps,
                     merged,
                     snp_id_current,
                     context,
                     intergenic_flag,
                     risk_allele_frequency,
                     pvalue,
                     pvalue_mlog,
                     pvalue_text,
                     or_or_beta,
                     confidence_interval_95,
                     platform_with_snps_passing_qc,
                     cnv_flag,
                     mapped_trait,
                     mapped_trait_uri,
                     study_accession,
                     GENOTYPING_TECHNOLOGY
                    ) = row

                    if self.test_mode:
                        continue

# 06-May-2015	25917933
#   Zai CC	20-Nov-2014	J Psychiatr Res	http://europepmc.org/abstract/MED/25917933
# A genome-wide association study of suicide severity scores in bipolar disorder.
# Suicide in bipolar disorder
# 959 European ancestry individuals	NA
# 10p11.22	10	32704340	C10orf68, CCDC7, ITGB1	CCDC7
# rs7079041-A	rs7079041	0	7079041	intron	0		2E-6	5.698970

                    variant_curie, variant_type = self._get_curie_and_type_from_id(
                        strongest_snp_risk_allele)

                    if strongest_snp_risk_allele.strip() == '':
                        LOG.debug(
                            "No strongest SNP risk allele for %s:\n%s",
                            pubmed_num, str(row))
                        # still consider adding in the EFO terms
                        # for what the study measured?
                        continue

                    if variant_type == 'snp':
                        self._add_snp_to_graph(
                            variant_curie, strongest_snp_risk_allele, chrom_num,
                            chrom_pos, context, risk_allele_frequency)

                        self._add_deprecated_snp(
                            variant_curie, snp_id_current, merged, chrom_num, chrom_pos)

                        self._add_snp_gene_relation(
                            variant_curie, snp_gene_nums, upstream_gene_num,
                            downstream_gene_num)
                    elif variant_type == 'haplotype':
                        self._process_haplotype(
                            variant_curie, strongest_snp_risk_allele, chrom_num,
                            chrom_pos, context, risk_allele_frequency, mapped_gene,
                            so_ontology)
                    elif variant_type is None:
                        LOG.warning(
                            "There's a snp id i can't manage: %s",
                            strongest_snp_risk_allele)
                        continue

                    description = self._make_description(
                        disease_or_trait, initial_sample_description,
                        replicate_sample_description,
                        platform_with_snps_passing_qc, pvalue)

                    self._add_variant_trait_association(
                        variant_curie, mapped_trait_uri, efo_ontology,
                        pubmed_num, description)

                    if not self.test_mode and (
                            limit is not None and filereader.line_num > limit):
                        break

        # TODO loop through the location hash,
        # and make all snps at that location equivalent
        for l in self.id_location_map:
            snp_ids = self.id_location_map[l]
            if len(snp_ids) > 1:
                LOG.info("%s has >1 snp id: %s", l, str(snp_ids))
        return
210823
private static String parseFilePath(String file)
    {
        final String excludePattern = "/META-INF/persistence.xml";
        file = file.substring(0, file.length() - excludePattern.length());

        // in case protocol is "file".
        file = file.endsWith("!") ? file.substring(0, file.length() - 1) : file;

        return file;
    }
194389
def parse_catalytic(insertion, gff):
    
    offset = insertion['offset']
    GeneStrand = insertion['strand']
    if type(insertion['intron']) is not str:
        return gff
    for intron in parse_fasta(insertion['intron'].split('|')):
        ID, annot, strand, pos = intron[0].split('>')[1].split()
        Start, End = [int(i) for i in pos.split('-')]
        if strand != GeneStrand:
            if strand == '+':
                strand = '-'
            else:
                strand = '+'
            Start, End = End - 2, Start - 2
        Start, End = abs(Start + offset) - 1, abs(End + offset) - 1
        gff['#seqname'].append(insertion['ID'])
        gff['source'].append('Rfam')
        gff['feature'].append('Catalytic RNA')
        gff['start'].append(Start)
        gff['end'].append(End)
        gff['score'].append('.')
        gff['strand'].append(strand)
        gff['frame'].append('.')
        gff['attribute'].append('ID=%s; Name=%s' % (ID, annot))
    return gff
85137
def load_just_in_time_genome_alignment(path, ref_spec, extensions=None,
                                       index_exts=None, fail_no_index=True,
                                       verbose=False):
    """"""
    if index_exts is None and fail_no_index:
      raise ValueError("Failure on no index specified for loading genome " +
                       "alignment, but no index extensions specified")

    extensions = __trim_extensions_dot(extensions)
    index_exts = __trim_extensions_dot(index_exts)

    partial_chrom_files = {}
    whole_chrom_files = {}
    for fn in os.listdir(path):
      pth = os.path.join(path, fn)
      if os.path.isfile(pth):
        base, ext = os.path.splitext(pth)
        ext = __trim_extension_dot(ext)
        if extensions is None or ext in extensions:
          idx_path = __find_index(pth, index_exts)
          if idx_path is None and fail_no_index:
            raise PyokitIOError("No index file for " + fn)
          chrom, start, end = __split_genomic_interval_filename(fn)
          assert((start is None and end is None) or
                 (start is not None and end is not None))
          if start is None:
            if chrom in whole_chrom_files:
              raise PyokitIOError("multiple files for chrom " + chrom)
            whole_chrom_files[chrom] = (pth, idx_path)
          else:
            k = (chrom, start, end)
            if k in partial_chrom_files:
              other_pth = partial_chrom_files[k][0]
              raise PyokitIOError("multiple files for " + str(k) + " --> " +
                                  pth + " and " + other_pth)
            partial_chrom_files[k] = (pth, idx_path)

    def factory(k):
      pth, idx = k
      return build_genome_alignment_from_file(pth, ref_spec, idx, verbose)

    return JustInTimeGenomeAlignment(whole_chrom_files, partial_chrom_files,
                                     factory)
164845
def parse_genome_results(self, f):
    
    regexes = {
        'bam_file': r"bam file = (.+)",
        'total_reads': r"number of reads = ([\d,]+)",
        'mapped_reads': r"number of mapped reads = ([\d,]+)",
        'mapped_bases': r"number of mapped bases = ([\d,]+)",
        'sequenced_bases': r"number of sequenced bases = ([\d,]+)",
        'mean_insert_size': r"mean insert size = ([\d,\.]+)",
        'median_insert_size': r"median insert size = ([\d,\.]+)",
        'mean_mapping_quality': r"mean mapping quality = ([\d,\.]+)",
        'general_error_rate': r"general error rate = ([\d,\.]+)",
    }
    d = dict()
    for k, r in regexes.items():
        r_search = re.search(r, f['f'], re.MULTILINE)
        if r_search:
            try:
                d[k] = float(r_search.group(1).replace(',',''))
            except ValueError:
                d[k] = r_search.group(1)
    # Check we have an input filename
    if 'bam_file' not in d:
        log.debug("Couldn't find an input filename in genome_results file {}".format(f['fn']))
        return None

    # Get a nice sample name
    s_name = self.clean_s_name(d['bam_file'], f['root'])

    # Add to general stats table & calculate a nice % aligned
    try:
        self.general_stats_data[s_name]['total_reads'] = d['total_reads']
        self.general_stats_data[s_name]['mapped_reads'] = d['mapped_reads']
        d['percentage_aligned'] = (d['mapped_reads'] / d['total_reads'])*100
        self.general_stats_data[s_name]['percentage_aligned'] = d['percentage_aligned']
        self.general_stats_data[s_name]['general_error_rate'] = d['general_error_rate']*100
    except KeyError:
        pass

    # Save results
    if s_name in self.qualimap_bamqc_genome_results:
        log.debug("Duplicate genome results sample name found! Overwriting: {}".format(s_name))
    self.qualimap_bamqc_genome_results[s_name] = d
    self.add_data_source(f, s_name=s_name, section='genome_results')
1349029
def parse_limits_list(path, limits):
    

    context = FilePathContext(path)

    for limit_def in limits:
        if 'include' in limit_def:
            include_context = context.resolve(limit_def['include'])
            for limit in parse_limits_file(include_context):
                yield limit
        else:
            yield parse_limit(limit_def)
848953
def process_associations(self, limit):
        

        myfile = '/'.join((self.rawdir, self.files['data']['file']))
        f = gzip.open(myfile, 'rb')
        filereader = io.TextIOWrapper(f, newline="")
        filereader.readline()  # remove the xml declaration line
        for event, elem in ET.iterparse(filereader):  # iterparse is not deprecated
            self.process_xml_table(
                elem, 'Article_Breed', self._process_article_breed_row, limit)
            self.process_xml_table(
                elem, 'Article_Phene', self._process_article_phene_row, limit)
            self.process_xml_table(
                elem, 'Breed_Phene', self._process_breed_phene_row, limit)
            self.process_xml_table(
                elem, 'Lida_Links', self._process_lida_links_row, limit)
            self.process_xml_table(
                elem, 'Phene_Gene', self._process_phene_gene_row, limit)
            self.process_xml_table(
                elem, 'Group_MPO', self._process_group_mpo_row, limit)
        f.close()
        return
161104
def read_10xgenomics(cls, tarball_fpath: str, prefix: str,
                         use_ensembl_ids: bool = False):
        """"""

        _LOGGER.info('Reading file: %s', tarball_fpath)

        with tarfile.open(tarball_fpath, mode='r:gz') as tf:
            ti = tf.getmember('%smatrix.mtx' % prefix)
            with tf.extractfile(ti) as fh:
                mtx = scipy.io.mmread(fh)

            ti = tf.getmember('%sgenes.tsv' % prefix)
            with tf.extractfile(ti) as fh:
                wrapper = io.TextIOWrapper(fh, encoding='ascii')
                i = 1
                if use_ensembl_ids:
                    i = 0
                gene_names = \
                        [row[i] for row in csv.reader(wrapper, delimiter='\t')]

            ti = tf.getmember('%sbarcodes.tsv' % prefix)
            with tf.extractfile(ti) as fh:
                wrapper = io.TextIOWrapper(fh, encoding='ascii')
                barcodes = \
                        [row[0] for row in csv.reader(wrapper, delimiter='\t')]


            assert mtx.shape[0] == len(gene_names)
            assert mtx.shape[1] == len(barcodes)
        
        _LOGGER.info('Matrix dimensions: %s', str(mtx.shape))
        X = mtx.todense()
        matrix = cls(X=X, genes=gene_names, cells=barcodes)
        
        return matrix
873587
def _process_pathway_pubmed(self, limit):
        
        LOG.info("Processing KEGG pathways to pubmed ids")
        if self.test_mode:
            graph = self.testgraph
        else:
            graph = self.graph
        line_counter = 0
        raw = '/'.join((self.rawdir, self.files['pathway_pubmed']['file']))
        with open(raw, 'r', encoding="iso-8859-1") as csvfile:
            filereader = csv.reader(csvfile, delimiter='\t', quotechar='\"')
            for row in filereader:
                line_counter += 1
                (pubmed_id, kegg_pathway_num) = row

                if self.test_mode and kegg_pathway_num not in self.test_ids['pathway']:
                    continue

                pubmed_id = pubmed_id.upper()
                # will look like KEGG-path:map04130
                kegg_id = 'KEGG-' + kegg_pathway_num

                r = Reference(graph, pubmed_id, self.globaltt['journal article'])
                r.addRefToGraph()
                graph.addTriple(pubmed_id, self.globaltt['is_about'], kegg_id)

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

        return