Estimates velocities in a gene-specific manner

    Arguments
    ---------
    data: :class:`~anndata.AnnData`
        Annotated data matrix.

    Returns
    -------
    Returns or updates `adata`
764215
def recover_dynamics(data, var_names='all', max_iter=100, learning_rate=None, add_key='fit', t_max=None, use_raw=False,
                     min_loss=True, fix_scaling=None, load_pars=None, return_model=False, plot_results=False, copy=False, **kwargs):
    
    adata = data.copy() if copy else data
    logg.info('recovering dynamics', r=True)

    if var_names is 'all':
        var_names = adata.var_names.tolist()
    else:
        var_names = make_unique_list(var_names, allow_array=True)
        var_names = [name for name in var_names if name in adata.var_names]

    if 'velocity_genes' in var_names and 'velocity_genes' in adata.var.keys():
        var_names = [name for name in var_names if adata[:, name].var.velocity_genes.values]

    alpha, beta, gamma, t_, scaling = read_pars(adata)
    idx = []
    L, P, T = [], [], adata.layers['fit_t'] if 'fit_t' in adata.layers.keys() else np.zeros(adata.shape) * np.nan

    for i, gene in enumerate(var_names):
        dm = DynamicsRecovery(adata, gene, use_raw=use_raw, load_pars=load_pars, fix_scaling=fix_scaling)
        if max_iter > 1:
            dm.fit(max_iter, learning_rate, **kwargs)

        ix = np.where(adata.var_names == gene)[0][0]
        idx.append(ix)

        alpha[ix], beta[ix], gamma[ix], t_[ix], scaling[ix] = dm.pars[:, np.argmin(dm.loss) if min_loss else -1]
        T[:, ix] = dm.t
        L.append(dm.loss)
        if plot_results and i < 4:
            P.append(dm.pars)

    T_max = np.percentile(T, 95, axis=0) - np.percentile(T, 5, axis=0)
    m = t_max / T_max if t_max is not None else np.ones(adata.n_vars)
    alpha, beta, gamma, T, t_ = alpha / m, beta / m, gamma / m, T * m, t_ * m

    write_pars(adata, [alpha, beta, gamma, t_, scaling])
    adata.layers['fit_t'] = T

    cur_len = adata.varm['loss'].shape[1] if 'loss' in adata.varm.keys() else 2
    max_len = max(np.max([len(l) for l in L]), cur_len)
    loss = np.ones((adata.n_vars, max_len)) * np.nan

    if 'loss' in adata.varm.keys():
        loss[:, :cur_len] = adata.varm['loss']

    loss[idx] = np.vstack([np.concatenate([l, np.ones(max_len-len(l)) * np.nan]) for l in L])
    adata.varm['loss'] = loss

    logg.info('    finished', time=True, end=' ' if settings.verbosity > 2 else '\n')
    logg.hint('added \n' 
              '    \'' + add_key + '_pars' + '\', fitted parameters for splicing dynamics (adata.var)')

    if plot_results:  # Plot Parameter Stats
        figsize = [12, 5]  # rcParams['figure.figsize']
        fontsize = rcParams['font.size']
        fig, axes = pl.subplots(nrows=len(var_names[:4]), ncols=6, figsize=figsize)
        pl.subplots_adjust(wspace=0.7, hspace=0.5)
        for i, gene in enumerate(var_names[:4]):
            P[i] *= np.array([1 / m[idx[i]], 1 / m[idx[i]], 1 / m[idx[i]], m[idx[i]], 1])[:, None]
            for j, pij in enumerate(P[i]):
                axes[i][j].plot(pij)
            axes[i][len(P[i])].plot(L[i])
            if i == 0:
                for j, name in enumerate(['alpha', 'beta', 'gamma', 't_', 'scaling', 'loss']):
                    axes[i][j].set_title(name, fontsize=fontsize)

    return dm if return_model else adata if copy else None
683988
def gffatributes2ids(s):
    
    Name,gene_id,transcript_id,protein_id,exon_id=np.nan,np.nan,np.nan,np.nan,np.nan
    if '=' in s:
        d=dict([i.split('=') for i in s.split(';')])
        if 'Parent' in d:
            d[d['Parent'].split(':')[0]+'_id']=d['Parent'].split(':')[1]
        Name,gene_id,transcript_id,protein_id,exon_id=np.nan,np.nan,np.nan,np.nan,np.nan
        if 'Name' in d:    
            Name=d['Name']
        if 'gene_id' in d:    
            gene_id=d['gene_id']
        if 'transcript_id' in d:    
            transcript_id=d['transcript_id']
        if 'protein_id' in d:    
            protein_id=d['protein_id']
        if 'exon_id' in d:    
            exon_id=d['exon_id']
    return Name,gene_id,transcript_id,protein_id,exon_id
