Vireo

https://vireoSNP.readthedocs.io/en/latest/
Admin
Vireo is a Bayesian method to demultiplex pooled single-cell RNA-seq data with or without genotype...
Tokens:19,695
Snippets:121
Trust Score:9.6
Update:1 month ago
Show doc for...
Tokens:
Raw
### Install VireoSNP from source

Source: https://vireosnp.readthedocs.io/en/latest/install.html

Download the source code from GitHub and run the setup script.

```bash
wget https://github.com/huangyh09/vireo/archive/master.zip
unzip master.zip
cd vireo-master

python setup.py install
```

--------------------------------

### Verify VireoSNP installation

Source: https://vireosnp.readthedocs.io/en/latest/install.html

Run the command to check if the package is correctly installed and accessible.

```bash
Welcome to vireoSNP v0.1.1!

use -h or --help for help on argument.
```

--------------------------------

### Install VireoSNP via PyPI

Source: https://vireosnp.readthedocs.io/en/latest/install.html

Use pip to install or upgrade the package from the Python Package Index.

```bash
pip install vireoSNP

pip install --upgrade --no-deps vireoSNP
```

--------------------------------

### Demultiplexing Strategies

Source: https://vireosnp.readthedocs.io/en/latest/manual.html

Command-line examples for running Vireo under different genotype availability scenarios.

```bash
vireo -c $CELL_DATA -N $n_donor -o $OUT_DIR
```

```bash
vireo -c $CELL_DATA -d $DONOR_GT_FILE -o $OUT_DIR
```

```bash
vireo -c $CELL_DATA -d $DONOR_GT_FILE -o $OUT_DIR -N $n_donor
```

```bash
vireo -c $CELL_DATA -d $DONOR_GT_FILE -o $OUT_DIR --forceLearnGT
```

--------------------------------

### Install VireoSNP via git

Source: https://vireosnp.readthedocs.io/en/latest/install.html

Install the latest version directly from the GitHub repository using pip.

```bash
pip install -U git+https://github.com/single-cell-genetics/vireo
```

--------------------------------

### Import vireoSNP and dependencies

Source: https://vireosnp.readthedocs.io/en/latest/vireoSNP_clones-mkn45.html

Initializes the environment by importing necessary libraries and checking the installed vireoSNP version.

```python
import vireoSNP
import numpy as np
from scipy import sparse
from scipy.io import mmread
import matplotlib.pyplot as plt

print(vireoSNP.__version__)
```

--------------------------------

### Troubleshoot command not found error

Source: https://vireosnp.readthedocs.io/en/latest/install.html

Example of the error message displayed when the executable is not in the system PATH.

```bash
vireo: command not found
```

--------------------------------

### Fit Model with Multiple Initializations

Source: https://vireosnp.readthedocs.io/en/latest/_modules/vireoSNP/utils/bmm_model.html

Runs the fitting process multiple times to find the optimal parameters, handling sparse matrix conversion if necessary.

```python
    def fit(self, AD, DP, n_init=10, max_iter=200, max_iter_pre=100, 
        random_seed=None, **kwargs):
        """Fit VB with multiple initializations

        Parameters
        ----------
        AD : scipy.sparse.csc_matrix (n_var, n_cell)
            Sparse count matrix for alternative allele
        DP : scipy.sparse.csc_matrix (n_var, n_cell)
            Sparse count matrix for depths, alternative + refeerence alleles
        n_inits : int
            Number of random initialisations to use
        max_iter : int
            Maximum number of iterations for _fit_BV() in best initial
        max_iter_pre : int
            Maximum number of iterations for _fit_BV() in multiple initials
        min_iter :
            Minimum number of iterations for _fit_BV()
        epsilon_conv : float
            Threshold for detecting convergence for _fit_BV()
        verbose : bool
            Whether print out log info for _fit_BV()
        random_seed : None or int
            Random seed in numpy.random for multiple initializations
        """
        if random_seed is not None:
            np.random.seed(random_seed)
            
        if type(DP) is np.ndarray and np.mean(DP > 0) < 0.3:
            print("Warning: input matrices is %.1f%% sparse, " 
                  %(100 - np.mean(DP > 0) * 100) +
                  "change to scipy.sparse.csc_matrix" )
            AD = csc_matrix(AD)
            DP = csc_matrix(DP)

        _binom_coeff = np.sum(get_binom_coeff(AD, DP, is_log = True))

        self.ELBO_inits = []
        for i in range(n_init):
            self.set_initial(
                self.beta_mu_init, self.beta_sum_init, self.ID_prob_init
            )
            self._fit_BV(AD, DP, max_iter=max_iter_pre, **kwargs)
            self.ELBO_inits.append(self.ELBO_iters[-1])
            
            ## first or better initialization
            if i == 0 or (self.ELBO_iters[-1] > np.max(self.ELBO_inits[:-1])):
                _ID_prob_best = self.ID_prob + 0
                _beta_mu_best = self.beta_mu + 0
                _beta_sum_best = self.beta_sum + 0
                _ELBO_iters_best = self.ELBO_iters + 0
                
        ## Re-fit with best parameters
        self.set_initial(_beta_mu_best, _beta_sum_best, _ID_prob_best)
        self.ELBO_iters = _ELBO_iters_best
        self._fit_BV(AD, DP, max_iter=max_iter, **kwargs)

        ## add binomial coefficient constants
        self.ELBO_iters = self.ELBO_iters + _binom_coeff
        self.ELBO_inits = np.array(self.ELBO_inits) + _binom_coeff
```

--------------------------------

### Vireo.__init__

Source: https://vireosnp.readthedocs.io/en/latest/API.html

Initializes the Vireo model with the specified dimensions and configuration parameters.

