Read Gene Association File (GAF). Return data.
332656
def read_gaf(fin_gaf, prt=sys.stdout, hdr_only=False, allow_missing_symbol=False, **kws):
    """"""
    return GafReader(fin_gaf, hdr_only,
                     prt=prt, allow_missing_symbol=allow_missing_symbol).get_id2gos(**kws)
1113020
def genotype_to_df(g, samples, as_string=False):
    
    name = g.variant.name if g.variant.name else "genotypes"
    df = pd.DataFrame(g.genotypes, index=samples, columns=[name])

    if as_string:
        df["alleles"] = None

        hard_calls = df[name].round()
        df.loc[hard_calls == 0, "alleles"] = "{0}/{0}".format(g.reference)
        df.loc[hard_calls == 1, "alleles"] = "{0}/{1}".format(g.reference,
                                                              g.coded)
        df.loc[hard_calls == 2, "alleles"] = "{0}/{0}".format(g.coded)

        df = df[["alleles"]]
        df.columns = [name]

    return df
1132454
def _read_datafile_entry(spg, no, symbol, setting, f):
    """"""
    spg._no = no
    spg._symbol = symbol.strip()
    spg._setting = setting
    spg._centrosymmetric = bool(int(f.readline().split()[1]))
    # primitive vectors
    f.readline()
    spg._scaled_primitive_cell = np.array([
        list(map(float, f.readline().split())) 
        for i in range(3)],
        dtype='float')
    
    # primitive reciprocal vectors
    f.readline()
    spg._reciprocal_cell = np.array([list(map(int, f.readline().split())) 
                                        for i in range(3)],
                                       dtype=np.int)
    # subtranslations
    spg._nsubtrans = int(f.readline().split()[0])
    spg._subtrans = np.array([list(map(float, f.readline().split())) 
                              for i in range(spg._nsubtrans)],
                             dtype=np.float)
    # symmetry operations
    nsym = int(f.readline().split()[0])
    symop = np.array([list(map(float, f.readline().split())) for i in range(nsym)],
                     dtype=np.float)
    spg._nsymop = nsym
    spg._rotations = np.array(symop[:,:9].reshape((nsym,3,3)), dtype=np.int)
    spg._translations = symop[:,9:]
396515
def df_kegg_metadata(self):
        """"""
        kegg_pre_df = []
        df_cols = ['gene', 'kegg', 'refseq', 'uniprot', 'num_pdbs', 'pdbs', 'seq_len', 'sequence_file', 'metadata_file']

        for g in self.genes:
            kegg_mappings = g.protein.filter_sequences(KEGGProp)
            for kegg_prop in kegg_mappings:
                kegg_dict = kegg_prop.get_dict(df_format=True, only_attributes=df_cols)
                kegg_dict['gene'] = g.id
                kegg_pre_df.append(kegg_dict)

        # Save a dataframe of the file mapping info
        df = pd.DataFrame.from_records(kegg_pre_df, columns=df_cols).set_index('gene')
        if df.empty:
            log.warning('Empty dataframe')
            return df
        else:
            return ssbio.utils.clean_df(df)
1205984
def read_sparse(cls, file_path: str):
        """"""
        
        with open(file_path) as fh:
            next(fh)  # skip header line
            genes = next(fh)[1:-1].split('\t')
            cells = next(fh)[1:-1].split('\t')
            next(fh)
            m, n, nnz = [int(s) for s in next(fh)[:-1].split(' ')]
        
        t = pd.read_csv(file_path, sep=' ', skiprows=5, header=None,
                        dtype={0: np.uint32, 1: np.uint32})
        
        i = t[0].values - 1
        j = t[1].values - 1
        data = t[2].values

        assert data.size == nnz

        X = sparse.coo_matrix((data, (i,j)), shape=[m, n]).todense()

        return cls(X=X, genes=genes, cells=cells)
1203411
def read_datafiles(files, dtype, column):
    '''
    '''
    pha = []
    pha_fpi = []
    for filename, filetype in zip(files, dtype):
        if filetype == 'cov':
            cov = load_cov(filename)
        elif filetype == 'mag':
            mag = load_rho(filename, column)
        elif filetype == 'pha':
            pha = load_rho(filename, 2)
        elif filetype == 'pha_fpi':
            pha_fpi = load_rho(filename, 2)

    return cov, mag, pha, pha_fpi