63401
def update_data(self, data):
        
        oldn_dim = self.n_dim
        if data.dtype != np.float32:
            print("Warning: data was not float32. A 32-bit copy was made")
            self._data = np.float32(data)
        else:
            self._data = data
        self.n_vectors, self.n_dim = data.shape
        if self.n_dim != oldn_dim and oldn_dim != 0:
            raise Exception("The dimension of the new data does not match!")
        self.bmus = np.zeros(self.n_vectors * 2, dtype=np.intc)
914731
def runningstd(t, data, width):
    
    ne = len(t) - width
    t_new = np.zeros(ne)
    std_r = np.zeros(ne)
    for i in range(ne):
        t_new[i] = np.mean(t[i:i+width+1])
        std_r[i] = scipy.stats.nanstd(data[i:i+width+1])
    return t_new, std_r
325066
def dmxData(self, data: tuple):
        
        newData = [0]*512
        for i in range(0, min(len(data), 512)):
            newData[i] = data[i]
        self._dmxData = tuple(newData)
        # in theory this class supports dynamic length, so the next line is correcting the length
        self.length = 126 + len(self._dmxData)
172824
def update(self):
        '''
        
        '''
        observed_arr = self.noise_sampler.generate()
        inferenced_arr = self.inference(observed_arr)

        error_arr = self.__convolutional_auto_encoder.computable_loss.compute_loss(
            observed_arr,
            inferenced_arr
        )

        delta_arr = self.__convolutional_auto_encoder.computable_loss.compute_delta(
            observed_arr,
            inferenced_arr
        )

        delta_arr = self.__convolutional_auto_encoder.back_propagation(delta_arr)
        self.__convolutional_auto_encoder.optimize(self.__learning_rate, 1)

        return error_arr
622901
def _update_likelihood_model(self, inst, partition_parameters, tree):
        
        # Build transition matrix from dict
        model = partition_parameters['model']
        freqs = partition_parameters.get('frequencies')
        if model == 'LG':
            subs_model = phylo_utils.models.LG(freqs)
        elif model == 'WAG':
            subs_model = phylo_utils.models.WAG(freqs)
        elif model == 'GTR':
            rates = partition_parameters.get('rates')
            subs_model = phylo_utils.models.GTR(rates, freqs, True)
        else:
            raise ValueError("Can't handle this model: {}".format(model))
        tm = phylo_utils.markov.TransitionMatrix(subs_model)
        
        # Read alpha value
        alpha = partition_parameters['alpha']
        inst.set_tree(tree)
        inst.update_alpha(alpha)
        inst.update_transition_matrix(tm)
357234
def load_gene_exp_to_df(inst_path):
    '''
    
    '''

    import pandas as pd
    from scipy import io
    from scipy import sparse
    from ast import literal_eval as make_tuple

    # matrix
    Matrix = io.mmread( inst_path + 'matrix.mtx')
    mat = Matrix.todense()

    # genes
    filename = inst_path + 'genes.tsv'
    f = open(filename, 'r')
    lines = f.readlines()
    f.close()

    # # add unique id to all genes
    # genes = []
    # unique_id = 0
    # for inst_line in lines:
    #     inst_line = inst_line.strip().split()

    #     if len(inst_line) > 1:
    #       inst_gene = inst_line[1]
    #     else:
    #       inst_gene = inst_line[0]

    #     genes.append(inst_gene + '_' + str(unique_id))
    #     unique_id = unique_id + 1

    # add unique id only to duplicate genes
    ini_genes = []
    for inst_line in lines:
        inst_line = inst_line.strip().split()
        if len(inst_line) > 1:
          inst_gene = inst_line[1]
        else:
          inst_gene = inst_line[0]
        ini_genes.append(inst_gene)

    gene_name_count = pd.Series(ini_genes).value_counts()
    duplicate_genes = gene_name_count[gene_name_count > 1].index.tolist()

    dup_index = {}
    genes = []
    for inst_row in ini_genes:

      # add index to non-unique genes
      if inst_row in duplicate_genes:

        # calc_non-unque index
        if inst_row not in dup_index:
          dup_index[inst_row] = 1
        else:
          dup_index[inst_row] = dup_index[inst_row] + 1

        new_row = inst_row + '_' + str(dup_index[inst_row])

      else:
        new_row = inst_row

      genes.append(new_row)

    # barcodes
    filename = inst_path + 'barcodes.tsv'
    f = open(filename, 'r')
    lines = f.readlines()
    f.close()

    cell_barcodes = []
    for inst_bc in lines:
        inst_bc = inst_bc.strip().split('\t')

        # remove dash from barcodes if necessary
        if '-' in inst_bc[0]:
          inst_bc[0] = inst_bc[0].split('-')[0]

        cell_barcodes.append(inst_bc[0])

    # parse tuples if necessary
    try:
        cell_barcodes = [make_tuple(x) for x in cell_barcodes]
    except:
        pass

    try:
        genes = [make_tuple(x) for x in genes]
    except:
        pass

    # make dataframe
    df = pd.DataFrame(mat, index=genes, columns=cell_barcodes)

    return df