```APIDOC
## Vireo.__init__

### Description
Initializes the Vireo model for variational inference. Multiple initializations are recommended to avoid local optima.

### Parameters
- **n_cell** (int) - Required - Number of cells
- **n_var** (int) - Required - Number of variants
- **n_donor** (int) - Required - Number of donors
- **n_GT** (int) - Optional - Number of genotype categories (default: 3)
- **learn_GT** (bool) - Optional - Whether to update GT_prob (default: True)
- **ASE_mode** (bool) - Optional - Whether to set allelic ratio theta to be variant specific (default: False)
- **fix_beta_sum** (bool) - Optional - Whether to fix the concentration parameter of theta's posterior (default: False)
- **beta_mu_init** (numpy array) - Optional - Initial value of beta_mu
- **beta_sum_init** (numpy array) - Optional - Initial value of beta_sum
- **ID_prob_init** (numpy array) - Optional - Initial value of ID_prob
- **GT_prob_init** (numpy array) - Optional - Initial value of GT_prob
```

--------------------------------

### VireoBulk.__init__

Source: https://vireosnp.readthedocs.io/en/latest/_modules/vireoSNP/utils/vireo_bulk.html

Initializes the VireoBulk model with the number of donors and genotype categories.

```APIDOC
## VireoBulk.__init__(n_donor, n_GT=3, psi_init=None, theta_init=[0.01, 0.5, 0.99])

### Description
Initializes the model parameters for donor abundance estimation.

### Parameters
- **n_donor** (int) - Required - Number of donors in the mixture.
- **n_GT** (int) - Optional - Number of genotype categories (default: 3).
- **psi_init** (numpy.array) - Optional - Initial fractional abundance of each donor.
- **theta_init** (list) - Optional - Initial alternative allele rate in each genotype category.
```

--------------------------------

### Initialize Vireo Model

Source: https://vireosnp.readthedocs.io/en/latest/_modules/vireoSNP/utils/vireo_model.html

The constructor initializes the model parameters. Multiple initializations are recommended to avoid local optima.

```python
def __init__(self, n_cell, n_var, n_donor, n_GT=3, learn_GT=True,  
    learn_theta=True, ASE_mode=False, fix_beta_sum=False, 
    beta_mu_init=None, beta_sum_init=None, ID_prob_init=None, 
    GT_prob_init=None):
    """Initialise Vireo model

    Note, multiple initializations are highly recomended to avoid local 
    optima.
    
    Parameters
    ----------
    n_cell : int. 
        Number of cells
    n_var : int. 
        Number of variants
    n_donor : int. 
        Number of donors
    n_GT : int. 
        Number of genotype categories
    learn_GT: bool. 
        Whether updating `GT_prob`; otherwise using the initial
    ASE_mode: bool. 
        Whether setting allelic ratio `theta` to be variant specific
    fix_beta_sum: bool. 
        Whether fixing the concetration parameter of theta's posterior
    beta_mu_init: numpy array (1, n_GT) or (n_var, n_GT)
        Initial value of beta_mu, the mean parameter of theta
    beta_sum_init: numpy array (1, n_GT) or (n_var, n_GT), same as beta_mu
        Initial value of beta_sum, the concetration parameter of theta
    ID_prob_init: numpy array (n_cell, n_donor)
        Initial value of ID_prob, cell assignment probability to each donor
    GT_prob_init: numpy array (n_var, n_donor, n_GT)
        Initial value of GT_prob, genotype probability per variant and donor
    """
    self.n_GT = n_GT
    self.n_var = n_var
    self.n_cell = n_cell
    self.n_donor = n_donor
    self.learn_GT = learn_GT
    self.ASE_mode = ASE_mode
    self.learn_theta = learn_theta
    self.fix_beta_sum = fix_beta_sum
    
    self.ELBO_ = np.zeros((0))
    
    # initial key parameters
    self.set_initial(beta_mu_init, beta_sum_init, ID_prob_init, GT_prob_init)
    
    # set hyper parameters for prior
    self.set_prior()
```

--------------------------------

### Initialize and fit the BinomMixtureVB model

Source: https://vireosnp.readthedocs.io/en/latest/vireoSNP_clones-mkn45.html

Imports the model class, initializes it with the loaded data, and performs variational inference fitting.

```python
from vireoSNP import BinomMixtureVB
```

```python
_model = BinomMixtureVB(n_var=AD.shape[0], n_cell=AD.shape[1], n_donor=5)
_model.fit(AD, DP, min_iter=30, n_init=50)
print(_model.ELBO_iters[-1])
```

--------------------------------

### Vireo Help Command

Source: https://vireosnp.readthedocs.io/en/latest/manual.html

Displays the full list of available arguments and options for the Vireo command-line tool.

```bash
Usage: vireo [options]

Options:
-h, --help            show this help message and exit
-c CELL_DATA, --cellData=CELL_DATA
                        The cell genotype file in VCF format or cellSNP folder
                        with sparse matrices.
-N N_DONOR, --nDonor=N_DONOR
                        Number of donors to demultiplex; can be larger than
                        provided in donor_file
-o OUT_DIR, --outDir=OUT_DIR
                        Dirtectory for output files [default:
                        $cellFilePath/vireo]

Optional input files:
   --vartrixData=VARTRIX_DATA
                        The cell genotype files in vartrix outputs (three/four
                        files, comma separated):
                        alt.mtx,ref.mtx,barcodes.tsv,SNPs.vcf.gz. This will
                        suppress cellData argument.
   -d DONOR_FILE, --donorFile=DONOR_FILE
                        The donor genotype file in VCF format. Please filter
                        the sample and region with bcftools -s and -R first!
   -t GENO_TAG, --genoTag=GENO_TAG
                        The tag for donor genotype: GT, GP, PL [default: PL]

Optional arguments:
   --noDoublet         If use, not checking doublets.
   -M N_INIT, --nInit=N_INIT
                        Number of random initializations, when GT needs to
                        learn [default: 50]
   --extraDonor=N_EXTRA_DONOR
                        Number of extra donor in pre-cluster, when GT needs to
                        learn [default: 0]
   --extraDonorMode=EXTRA_DONOR_MODE
                        Method for searching from extra donors. size: n_cell
                        per donor; distance: GT distance between donors
                        [default: distance]
   --forceLearnGT      If use, treat donor GT as prior only.
   --ASEmode           If use, turn on SNP specific allelic ratio.
   --noPlot            If use, turn off plotting GT distance.
   --randSeed=RAND_SEED
                        Seed for random initialization [default: none]
   --cellRange=CELL_RANGE
                        Range of cells to process, eg. 0-10000 [default: all]
   --callAmbientRNAs   If use, detect ambient RNAs in each cell (under
                        development)
   -p NPROC, --nproc=NPROC
                        Number of subprocesses for computing - this sacrifices
                        memory for speedups [default: 1]
```

--------------------------------

### Set Initial Values

Source: https://vireosnp.readthedocs.io/en/latest/_modules/vireoSNP/utils/vireo_model.html

Configures initial values for beta parameters, cell assignment probabilities, and genotype probabilities.

```python
def set_initial(self, beta_mu_init=None, beta_sum_init=None, 
    ID_prob_init=None, GT_prob_init=None):
    """Set initial values
    """
    theta_len = self.n_var if self.ASE_mode else 1

    if beta_mu_init is not None:
        self.beta_mu = beta_mu_init
    else:
        self.beta_mu = (np.ones((theta_len, self.n_GT)) * 
            np.linspace(0.01, 0.99, self.n_GT).reshape(1, -1))

    if beta_sum_init is not None:
        self.beta_sum = beta_sum_init
    else:
        self.beta_sum = np.ones((theta_len, self.n_GT)) * 50

    if ID_prob_init is not None:
        self.ID_prob = normalize(ID_prob_init, axis=1)
    else:
        self.ID_prob = normalize(np.random.rand(self.n_cell, self.n_donor))

    if GT_prob_init is not None:
        self.GT_prob = normalize(GT_prob_init)
    else:        
        _GT_val = np.random.rand(self.n_var, self.n_donor, self.n_GT)
        self.GT_prob = normalize(_GT_val)
```

--------------------------------

### BinomMixtureVB.__init__

Source: https://vireosnp.readthedocs.io/en/latest/API.html

Initializes the BinomMixtureVB model with specified dimensions and optional initial parameters.

```APIDOC
## BinomMixtureVB(n_cell, n_var, n_donor, fix_beta_sum=False, beta_mu_init=None, beta_sum_init=None, ID_prob_init=None)

### Description
Initializes the Vireo model. Multiple initializations are highly recommended to avoid local optima.

### Parameters
- **n_cell** (int) - Required - Number of cells
- **n_var** (int) - Required - Number of variants
- **n_donor** (int) - Required - Number of donors
- **fix_beta_sum** (bool) - Optional - Whether fixing the concentration parameter of theta’s posterior
- **beta_mu_init** (numpy array) - Optional - Initial value of beta_mu
- **beta_sum_init** (numpy array) - Optional - Initial value of beta_sum
- **ID_prob_init** (numpy array) - Optional - Initial value of ID_prob
```

--------------------------------

### Add VireoSNP to PATH

Source: https://vireosnp.readthedocs.io/en/latest/install.html

Append the local binary directory to the PATH environment variable in your shell configuration file.

```bash
export PATH="~/.local/bin:$PATH"
```

--------------------------------

### vireo CLI

Source: https://vireosnp.readthedocs.io/en/latest/manual.html

The main command-line interface for running Vireo demultiplexing.

```APIDOC
## vireo

### Description
Runs the Vireo demultiplexing pipeline on cell genotype data.

### Usage
`vireo [options]`

### Parameters
- **-c, --cellData** (string) - Required - The cell genotype file in VCF format or cellSNP folder with sparse matrices.
- **-N, --nDonor** (integer) - Optional - Number of donors to demultiplex.
- **-o, --outDir** (string) - Optional - Directory for output files.
- **--vartrixData** (string) - Optional - Comma separated list of Vartrix output files (alt.mtx,ref.mtx,barcodes.tsv,SNPs.vcf.gz).
- **-d, --donorFile** (string) - Optional - The donor genotype file in VCF format.
- **-t, --genoTag** (string) - Optional - The tag for donor genotype: GT, GP, PL [default: PL].
- **--noDoublet** (flag) - Optional - If used, disables doublet checking.
- **-M, --nInit** (integer) - Optional - Number of random initializations [default: 50].
- **--forceLearnGT** (flag) - Optional - If used, treats donor GT as prior only.
- **-p, --nproc** (integer) - Optional - Number of subprocesses for computing [default: 1].
```

--------------------------------

### Initialize Model Priors

Source: https://vireosnp.readthedocs.io/en/latest/_modules/vireoSNP/utils/vireo_model.html

Initializes ID and GT priors, ensuring they are properly shaped and normalized.

```python
        if ID_prior is not None:
            if len(ID_prior.shape) == 1:
                ID_prior = np.expand_dims(ID_prior, axis=0)
            self.ID_prior = ID_prior
        else:
            self.ID_prior = normalize(np.ones(self.ID_prob.shape))

        if GT_prior is not None:
            if len(GT_prior.shape) == 2:
                GT_prior = np.expand_dims(GT_prior, axis=0)
            GT_prior[GT_prior < min_GP] = min_GP
            GT_prior[GT_prior > 1 - min_GP] = 1 - min_GP
            GT_prior = normalize(GT_prior)
            self.GT_prior = GT_prior
        else:        
            self.GT_prior = normalize(np.ones(self.GT_prob.shape))
```

--------------------------------

### Internal VB fitting method

Source: https://vireosnp.readthedocs.io/en/latest/_modules/vireoSNP/utils/vireo_model.html

The _fit_VB method performs coordinate ascent to fit the Vireo model, tracking the ELBO for convergence.

```python
def _fit_VB(self, AD, DP, max_iter=200, min_iter=5, epsilon_conv=1e-2,
    delay_fit_theta=0, verbose=True):
    """Fit Vireo model with coordinate ascent
    """
    ELBO = np.zeros(max_iter)
    numerical_minimal = 1e-6
    for it in range(max_iter):
        if self.learn_theta and it >= delay_fit_theta:
            self.update_theta_size(AD, DP)
        if self.learn_GT:
            self.update_GT_prob(AD, DP)

        _logLik_ID = self.update_ID_prob(AD, DP)
        ELBO[it] = self.get_ELBO(_logLik_ID) #+ _binom_coeff_log
        
        if it > min_iter:
            if ELBO[it] < ELBO[it - 1] - numerical_minimal:
                if verbose:
                    print("Warning: Lower bound decreases!\n")
            elif it == max_iter - 1:
                if verbose:
                    print("Warning: VB did not converge!\n")
            elif ELBO[it] - ELBO[it - 1] < epsilon_conv:
                break
    
    return ELBO[:it]
```

--------------------------------

### Verify Global Optima with Multiple Runs

Source: https://vireosnp.readthedocs.io/en/latest/vireoSNP_clones-mkn45.html

Run the model multiple times to check if the same ELBO is consistently achieved, indicating a likely global optimum.

```python
n_init = 50
for i in range(3):
    _model = BinomMixtureVB(n_var=AD.shape[0], n_cell=AD.shape[1], n_donor=3)
    _model.fit(AD, DP, min_iter=30, n_init=n_init)
    print("rerun %d:" %i, _model.ELBO_iters[-1])
```

--------------------------------

### Visualize model fitting diagnostics

Source: https://vireosnp.readthedocs.io/en/latest/vireoSNP_clones-mkn45.html

Plots the ELBO values across multiple initializations and iterations to assess model convergence.

```python
fig = plt.figure(figsize=(11, 4))
plt.subplot(1, 2, 1)
plt.hist(_model.ELBO_inits)
plt.ylabel("Frequency")
plt.xlabel("ELBO in multiple initializations")

plt.subplot(1, 2, 2)
plt.plot(_model.ELBO_iters)
plt.xlabel("Iterations")
plt.ylabel("ELBO in a single initialization")

plt.tight_layout()
plt.show()
```

--------------------------------

### vireoSNP.vcf.load_VCF

Source: https://vireosnp.readthedocs.io/en/latest/API.html

Loads a whole VCF file into matrices.

```APIDOC
## vireoSNP.vcf.load_VCF(vcf_file, biallelic_only=False, load_sample=True, sparse=True, format_list=None)

### Description
Load whole VCF file. Initially designed to load VCF from cellSNP output.

### Parameters
- **vcf_file** (str) - Required - Path to the VCF file.
- **biallelic_only** (bool) - Optional - Whether to keep only biallelic variants.
- **load_sample** (bool) - Optional - Whether to load sample IDs.
- **sparse** (bool) - Optional - Whether to use sparse matrices.
- **format_list** (list) - Optional - List of formats to extract (e.g., ['GT', 'AD', 'DP']).

### Returns
- A dictionary containing variants, samples, and GenoINFO.
```

--------------------------------

### Import vireoSNP

Source: https://vireosnp.readthedocs.io/en/latest/API.html

Standard import statement for the vireoSNP package.

```python
import vireoSNP
```

--------------------------------

### Compress probability files using gzip

Source: https://vireosnp.readthedocs.io/en/latest/_modules/vireoSNP/utils/io_utils.html

Executes a system gzip command to compress the generated probability TSV files.

```python
bashCommand = "gzip -f %s %s" %(out_dir + "/prob_singlet.tsv", 
    out_dir + "/prob_doublet.tsv")
pro = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE)
pro.communicate()[0]
```

--------------------------------

### Public fit method

Source: https://vireosnp.readthedocs.io/en/latest/_modules/vireoSNP/utils/vireo_model.html

The fit method serves as the primary entry point for fitting the Vireo model, handling input validation and sparse matrix conversion.

```python
def fit(self, AD, DP, max_iter=200, min_iter=5, epsilon_conv=1e-2,
    delay_fit_theta=0, verbose=True, n_inits=50, nproc=1):
    """Fit Vireo model with coordinate ascent

    Parameters
    ----------
    AD : scipy.sparse.csc_matrix (n_var, n_cell)
        Sparse count matrix for alternative allele
    DP : scipy.sparse.csc_matrix (n_var, n_cell)
        Sparse count matrix for depths, alternative + refeerence alleles
    max_iter : int
        Maximum number of iterations
    min_iter :
        Minimum number of iterations
    epsilon_conv : float
        Threshold for detecting convergence
    delay_fit_theta : int
        Number of steps to delay updating theta. This can be very useful 
        for common genetics when there is good prior on allelic ratio.
    verbose : bool
        Whether print out log info
    """
    if type(DP) is np.ndarray and np.mean(DP > 0) < 0.3:
        print("Warning: input matrices is %.1f%% sparse, " 
              %(100 - np.mean(DP > 0) * 100) +
              "change to scipy.sparse.csc_matrix" )
        AD = csc_matrix(AD)
        DP = csc_matrix(DP)

    ELBO = self._fit_VB(AD, DP, max_iter, min_iter, epsilon_conv, 
                        delay_fit_theta, verbose)

    # _binom_coeff_log = np.sum(logbincoeff(DP, AD, is_sparse=True))
    # _binom_coeff_log = np.sum(get_binom_coeff(AD, DP))
    
    ELBO += np.sum(get_binom_coeff(AD, DP))

    self.ELBO_ = np.append(self.ELBO_, ELBO)
```

--------------------------------

### Visualize assignment probability and allelic ratio

Source: https://vireosnp.readthedocs.io/en/latest/vireoSNP_clones-mkn45.html

Generates heatmaps for cell assignment probabilities and mean allelic ratios per clone.

```python
from vireoSNP.plot import heat_matrix

fig = plt.figure(figsize=(7, 4), dpi=100)
plt.subplot(1, 2, 1)
im = heat_matrix(_model.ID_prob, cmap="Blues", alpha=0.8,
                 display_value=False, row_sort=True)
plt.colorbar(im, fraction=0.046, pad=0.04)
plt.title("Assignment probability")
plt.xlabel("Clone")
plt.ylabel("%d cells" %(_model.n_cell))
plt.xticks(range(_model.n_donor))


plt.subplot(1, 2, 2)
im = heat_matrix(_model.beta_mu, cmap=segpink, alpha=0.8,
                 display_value=False, row_sort=True)
plt.colorbar(im, fraction=0.046, pad=0.04)
plt.title("Mean allelic ratio")
plt.xlabel("Clone")
plt.ylabel("%d SNPs" %(_model.n_var))
plt.xticks(range(_model.n_donor))

plt.tight_layout()
plt.show()
# plt.savefig("you_favorate_path with png or pdf")
```

--------------------------------

### Load VCF file

Source: https://vireosnp.readthedocs.io/en/latest/API.html

Loads a VCF file into matrices. Ensure the VCF is filtered with bcftools first to manage memory usage.

```python
>>> import vireoSNP
>>> import numpy as np
>>> vcf_dat = vireoSNP.vcf.load_VCF("cellSNP.cells.vcf.gz", sparse=False,
>>>     biallelic_only=False, format_list=['GT', 'AD', 'DP', 'ALL'])
>>> var_ids = np.array(vcf_dat['variants'])
>>> samples = np.array(vcf_dat['samples'])
>>> GT_mat = np.array(vcf_dat['GenoINFO']['GT'])
>>> AD_mat = np.array(vcf_dat['GenoINFO']['AD']).astype(float)
>>> DP_mat = np.array(vcf_dat['GenoINFO']['DP']).astype(float)
>>> ALL_bases_mat = np.array(vcf_dat['GenoINFO']['ALL'])
```

--------------------------------

### Generating Heatmap

Source: https://vireosnp.readthedocs.io/en/latest/vireoSNP_clones-mkn45.html

Creates a heatmap visualization of allele frequencies annotated by clone labels.

```python
vireoSNP.plot.anno_heat((AD/DP)[var_idx, :], col_anno=cell_label, col_order_ids=id_uniq,
                        cmap=segpink, yticklabels=mtSNP_ids[var_idx])
```

```python
<seaborn.matrix.ClusterGrid at 0x7428b2c8ab90>
```

--------------------------------

### Visualize ELBO Results

Source: https://vireosnp.readthedocs.io/en/latest/vireoSNP_clones-mkn45.html

Plot the maximum ELBO and distribution for different numbers of clones to assist in selection.

```python
plt.plot(np.arange(1, len(n_clone_list)+1), np.max(_ELBO_mat, axis=1))
plt.boxplot(_ELBO_mat)
plt.xticks(np.arange(1, len(n_clone_list)+1), n_clone_list)
plt.ylabel("ELBO")
plt.xlabel("n_clones")
plt.show()
```

--------------------------------

### Configure numpy print options

Source: https://vireosnp.readthedocs.io/en/latest/vireoSNP_clones-mkn45.html

Sets the formatting for numpy array printing.

```python
# np.set_printoptions(formatter={'float': lambda x: format(x, '.5f')})
```

--------------------------------

### Read VarTrix output

Source: https://vireosnp.readthedocs.io/en/latest/_modules/vireoSNP/utils/io_utils.html

Reads alternative and reference allele matrices from VarTrix output files.

```python
def read_vartrix(alt_mtx, ref_mtx, cell_file, vcf_file=None):
    """Read data from VarTrix

    Parameters
    ----------
    alt_mtx:
        sparse matrix file for alternative alleles
    ref_mtx:
        sparse matrix file for reference alleles
    cell_file:
        file for cell barcodes, each per line
    vcf_file:
        the vcf file used for fetch variants in VarTrix
    
    Return
    ------
    A disctionary containing AD, DP, cells and optionally variants
    """
    if vcf_file is not None:
        cell_dat = load_VCF(vcf_file, load_sample=False, biallelic_only=False)
        cell_dat['variants'] = np.array(cell_dat['variants'])
    else:
        cell_dat = {}
    cell_dat['AD'] = mmread(alt_mtx).tocsc()
    cell_dat['DP'] = mmread(ref_mtx).tocsc() + cell_dat['AD']
    cell_dat['samples'] = np.genfromtxt(cell_file, dtype=str)
    return cell_dat
```

--------------------------------

### Vireo.fit

Source: https://vireosnp.readthedocs.io/en/latest/API.html

Fits the Vireo model using coordinate ascent.

```APIDOC
## Vireo.fit

### Description
Fits the Vireo model to the provided count matrices using coordinate ascent.

### Parameters
- **AD** (scipy.sparse.csc_matrix) - Required - Sparse count matrix for alternative allele
- **DP** (scipy.sparse.csc_matrix) - Required - Sparse count matrix for depths
- **max_iter** (int) - Optional - Maximum number of iterations (default: 200)
- **min_iter** (int) - Optional - Minimum number of iterations (default: 5)
- **epsilon_conv** (float) - Optional - Threshold for detecting convergence (default: 0.01)
- **delay_fit_theta** (int) - Optional - Steps to delay updating theta (default: 0)
- **verbose** (bool) - Optional - Whether to print log info (default: True)
```

--------------------------------

### VireoBulk Class Definition and Initialization

Source: https://vireosnp.readthedocs.io/en/latest/_modules/vireoSNP/utils/vireo_bulk.html

The VireoBulk class estimates donor abundance in a multiplexed bulk sample. The constructor initializes the model with the number of donors and genotype categories.

```python
class VireoBulk():
    """
    Estimate of donor abundance in a multipexed bulk sample

    Varibale to infer
    -----------------
    psi: numpy.array (n_donor, )
        The fractional abundance of each donor in the mixture
    theta: numpy.array (n_GT, )
        The alternative allele rate in each genotype category

    Parameters
    ----------
    n_GT: int, number of genotype categories
    n_donor: int, number of donors in the mixture
    """


    def __init__(self, n_donor, n_GT=3, psi_init=None, 
                 theta_init=[0.01, 0.5, 0.99]):
        self.n_GT = n_GT
        self.n_donor = n_donor
        
        self.psi = np.random.dirichlet([1] * n_donor)
        self.theta = np.random.rand(n_GT)
        
        if psi_init is not None:
            if n_donor != len(psi_init):
                print("Warning: n_donor != len(psi_init)")
            else:
                self.psi = np.random.dirichlet([1] * n_donor)

        if theta_init is not None:
            if n_GT != len(theta_init):
                print("Warning: n_GT != len(theta_init)")
            else:
                self.theta = theta_init
```

--------------------------------

### BinomMixtureVB.fit

Source: https://vireosnp.readthedocs.io/en/latest/API.html

Fits the variational inference model using multiple initializations.

```APIDOC
## fit(AD, DP, n_init=10, max_iter=200, max_iter_pre=100, random_seed=None, **kwargs)

### Description
Fits the VB model with multiple initializations to find the best fit.

### Parameters
- **AD** (scipy.sparse.csc_matrix) - Required - Sparse count matrix for alternative allele
- **DP** (scipy.sparse.csc_matrix) - Required - Sparse count matrix for depths
- **n_inits** (int) - Optional - Number of random initialisations
- **max_iter** (int) - Optional - Maximum iterations for best initial
- **max_iter_pre** (int) - Optional - Maximum iterations for multiple initials
- **random_seed** (int/None) - Optional - Random seed for reproducibility
```

--------------------------------

### Define custom colormap for allelic ratio

Source: https://vireosnp.readthedocs.io/en/latest/vireoSNP_clones-mkn45.html

Creates a custom colormap to highlight specific allelic frequency ranges.

```python
from matplotlib import cm
from matplotlib.colors import ListedColormap, LinearSegmentedColormap

raw_col = cm.get_cmap('pink_r', 200)
new_col = np.vstack((raw_col(np.linspace(0, 0.7, 10)),
                     raw_col(np.linspace(0.7, 1, 90))))
segpink = ListedColormap(new_col, name='segpink')
```

--------------------------------

### Save ambient RNA profile to TSV

Source: https://vireosnp.readthedocs.io/en/latest/_modules/vireoSNP/utils/io_utils.html

Writes ambient RNA proportions and log-likelihood ratios to a TSV file if available.

```python
if res_vireo['ambient_Psi'] is not None:
    fid = open(out_dir + "/prop_ambient.tsv", "w")
    fid.writelines("\t".join(["cell"] + donor_names + ['logLik_ratio']) + "\n")
    for i in range(len(cell_names)):
        line = ["%.4e" %x for x in res_vireo['ambient_Psi'][i, :]]
        line += ['%.2f' %res_vireo['Psi_LLRatio'][i]]
        fid.writelines("\t".join([cell_names[i]] + line) + "\n")
    fid.close()
```

--------------------------------

### Load mitochondrial variant data

Source: https://vireosnp.readthedocs.io/en/latest/vireoSNP_clones-mkn45.html

Downloads and reads sparse matrices for allele depth (AD) and depth (DP) from a remote repository.

```python
import requests
from io import BytesIO

data_path = "https://github.com/aaronkwc/MQuad_paper_reproduced_results/raw/main/mkn45/data/mquad/"

AD = mmread(BytesIO(requests.get(data_path + "passed_ad.mtx").content)).tocsc()
DP = mmread(BytesIO(requests.get(data_path + "passed_dp.mtx").content)).tocsc()

mtSNP_ids = np.genfromtxt(data_path + "passed_variant_names.txt", dtype='str')

# AD = mmread("data/passed_ad.mtx").tocsc()
# DP = mmread("data/passed_dp.mtx").tocsc()
# mtSNP_ids = np.array(['mt_variant%d' %x for x in range(AD.shape[0])])
```

--------------------------------

### Calculate chi-squared statistic for bulk data

Source: https://vireosnp.readthedocs.io/en/latest/_modules/vireoSNP/utils/vireo_bulk.html

Computes the likelihood ratio test statistic and p-value based on genotype probabilities and allele counts.