237564
def tx2genefile(gtf, out_file=None):
    
    installed_tx2gene = os.path.join(os.path.dirname(gtf), "tx2gene.csv")
    if file_exists(installed_tx2gene):
        return installed_tx2gene
    if file_exists(out_file):
        return out_file
    with file_transaction(out_file) as tx_out_file:
        with open(tx_out_file, "w") as out_handle:
            for k, v in transcript_to_gene(gtf).items():
                out_handle.write(",".join([k, v]) + "\n")
    return out_file
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)
479847
def get_gene2gos(go2nt):
        """"""
        gene2gos = cx.defaultdict(set)
        nt0 = next(iter(go2nt.values()))
        b_str = isinstance(nt0.study_items, str)
        # print("NNNNTTTT", nt0)
        for goid, ntgo in go2nt.items():
            study_items = ntgo.study_items.split(', ') if b_str else ntgo.study_items
            for geneid in study_items:
                gene2gos[geneid].add(goid)
        if b_str:
          b_set = set(isinstance(g.isdigit(), int) for g in nt0.study_items.split(', '))
          if b_set == set([True]):
            return {int(g):gos for g, gos in gene2gos.items()}
        return {g:gos for g, gos in gene2gos.items()}
634517
def rnaseq2ga(quantificationFilename, sqlFilename, localName, rnaType,
              dataset=None, featureType="gene",
              description="", programs="", featureSetNames="",
              readGroupSetNames="", biosampleId=""):
    
    readGroupSetName = ""
    if readGroupSetNames:
        readGroupSetName = readGroupSetNames.strip().split(",")[0]
    featureSetIds = ""
    readGroupIds = ""
    if dataset:
        featureSetIdList = []
        if featureSetNames:
            for annotationName in featureSetNames.split(","):
                featureSet = dataset.getFeatureSetByName(annotationName)
                featureSetIdList.append(featureSet.getId())
            featureSetIds = ",".join(featureSetIdList)
        # TODO: multiple readGroupSets
        if readGroupSetName:
            readGroupSet = dataset.getReadGroupSetByName(readGroupSetName)
            readGroupIds = ",".join(
                [x.getId() for x in readGroupSet.getReadGroups()])
    if rnaType not in SUPPORTED_RNA_INPUT_FORMATS:
        raise exceptions.UnsupportedFormatException(rnaType)
    rnaDB = RnaSqliteStore(sqlFilename)
    if rnaType == "cufflinks":
        writer = CufflinksWriter(rnaDB, featureType, dataset=dataset)
    elif rnaType == "kallisto":
        writer = KallistoWriter(rnaDB, featureType, dataset=dataset)
    elif rnaType == "rsem":
        writer = RsemWriter(rnaDB, featureType, dataset=dataset)
    writeRnaseqTable(rnaDB, [localName], description, featureSetIds,
                     readGroupId=readGroupIds, programs=programs,
                     biosampleId=biosampleId)
    writeExpressionTable(writer, [(localName, quantificationFilename)])
192114
def read_cdf(self, infile):
        """"""
        if not self.return_cdf:
            self.load_cdf_module()

        if self.cdf_module == "scipy":
            # making it compatible to older scipy versions
            file_obj = self.cdf.netcdf_file(infile, mode="r")
        elif self.cdf_module == "netcdf4":
            file_obj = self.cdf.Dataset(infile)
        else:
            raise ImportError(
                "Could not import data \
                              from file {0}".format(
                    infile
                )
            )
        return file_obj
