<jupyter_start><jupyter_text>Cell type differences and effects of interferon stimulation on immune cellsDemonstrating differential expression between cell types and the effect of interferon stimulation within a cell type (CD4 T cells).<jupyter_code>import pandas as pd
import matplotlib.pyplot as plt
import scanpy as sc
import scipy as sp
import itertools
import numpy as np
import scipy.stats as stats
from scipy.integrate import dblquad
import seaborn as sns
from statsmodels.stats.multitest import fdrcorrection
import imp
pd.options.display.max_rows = 999
pd.set_option('display.max_colwidth', -1)
import pickle as pkl
import time
#plt.style.use('ggplot')
%matplotlib inline
import matplotlib as mpl
mpl.rcParams['pdf.fonttype'] = 42
mpl.rcParams['ps.fonttype'] = 42

import matplotlib.pylab as pylab
params = {'legend.fontsize': 'medium',
         'axes.labelsize': 'medium',
         'axes.titlesize':'medium',
         'figure.titlesize':'medium',
         'xtick.labelsize':'small',
         'ytick.labelsize':'small'}
pylab.rcParams.update(params)

import warnings
warnings.filterwarnings('ignore')
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import sys
sys.path.append('/data/home/Github/scrna-parameter-estimation/schypo')
import estimator, simulate, schypo, bootstrap, util
data_path = '/data/parameter_estimation/'
fig_path = '/data/home/Github/scrna-parameter-estimation/figures/fig3/'<jupyter_output><empty_output><jupyter_text>Perform some basic filtering<jupyter_code>adata = sc.read(data_path + 'interferon_filtered.h5ad')
adata = adata[(adata.obs.multiplets == 'singlet') & (adata.obs.cell != 'nan'), :].copy()
adata.obs.cell.value_counts()
adata.obs['cell_type'] = (adata.obs['cell'].astype(str) + ' - ' + adata.obs['stim'].astype(str)).astype('category')
sc.pl.tsne(adata, color='cell_type')
plt.hist(np.log(adata[adata.obs.cell_type == 'CD14+ Monocytes - ctrl'].X.sum(axis=1).A1))
x = adata[adata.obs.cell_type == 'CD14+ Monocytes - ctrl'].X.mean(axis=0).A1
plt.hist(np.log(x[x > 1e-10]))<jupyter_output><empty_output><jupyter_text>Cell sizes differences<jupyter_code>data = adata.obs.copy()
data['gene'] = adata.X.sum(axis=1)
sns.boxplot(x='gene', y='cell_type', data=data)<jupyter_output><empty_output><jupyter_text>Number of unique expression pairs for random pairs of genes<jupyter_code>mat = adata.X.tocsc()
num_cells = (mat.toarray() > 0).sum(axis=0)
num_uniques = []
temp = []
num_trial = 100000
for trial in range(num_trial):

    idxs = np.random.choice(mat.shape[1], 2)
    if num_cells[idxs[0]] < 100 or num_cells[idxs[1]] < 100:
        continue
    pair_data = mat[:, idxs].toarray()

    num_unique = []

    for num_cell in np.logspace(2, 4.45, 10):

        subsample = pair_data[np.random.randint(0, mat.shape[0], size=int(num_cell)),]
        num_unique.append(np.unique(subsample).shape[0])

        del subsample
    num_uniques.append(num_unique)

np.percentile([i[-1] for i in num_uniques], 90)
fig = plt.figure(figsize=(5, 3))
plt.gca().xaxis.set_major_locator(plt.MultipleLocator(8000))
plt.xlabel('Number of cells')
plt.ylabel('Number of unique \npairs of expression')
for trial in range(len(num_uniques)):
    plt.plot(list(np.logspace(2, 4.45, 10)), num_uniques[trial], '-o')
plt.savefig(fig_path + 'ifn_unique_pairs.png', bbox_inches='tight', dpi=1200)<jupyter_output><empty_output><jupyter_text>Number of nonzero entries as a function of number of cells<jupyter_code>mat = adata.X#.toarray()
num_unique = []
for num_cell in np.logspace(2, 4.45, 10):
    print(num_cell)
    
    subsample = mat[np.random.randint(0, mat.shape[0], size=int(num_cell)), ].toarray()
    subsample_nonzero = subsample[subsample != 0]
    num_unique.append(np.unique(subsample_nonzero).shape[0]+1)
    
    del subsample
    del subsample_nonzero
fig = plt.figure(figsize=(2, 2))
plt.plot(list(np.logspace(2, 4.45, 10)), num_unique, '-o', color='k')
plt.gca().xaxis.set_major_locator(plt.MultipleLocator(8000))
plt.xlabel('Number of cells')
plt.ylabel('Number of unique \nexpression values')
plt.savefig(fig_path + 'all_ifn_unique_vals.pdf', bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>Further dropout<jupyter_code>def compare_esimators(q, plot=False, true_data=None, var_q=1e-10):
    
    q_sq = var_q + q**2
    
    true_data = simulate.simulate_transcriptomes(1000, 1000, correlated=True) if true_data is None else true_data
    true_relative_data = true_data / true_data.sum(axis=1).reshape(-1, 1)

    qs, captured_data = simulate.capture_sampling(true_data, q, q_sq)
    Nr = captured_data.sum(axis=1).mean()
    captured_relative_data = captured_data/captured_data.sum(axis=1).reshape(-1, 1)
    adata = sc.AnnData(sp.sparse.csr_matrix(captured_data))
    sf = estimator._estimate_size_factor(adata.X)

    good_idx = (captured_data.mean(axis=0) > q)

    # True moments
    m_true, v_true, corr_true = true_relative_data.mean(axis=0), true_relative_data.var(axis=0), np.corrcoef(true_relative_data, rowvar=False)

    # Compute 1D moments
    m_obs, v_obs = captured_relative_data.mean(axis=0), captured_relative_data.var(axis=0)
    m_poi, v_poi = estimator._poisson_1d(adata.X, size_factor=sf, n_obs=true_data.shape[0])
    m_hyp, v_hyp = estimator._hyper_1d(adata.X, size_factor=sf, n_obs=true_data.shape[0], q=q)

    # Compute 2D moments
    corr_obs = np.corrcoef(captured_relative_data, rowvar=False)
    cov_poi = estimator._poisson_cov(adata.X, n_obs=true_data.shape[0], size_factor=sf, idx1=None, idx2=None)
    cov_hyp = estimator._hyper_cov(adata.X, n_obs=true_data.shape[0], size_factor=sf, idx1=None, idx2=None, q=q)
    corr_poi = estimator._corr_from_cov(cov_poi, v_poi, v_poi)
    corr_hyp = estimator._corr_from_cov(cov_hyp, v_hyp, v_hyp)

    corr_poi[np.abs(corr_poi) > 1] = np.nan
    corr_hyp[np.abs(corr_hyp) > 1] = np.nan

    mean_list = [m_obs, m_poi, m_hyp]
    var_list = [v_obs, v_poi, v_hyp]
    corr_list = [corr_obs, corr_poi, corr_hyp]
    estimated_list = [mean_list, var_list, corr_list]
    true_list = [m_true, v_true, corr_true]

    if plot:
        count = 0
        for j in range(3):
            for i in range(3):

                plt.subplot(3, 3, count+1)


                if i != 2:
                    plt.scatter(
                        np.log(true_list[i][good_idx]),
                        np.log(estimated_list[i][j][good_idx]),
                        s=0.1)
                    plt.plot(np.log(true_list[i][good_idx]), np.log(true_list[i][good_idx]), linestyle='--', color='m')
                    plt.xlim(np.log(true_list[i][good_idx]).min(), np.log(true_list[i][good_idx]).max())
                    plt.ylim(np.log(true_list[i][good_idx]).min(), np.log(true_list[i][good_idx]).max())

                else:

                    x = true_list[i][good_idx, :][:, good_idx].ravel()
                    y = estimated_list[i][j][good_idx, :][:, good_idx].ravel()

                    sample_idx = np.random.choice(x.shape[0], 500000)

                    plt.scatter(
                        x[sample_idx],
                        y[sample_idx],
                        s=0.1)
                    plt.plot([-1, 1], [-1, 1],linestyle='--', color='m')
                    plt.xlim(-1, 1);
                    plt.ylim(-1, 1);
                    
                if i != 0:
                    plt.yticks([])
                    
                if i == 1 or i == 0:
                    
                    print((np.log(true_list[i][good_idx]) > np.log(estimated_list[i][j][good_idx])).mean())

                count += 1
    else:
        return qs, good_idx, estimated_list, true_list
adata.var
overall_mean = adata.X.mean(axis=0).A1
true_data =adata[:, overall_mean > 0.2].X.toarray()
print(true_data.shape)
q = 0.3
plt.figure(figsize=(4, 4))
plt.subplots_adjust(wspace=0.1, hspace=0.3)
compare_esimators(q, plot=True, true_data=true_data.astype(int))
# plt.savefig(figpath + 'poi_vs_hyper_scatter_1.png', bbox_inches='tight')<jupyter_output>0.5025380710659898
0.0
0.5025380710659898
0.9898477157360406
0.5025380710659898
0.5025380710659898
{'fraction_non_alphanumeric': 0.14282187384273545, 'fraction_numerical': 0.03209480311072707, 'mean_word_length': 3.598183881952327, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '177777', 'license': 'MIT', 'n_tokens_mistral': 3229, 'n_tokens_neox': 2944, 'n_words': 514, 'path': 'analysis/interferon_stimulation/interferon_simple_stats.ipynb', 'repo_name': 'yelabucsf/scrna-parameter-estimation', 'size': '8101'}