```python
    from scipy.stats import chi2

    BD = DP - AD
    theta_vct_alt = np.dot(np.dot(GT_prob, theta), psi)
    logLik_alt = np.sum(
        AD * np.log(theta_vct_alt) + BD * np.log(1 - theta_vct_alt))

    theta_vct_null = np.dot(np.dot(GT_prob, theta), psi_null)
    logLik_null = np.sum(
        AD * np.log(theta_vct_null) + BD * np.log(1 - theta_vct_null))

    LR = 2 * (logLik_alt - logLik_null)
    df = len(psi_null) - 1
    if log:
        pval = chi2.logsf(LR, df)
    else:
        pval = chi2.sf(LR, df)
    
    return LR, pval
```

--------------------------------

### Binomial Mixture Model Implementation

Source: https://vireosnp.readthedocs.io/en/latest/_modules/vireoSNP/utils/bmm_model.html

The BinomMixtureVB class handles variational inference for binomial mixture models. It requires initialization of cell, variant, and donor counts, and supports custom priors and initial parameter values.

```python
import itertools
import numpy as np
from scipy.stats import entropy
from scipy.sparse import csc_matrix
from scipy.special import logsumexp, digamma, betaln
from .vireo_base import normalize, loglik_amplify, beta_entropy, get_binom_coeff



[docs]
class BinomMixtureVB():
    """Binomial mixture model with variational inference

    The prior can be set via set_prior() before fitting the model.

    Key properties
    --------------
    beta_mu: numpy array (n_var, n_donor)
        Beta mean parameter of theta's posterior
    beta_sum: numpy array (n_var, n_donor)
        Beta concetration parameter of theta's posterior
    ID_prob: numpy array (n_cell, n_donor)
        Posterior cell assignment probability to each donor
    """



[docs]
    def __init__(self, n_cell, n_var, n_donor, fix_beta_sum=False, 
        beta_mu_init=None, beta_sum_init=None, ID_prob_init=None):
        """Initialise Vireo model

        Note, multiple initializations are highly recomended to avoid local 
        optima.
        
        Parameters
        ----------
        n_cell : int. 
            Number of cells
        n_var : int. 
            Number of variants
        n_donor : int. 
            Number of donors
        fix_beta_sum: bool. 
            Whether fixing the concetration parameter of theta's posterior
        beta_mu_init: numpy array (n_var, n_donor)
            Initial value of beta_mu, the mean parameter of theta
        beta_sum_init: numpy array (n_var, n_donor)
            Initial value of beta_sum, the concetration parameter of theta
        ID_prob_init: numpy array (n_cell, n_donor)
            Initial value of ID_prob, cell assignment probability to each donor
        """
        self.n_var = n_var
        self.n_cell = n_cell
        self.n_donor = n_donor
        self.fix_beta_sum = fix_beta_sum
        self.ID_prob_init = ID_prob_init
        self.beta_mu_init = beta_mu_init
        self.beta_sum_init = beta_sum_init
        
        # set priors; you can re-set by run this function
        self.set_prior()

        # initial key parameters
        self.set_initial(
            self.beta_mu_init, self.beta_sum_init, self.ID_prob_init
        )


        

[docs]
    def set_initial(self, beta_mu_init=None, beta_sum_init=None, 
        ID_prob_init=None):
        """Random initialization
        """
        # initial key parameters
        if beta_mu_init is not None:
            self.beta_mu = beta_mu_init
        else:
            self.beta_mu = np.ones((self.n_var, self.n_donor)) * 0.5

        if beta_sum_init is not None:
            self.beta_sum = beta_sum_init
        else:
            self.beta_sum = np.ones(self.beta_mu.shape) * 30

        if ID_prob_init is not None:
            self.ID_prob = normalize(ID_prob_init, axis=1)
        else:
            self.ID_prob = normalize(np.random.rand(self.n_cell, self.n_donor))

        self.ELBO_iters = np.array([])


    

[docs]
    def set_prior(self, ID_prior=None, beta_mu_prior=None, 
        beta_sum_prior=None):
        """Set prior for key variables: theta and ID_prob.
        The priors are in the same shape as its according variables.
        """
        if beta_mu_prior is None:
            beta_mu_prior = np.ones((self.n_var, self.n_donor)) * 0.5
        if beta_sum_prior is None:
            beta_sum_prior = np.ones(beta_mu_prior.shape) * 2.0

        self.theta_s1_prior = beta_mu_prior * beta_sum_prior
        self.theta_s2_prior = (1 - beta_mu_prior) * beta_sum_prior

        if ID_prior is not None:
            if len(ID_prior.shape) == 1:
                ID_prior = np.expand_dims(ID_prior, axis=0)
            self.ID_prior = ID_prior
        else:
            self.ID_prior = normalize(np.ones((self.n_cell, self.n_donor)))



    @property
    def theta_s1(self):
        """Beta concetration1 parameter for theta posterior"""
        return self.beta_mu * self.beta_sum

    @property
    def theta_s2(self):
        """Beta concetration2 parameter for theta posterior"""
        return (1 - self.beta_mu) * self.beta_sum


    def get_E_logLik(self, AD, DP):
        """Get the expecation of logLikelihood
        E_theta [P(AD|DP, theta, Z)]
        """
        BD = DP - AD

        # shape: (n_cell, n_donor)
        _E_logLik_mat = (
            AD.T @ digamma(self.theta_s1) + 
            BD.T @ digamma(self.theta_s2) -
            DP.T @ digamma(self.theta_s1 + self.theta_s2)
        )
        return _E_logLik_mat


    def update_theta_size(self, AD, DP):
        """Coordinate ascent for updating theta posterior parameters
        """
        BD = DP - AD
        _theta_s1 = AD @ self.ID_prob  #(n_var, n_donor)
        _theta_s2 = BD @ self.ID_prob  #(n_var, n_donor)
        _theta_s1 += self.theta_s1_prior
        _theta_s2 += self.theta_s2_prior
```

--------------------------------

### Fitting VireoBulk Model with EM Algorithm