299884
def _load_gtf_as_dataframe(self, usecols=None, features=None):
        
        logger.info("Reading GTF from %s", self.gtf_path)
        df = read_gtf(
            self.gtf_path,
            column_converters={
                "seqname": normalize_chromosome,
                "strand": normalize_strand,
            },
            infer_biotype_column=True,
            usecols=usecols,
            features=features)

        column_names = set(df.keys())
        expect_gene_feature = features is None or "gene" in features
        expect_transcript_feature = features is None or "transcript" in features
        observed_features = set(df["feature"])

        # older Ensembl releases don't have "gene" or "transcript"
        # features, so fill in those rows if they're missing
        if expect_gene_feature and "gene" not in observed_features:
            # if we have to reconstruct gene feature rows then
            # fill in values for 'gene_name' and 'gene_biotype'
            # but only if they're actually present in the GTF
            logger.info("Creating missing gene features...")
            df = create_missing_features(
                dataframe=df,
                unique_keys={"gene": "gene_id"},
                extra_columns={
                    "gene": {
                        "gene_name",
                        "gene_biotype"
                    }.intersection(column_names),
                },
                missing_value="")
            logger.info("Done.")

        if expect_transcript_feature and "transcript" not in observed_features:
            logger.info("Creating missing transcript features...")
            df = create_missing_features(
                dataframe=df,
                unique_keys={"transcript": "transcript_id"},
                extra_columns={
                    "transcript": {
                        "gene_id",
                        "gene_name",
                        "gene_biotype",
                        "transcript_name",
                        "transcript_biotype",
                        "protein_id",
                    }.intersection(column_names)
                },
                missing_value="")
            logger.info("Done.")

        return df
1088047
def get_assoc_ncbi_taxids(taxids, force_dnld=False, loading_bar=True, **kws):
    """"""
    fin = kws['gene2go'] if 'gene2go' in kws else os.path.join(os.getcwd(), "gene2go")
    dnld_ncbi_gene_file(fin, force_dnld, loading_bar=loading_bar)
    return read_ncbi_gene2go(fin, taxids, **kws)
347100
def _biotype_lookup_fn(gtf):
    
    db = get_gtf_db(gtf)
    sources = set([feature.source for feature in db.all_features()])
    gene_biotypes = set([feature.attributes.get("gene_biotype", [None])[0]
                         for feature in db.all_features()])
    biotypes = set([feature.attributes.get("biotype", [None])[0]
                    for feature in db.all_features()])
    if "protein_coding" in sources:
        return lambda feature: feature.source
    elif "protein_coding" in biotypes:
        return lambda feature: feature.attributes.get("biotype", [None])[0]
    elif "protein_coding" in gene_biotypes:
        return lambda feature: feature.attributes.get("gene_biotype", [None])[0]
    else:
        return None
665948
def _chk_gene2go(assoc):
    """"""
    if not assoc:
        raise RuntimeError("NO ITEMS FOUND IN ASSOCIATIONS {A}".format(A=assoc))
    for key in assoc:
        if isinstance(key, str) and key[:3] == "GO:":
            raise Exception("ASSOCIATIONS EXPECTED TO BE gene2go, NOT go2gene: {EX}".format(
                EX=assoc.items()[:2]))
        return
117509
def read_ma_array(self, infile, var_name):
        """"""
        file_obj = self.read_cdf(infile)

        # .data is not backwards compatible to old scipy versions, [:] is
        data = file_obj.variables[var_name][:]

        # load numpy if available
        try:
            import numpy as np
        except Exception:
            raise ImportError("numpy is required to return masked arrays.")

        if hasattr(file_obj.variables[var_name], "_FillValue"):
            # return masked array
            fill_val = file_obj.variables[var_name]._FillValue
            retval = np.ma.masked_where(data == fill_val, data)
        else:
            # generate dummy mask which is always valid
            retval = np.ma.array(data)

        return retval