98860
def _get_data_versions(data):
    
    genome_dir = install.get_genome_dir(data["genome_build"], data["dirs"].get("galaxy"), data)
    if genome_dir:
        version_file = os.path.join(genome_dir, "versions.csv")
        if version_file and os.path.exists(version_file):
            return version_file
    return None
912641
def train(self, traindata: np.ndarray) -> None:
        
        self.clf.fit(traindata[:, 1:5], traindata[:, 5])
325134
def run(_, data, out_dir):
    
    out = []
    vcinfo = get_active_vcinfo(data)
    if vcinfo:
        if dd.get_phenotype(data) == "normal" and "germline" in vcinfo:
            out.append(_bcftools_stats(data, out_dir, "germline", germline=True))
        elif dd.get_phenotype(data) != "germline":
            out.append(_bcftools_stats(data, out_dir))
            if "germline" in vcinfo:
                out.append(_bcftools_stats(data, out_dir, "germline", germline=True))
        else:
            out.append(_bcftools_stats(data, out_dir, germline=True))
    out.append(_snpeff_stats(data, out_dir))

    out = [item for item in out if item]
    if out:
        return {"base": out[0], "secondary": out[1:]}
158052
def run(data):
    
    hlas = []
    for hla_fq in tz.get_in(["hla", "fastq"], data, []):
        hla_type = re.search("[.-](?P<hlatype>HLA-[\w-]+).fq", hla_fq).group("hlatype")
        if hla_type in SUPPORTED_HLAS:
            if utils.file_exists(hla_fq):
                hlas.append((hla_type, hla_fq))
    if len(hlas) > 0:
        out_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data), "align",
                                                  dd.get_sample_name(data), "hla",
                                                  "OptiType-HLA-A_B_C"))
        # When running UMIs and hla typing we want to pick the original fastqs
        if len(hlas) > len(SUPPORTED_HLAS):
            hlas = [x for x in hlas if os.path.basename(x[1]).find("-cumi") == -1]
        if len(hlas) == len(SUPPORTED_HLAS):
            hla_fq = combine_hla_fqs(hlas, out_dir + "-input.fq", data)
            if utils.file_exists(hla_fq):
                out_file = glob.glob(os.path.join(out_dir, "*", "*_result.tsv"))
                if len(out_file) > 0:
                    out_file = out_file[0]
                else:
                    out_file = _call_hla(hla_fq, out_dir, data)
                out_file = _prepare_calls(out_file, os.path.dirname(out_dir), data)
                data["hla"].update({"call_file": out_file,
                                    "hlacaller": "optitype"})
    return data
217730
def save_sparse_data(self, fname):
        
        data = self.adata_raw.X.T
        if data.getformat()=='csr':
            data=data.tocsc()

        cell_names = np.array(list(self.adata_raw.obs_names))
        gene_names = np.array(list(self.adata_raw.var_names))
        
        pickle.dump((data, cell_names, gene_names), open(fname, 'wb'))
745894
def matrix(fasta_path: 'path to tictax annotated fasta input',
           scafstats_path: 'path to BBMap scaftstats file'):
    '''
    
    '''
    records = SeqIO.parse(fasta_path, 'fasta')
    df = tictax.matrix(records, scafstats_path)
    df.to_csv(sys.stdout)
678364
def npar(self):
        
        self.control_data.npar = self.parameter_data.shape[0]
        return self.control_data.npar
240081
def fasta(self):
        
        logging.info('Extracting FASTA sequences matching PointFinder database')
        for sample in self.runmetadata.samples:
            # Ensure that there are sequence data to extract from the GenObject
            if GenObject.isattr(sample[self.analysistype], 'sequences'):
                # Set the name of the FASTA file
                sample[self.analysistype].pointfinderfasta = \
                    os.path.join(sample[self.analysistype].outputdir,
                                 '{seqid}_pointfinder.fasta'.format(seqid=sample.name))
                # Create a list to store all the SeqRecords created
                sequences = list()
                with open(sample[self.analysistype].pointfinderfasta, 'w') as fasta:
                    for gene, sequence in sample[self.analysistype].sequences.items():
                        # Create a SeqRecord using a Seq() of the sequence - both SeqRecord and Seq are from BioPython
                        seq = SeqRecord(seq=Seq(sequence),
                                        id=gene,
                                        name=str(),
                                        description=str())
                        sequences.append(seq)
                    # Write all the SeqRecords to file
                    SeqIO.write(sequences, fasta, 'fasta')
667372
def find_nb_genes(data):
    
    data_means = data.mean(1)
    data_vars = data.var(1)
    nb_indices = data_means < 0.9*data_vars
    return nb_indices