Source: https://vireosnp.readthedocs.io/en/latest/_modules/vireoSNP/utils/vireo_bulk.html

The fit method uses the Expectation-Maximization algorithm to estimate psi and theta parameters based on allele counts and genotype probabilities.

```python
def fit(self, AD, DP, GT_prob, max_iter=200, min_iter=5, epsilon_conv=1e-3,
        learn_theta=True, delay_fit_theta=0, model="EM", verbose=False):
    """Fit the unknown variable psi and theta with EM algorithm

    Parameters
    ----------
    AD: numpy.array, (n_variant, ), int
        The count vector for alternative allele in all variants
    DP: numpy.array (n_variant, ), int
        The count vector for depths in all variants (i.e., two alleles)
    GT_prob: numpy.array, (n_variants, n_donor, n_GT)
        The probability tensor for each genotype in each donor
    learn_theta: bool
        Whether learn theta, otherwise use theta_init
    delay_fit_theta: int
        The number of steps to delay in updating theta
    max_iter : int
        Maximum number of iterations
    min_iter :
        Minimum number of iterations
    epsilon_conv : float
        Threshold for detecting convergence
    model: string
        The algorithm used to fit the model. Only "EM" is supported for
        Expectation-Maximumization algorithm
    verbose : bool
        Whether print out log info
    """
    BD = DP - AD
    
    logLik = np.zeros(max_iter)
    for it in range(max_iter):
        # E step: expectation of count assignment probability
        theta_mat = np.tensordot(GT_prob, self.theta, axes=(2, 0)) #np.dot
        Z1 = theta_mat * np.expand_dims(self.psi, 0)
        Z1 = Z1 / np.sum(Z1, axis=1, keepdims=True)
        
        Z0 = (1 - theta_mat) * np.expand_dims(self.psi, 0)
        Z0 = Z0 / np.sum(Z0, axis=1, keepdims=True)

        # M step: maximize logLikehood over psi and theta
        psi_raw = np.dot(AD, Z1) + np.dot(BD, Z0)
        self.psi = psi_raw / np.sum(psi_raw)
        
        if learn_theta and it >= delay_fit_theta:
            theta_s1 = np.dot(AD, np.sum(GT_prob * np.expand_dims(Z1, 2), axis = 1))
            theta_s2 = np.dot(BD, np.sum(GT_prob * np.expand_dims(Z0, 2), axis = 1))
            self.theta = theta_s1 / (theta_s1 + theta_s2)
        
        # Likelihood and check convergence
        theta_vct = np.dot(np.dot(GT_prob, self.theta), self.psi)
        logLik[it] = np.sum(
            AD * np.log(theta_vct) + BD * np.log(1 - theta_vct))
        if it > min_iter:
            if logLik[it] < logLik[it - 1]:
                if verbose:
                    print("Warning: logLikelihood decreases!\n")
            elif it == max_iter - 1:
                if verbose:
                    print("Warning: VB did not converge!\n")
            elif logLik[it] - logLik[it - 1] < epsilon_conv:
                break
    
    self.logLik = logLik[it]
    self.logLik_all = logLik[:it]
```

--------------------------------

### Fit Model with Coordinate Ascent

Source: https://vireosnp.readthedocs.io/en/latest/_modules/vireoSNP/utils/bmm_model.html

Iteratively updates model parameters until convergence or maximum iterations are reached.

```python
    def _fit_BV(self, AD, DP, max_iter=200, min_iter=20, epsilon_conv=1e-2,
        verbose=True):
        """Fit Vireo model with coordinate ascent
        """
        ELBO = np.zeros(max_iter)
        for it in range(max_iter):
            self.update_theta_size(AD, DP)
            _logLik_ID = self.get_E_logLik(AD, DP)

            self.update_ID_prob(logLik_ID = _logLik_ID)
            ELBO[it] = self.get_ELBO(logLik_ID = _logLik_ID)

            if it > min_iter:
                if ELBO[it] - ELBO[it - 1] < -1e-6:
                    if verbose:
                        print("Warning: ELBO decreases %.8f to %.8f!\n"
                              %(ELBO[it - 1], ELBO[it]))
                elif it == max_iter - 1:
                    if verbose:
                        print("Warning: VB did not converge!\n")
                elif ELBO[it] - ELBO[it - 1] < epsilon_conv:
                    break

        self.ELBO_iters = np.append(self.ELBO_iters, ELBO[:it])
```

--------------------------------

### vireoSNP.vcf.match_VCF_samples

Source: https://vireosnp.readthedocs.io/en/latest/API.html

Matches donors in two VCF files.

```APIDOC
## vireoSNP.vcf.match_VCF_samples(VCF_file1, VCF_file2, GT_tag1, GT_tag2)

### Description
Match donors in two VCF files.

### Parameters
- **VCF_file1** (str) - Required - Full path of first VCF file.
- **VCF_file2** (str) - Required - Full path of second VCF file.
- **GT_tag1** (str) - Required - Tag for extracting genotype probability in VCF1.
- **GT_tag2** (str) - Required - Tag for extracting genotype probability in VCF2.
```

--------------------------------

### Evaluate ELBO for Different Clone Counts

Source: https://vireosnp.readthedocs.io/en/latest/vireoSNP_clones-mkn45.html

Iterate through a range of clone counts to collect ELBO values for model selection.

```python
n_init = 50
n_clone_list = np.arange(2, 10)

_ELBO_mat = []
for k in n_clone_list:
    _model = BinomMixtureVB(n_var=AD.shape[0], n_cell=AD.shape[1], n_donor=k)
    _model.fit(AD, DP, min_iter=30, n_init=n_init)
    _ELBO_mat.append(_model.ELBO_inits)
```
Tokens:
Raw