54099
def _importGenomeObjects(gtfFilePath, chroSet, genome, batchSize, verbose = 0) :
    """"""

    class Store(object) :
        
        def __init__(self, conf) :
            self.conf = conf
            self.chromosomes = {}
            
            self.genes = {}
            self.transcripts = {}
            self.proteins = {}
            self.exons = {}

        def batch_save(self) :
            self.conf.db.beginTransaction()
            
            for c in self.genes.itervalues() :
                c.save()
                conf.removeFromDBRegistery(c)
                
            for c in self.transcripts.itervalues() :
                c.save()
                conf.removeFromDBRegistery(c.exons)
                conf.removeFromDBRegistery(c)
            
            for c in self.proteins.itervalues() :
                c.save()
                conf.removeFromDBRegistery(c)
            
            self.conf.db.endTransaction()
            
            del(self.genes)
            del(self.transcripts)
            del(self.proteins)
            del(self.exons)
            
            self.genes = {}
            self.transcripts = {}
            self.proteins = {}
            self.exons = {}

            gc.collect()

        def save_chros(self) :
            pBar = ProgressBar(nbEpochs = len(self.chromosomes))
            for c in self.chromosomes.itervalues() :
                pBar.update(label = 'Chr %s' % c.number)
                c.save()
            pBar.close()
        
    printf('Importing gene set infos from %s...' % gtfFilePath)
    
    printf('Backuping indexes...')
    indexes = conf.db.getIndexes()
    printf("Droping all your indexes, (don't worry i'll restore them later)...")
    Genome_Raba.flushIndexes()
    Chromosome_Raba.flushIndexes()
    Gene_Raba.flushIndexes()
    Transcript_Raba.flushIndexes()
    Protein_Raba.flushIndexes()
    Exon_Raba.flushIndexes()
    
    printf("Parsing gene set...")
    gtf = GTFFile(gtfFilePath, gziped = True)
    printf('Done. Importation begins!')
    
    store = Store(conf)
    chroNumber = None
    pBar = ProgressBar(nbEpochs = len(gtf))
    for line in gtf :
        chroN = line['seqname']
        pBar.update(label = "Chr %s" % chroN)
        
        if (chroN.upper() in chroSet or chroN.lower() in chroSet):
            strand = line['strand']
            gene_biotype = line['gene_biotype']
            regionType = line['feature']
            frame = line['frame']

            start = int(line['start']) - 1
            end = int(line['end'])
            if start > end :
                start, end = end, start

            chroNumber = chroN.upper()
            if chroNumber not in store.chromosomes :
                store.chromosomes[chroNumber] = Chromosome_Raba()
                store.chromosomes[chroNumber].set(genome = genome, number = chroNumber)
            
            try :
                geneId = line['gene_id']
                geneName =  line['gene_name']
            except KeyError :
                geneId = None
                geneName = None
                if verbose :
                    printf('Warning: no gene_id/name found in line %s' % gtf[i])
            
            if geneId is not None :
                if geneId not in store.genes :
                    if len(store.genes) > batchSize :
                        store.batch_save()
                    
                    if verbose > 0 :
                        printf('\tGene %s, %s...' % (geneId, geneName))
                    store.genes[geneId] = Gene_Raba()
                    store.genes[geneId].set(genome = genome, id = geneId, chromosome = store.chromosomes[chroNumber], name = geneName, strand = strand, biotype = gene_biotype)
                if start < store.genes[geneId].start or store.genes[geneId].start is None :
                        store.genes[geneId].start = start
                if end > store.genes[geneId].end or store.genes[geneId].end is None :
                    store.genes[geneId].end = end
            try :
                transId = line['transcript_id']
                transName = line['transcript_name']
                try :
                    transcript_biotype = line['transcript_biotype']
                except KeyError :
                    transcript_biotype = None
            except KeyError :
                transId = None
                transName = None
                if verbose > 2 :
                    printf('\t\tWarning: no transcript_id, name found in line %s' % gtf[i])
            
            if transId is not None :
                if transId not in store.transcripts :
                    if verbose > 1 :
                        printf('\t\tTranscript %s, %s...' % (transId, transName))
                    store.transcripts[transId] = Transcript_Raba()
                    store.transcripts[transId].set(genome = genome, id = transId, chromosome = store.chromosomes[chroNumber], gene = store.genes.get(geneId, None), name = transName, biotype=transcript_biotype)
                if start < store.transcripts[transId].start or store.transcripts[transId].start is None:
                    store.transcripts[transId].start = start
                if end > store.transcripts[transId].end or store.transcripts[transId].end is None:
                    store.transcripts[transId].end = end
            
                try :
                    protId = line['protein_id']
                except KeyError :
                    protId = None
                    if verbose > 2 :
                        printf('Warning: no protein_id found in line %s' % gtf[i])

                # Store selenocysteine positions in transcript
                if regionType == 'Selenocysteine':
                    store.transcripts[transId].selenocysteine.append(start)
                        
                if protId is not None and protId not in store.proteins :
                    if verbose > 1 :
                        printf('\t\tProtein %s...' % (protId))
                    store.proteins[protId] = Protein_Raba()
                    store.proteins[protId].set(genome = genome, id = protId, chromosome = store.chromosomes[chroNumber], gene = store.genes.get(geneId, None), transcript = store.transcripts.get(transId, None), name = transName)
                    store.transcripts[transId].protein = store.proteins[protId]

                try :
                    exonNumber = int(line['exon_number']) - 1
                    exonKey = (transId, exonNumber)
                except KeyError :
                    exonNumber = None
                    exonKey = None
                    if verbose > 2 :
                        printf('Warning: no exon number or id found in line %s' % gtf[i])
                
                if exonKey is not None :
                    if verbose > 3 :
                        printf('\t\t\texon %s...' % (exonId))
                    
                    if exonKey not in store.exons and regionType == 'exon' :
                        store.exons[exonKey] = Exon_Raba()
                        store.exons[exonKey].set(genome = genome, chromosome = store.chromosomes[chroNumber], gene = store.genes.get(geneId, None), transcript = store.transcripts.get(transId, None), protein = store.proteins.get(protId, None), strand = strand, number = exonNumber, start = start, end = end)
                        store.transcripts[transId].exons.append(store.exons[exonKey])
                    
                    try :
                        store.exons[exonKey].id = line['exon_id']
                    except KeyError :
                        pass
                    
                    if regionType == 'exon' :
                        if start < store.exons[exonKey].start or store.exons[exonKey].start is None:
                            store.exons[exonKey].start = start
                        if end > store.transcripts[transId].end or store.exons[exonKey].end is None:
                            store.exons[exonKey].end = end
                    elif regionType == 'CDS' :
                        store.exons[exonKey].CDS_start = start
                        store.exons[exonKey].CDS_end = end
                        store.exons[exonKey].frame = frame
                    elif regionType == 'stop_codon' :
                        if strand == '+' :
                            if store.exons[exonKey].CDS_end != None :
                                store.exons[exonKey].CDS_end += 3
                                if store.exons[exonKey].end < store.exons[exonKey].CDS_end :
                                    store.exons[exonKey].end = store.exons[exonKey].CDS_end
                                if store.transcripts[transId].end < store.exons[exonKey].CDS_end :
                                    store.transcripts[transId].end = store.exons[exonKey].CDS_end
                                if store.genes[geneId].end < store.exons[exonKey].CDS_end :
                                    store.genes[geneId].end = store.exons[exonKey].CDS_end
                        if strand == '-' :
                            if store.exons[exonKey].CDS_start != None :
                                store.exons[exonKey].CDS_start -= 3
                                if store.exons[exonKey].start > store.exons[exonKey].CDS_start :
                                    store.exons[exonKey].start = store.exons[exonKey].CDS_start
                                if store.transcripts[transId].start > store.exons[exonKey].CDS_start :
                                    store.transcripts[transId].start = store.exons[exonKey].CDS_start
                                if store.genes[geneId].start > store.exons[exonKey].CDS_start :
                                    store.genes[geneId].start = store.exons[exonKey].CDS_start
    pBar.close()
    
    store.batch_save()
    
    conf.db.beginTransaction()
    printf('almost done saving chromosomes...')
    store.save_chros()
    
    printf('saving genome object...')
    genome.save()
    conf.db.endTransaction()
    
    conf.db.beginTransaction()
    printf('restoring core indexes...')
    # Genome.ensureGlobalIndex(('name', 'species'))
    # Chromosome.ensureGlobalIndex('genome')
    # Gene.ensureGlobalIndex('genome')
    # Transcript.ensureGlobalIndex('genome')
    # Protein.ensureGlobalIndex('genome')
    # Exon.ensureGlobalIndex('genome')
    Transcript.ensureGlobalIndex('exons')
    
    printf('commiting changes...')
    conf.db.endTransaction()
    
    conf.db.beginTransaction()
    printf('restoring user indexes')
    pBar = ProgressBar(label = "restoring", nbEpochs = len(indexes))
    for idx in indexes :
        pBar.update()
        conf.db.execute(idx[-1].replace('CREATE INDEX', 'CREATE INDEX IF NOT EXISTS'))
    pBar.close()
    
    printf('commiting changes...')
    conf.db.endTransaction()
    
    return store.chromosomes.values()