In [1]:

%pylab inline
from glob import glob

import anndata
import pandas as pd
from scipy import optimize
from scipy.special import gammaln
from scipy.special import psi
from scipy.special import factorial
from scipy.optimize import fmin_l_bfgs_b as optim
from tqdm import tqdm

Populating the interactive namespace from numpy and matplotlib

In [2]:

datasets =  [anndata.read(fname) for fname in glob('../Data/output/*.h5ad')]

In [3]:

datasets[4].var.head().T

Out[3]:

.dataframe tbody tr th:only-of-type { vertical-align: middle; } .dataframe tbody tr th { vertical-align: top; } .dataframe thead th { text-align: right; } index ERCC-00002 ERCC-00003 ERCC-00004 ERCC-00009 ERCC-00012 empirical_mean 1.796595e+03 2.336190e+02 6.908095e+02 1.139643e+02 0.047619 empirical_variance 1.753958e+05 6.245618e+03 3.797782e+04 1.969701e+03 0.116780 empirical_zero_fraction 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.976190 ml_mean 1.796595e+03 2.336190e+02 6.908095e+02 1.139643e+02 0.047412 genewise_dispersion 5.378328e-02 1.101545e-01 7.813416e-02 1.428824e-01 30.633224 global_zero_fraction 4.582850e-41 2.508432e-23 1.318163e-32 1.357982e-17 0.953549 genewise_zero_fraction 1.017812e-37 1.108822e-13 5.350741e-23 2.178272e-09 0.971059 scaled_count_mean 1.478226e-01 1.922201e-02 5.683932e-02 9.376901e-03 0.000004 poisson_zero_fraction 0.000000e+00 7.780822e-39 0.000000e+00 3.329074e-20 0.953562

In [4]:

for adata in datasets:
    mean_of_scaled = np.array(adata.X.sum(0) / adata.X.sum())
    adata.var['scaled_count_mean'] = mean_of_scaled[0]
    adata.obs['total_counts'] = np.array(adata.X.sum(1))

In [5]:

def prob_zero_fun(mu, counts):
    return np.exp(-(mu[:,np.newaxis]).dot(counts.values[np.newaxis])).sum(1) / len(counts)

In [6]:

datasets[4].var.head().T

Out[6]:

.dataframe tbody tr th:only-of-type { vertical-align: middle; } .dataframe tbody tr th { vertical-align: top; } .dataframe thead th { text-align: right; } index ERCC-00002 ERCC-00003 ERCC-00004 ERCC-00009 ERCC-00012 empirical_mean 1.796595e+03 2.336190e+02 6.908095e+02 1.139643e+02 0.047619 empirical_variance 1.753958e+05 6.245618e+03 3.797782e+04 1.969701e+03 0.116780 empirical_zero_fraction 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.976190 ml_mean 1.796595e+03 2.336190e+02 6.908095e+02 1.139643e+02 0.047412 genewise_dispersion 5.378328e-02 1.101545e-01 7.813416e-02 1.428824e-01 30.633224 global_zero_fraction 4.582850e-41 2.508432e-23 1.318163e-32 1.357982e-17 0.953549 genewise_zero_fraction 1.017812e-37 1.108822e-13 5.350741e-23 2.178272e-09 0.971059 scaled_count_mean 1.478226e-01 1.922201e-02 5.683932e-02 9.376901e-03 0.000004 poisson_zero_fraction 0.000000e+00 7.780822e-39 0.000000e+00 3.329074e-20 0.953562

In [7]:

for adata in datasets:
    adata.var['poisson_zero_fraction'] = prob_zero_fun(adata.var['scaled_count_mean'],
                                                       adata.obs['total_counts'])

In [8]:

for adata in datasets:
    adata.write('../Data/output/' + adata.uns['name'] + '.h5ad')

In [ ]: