<jupyter_start><jupyter_text>Variability analysis for HBEC IFN experiment<jupyter_code>import scanpy as sc
import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats
from pybedtools import BedTool
import pickle as pkl
%matplotlib inline
import sys
sys.path.append('/data/home/Github/scrna-parameter-estimation/dist/memento-0.0.2-py3.7.egg')
sys.path.append('/data/home/Github/misc-seq/miscseq/')
import encode
import memento
data_path = '/data_volume/ifn_hbec/'<jupyter_output><empty_output><jupyter_text>Get genes that are likely ISGF3 binding<jupyter_code>enc = encode.Encode('/data/home/Github/misc-seq/miscseq/GRCh38Genes.bed')
stat1_bed = enc.get_encode_peaks_union([
    'https://www.encodeproject.org/files/ENCFF747ICD/@@download/ENCFF747ICD.bed.gz',
    'https://www.encodeproject.org/files/ENCFF646MXG/@@download/ENCFF646MXG.bed.gz'
])
stat2_bed = enc.get_encode_peaks_union([
    'https://www.encodeproject.org/files/ENCFF204VQS/@@download/ENCFF204VQS.bed.gz',
    'https://www.encodeproject.org/files/ENCFF545JVD/@@download/ENCFF545JVD.bed.gz'
])
intersect_bed = stat1_bed.intersect(stat2_bed).sort()
isgf3_genes = enc.get_peak_genes_bed(intersect_bed, 0).query('distance == 0 & ~gene.str.contains("MT").values').gene.tolist()<jupyter_output>/data/home/anaconda3/envs/single_cell/lib/python3.7/site-packages/pybedtools/bedtool.py:3681: UserWarning: Default names for filetype bed are:
['chrom', 'start', 'end', 'name', 'score', 'strand', 'thickStart', 'thickEnd', 'itemRgb', 'blockCount', 'blockSizes', 'blockStarts']
but file has 16 fields; you can supply custom names with the `names` kwarg
  "`names` kwarg" % (self.file_type, _names, self.field_count())<jupyter_text>Get genes that are likely GAF binding<jupyter_code>stat1_bed = enc.get_encode_peaks_union([
    'https://www.encodeproject.org/files/ENCFF431NLF/@@download/ENCFF431NLF.bed.gz',
    'https://www.encodeproject.org/files/ENCFF646MXG/@@download/ENCFF646MXG.bed.gz'
])
potential_gaf_genes = enc\
    .get_peak_genes_encode('https://www.encodeproject.org/files/ENCFF431NLF/@@download/ENCFF431NLF.bed.gz', 0)\
    .query('distance == 0 & ~gene.str.contains("MT").values').gene.tolist()
baseline_stat1_genes = enc\
    .get_peak_genes_encode('https://www.encodeproject.org/files/ENCFF323QQU/@@download/ENCFF323QQU.bed.gz', 0)\
    .query('distance == 0 & ~gene.str.contains("MT").values').gene.tolist()
gaf_genes =potential_gaf_genes #list(set(potential_gaf_genes) - set(baseline_stat1_genes))<jupyter_output><empty_output><jupyter_text>Read the REACTOME gene sets<jupyter_code># isgf3_genes = pd.read_csv('alpha_beta_reactome.txt', header=None)[0].tolist()
# gaf_genes = pd.read_csv('gamma_reactome.txt', header=None)[0].tolist()<jupyter_output><empty_output><jupyter_text>Read the processed RNA dataFocus on the club and bc/club cells and type I interferons for now.Encode the timestamps to integers.<jupyter_code>adata_processed = sc.read(data_path + 'HBEC_type_I_processed.h5ad')
adata = sc.read(data_path + 'HBEC_type_I_filtered_counts.h5ad')
adata.obs.donor.value_counts()
sc.pl.umap(adata_processed, color=['cell_type', 'time', 'stim'])
# adata = adata[
#     adata.obs.cell_type.isin(['club', 'basal/club']) & \
#     adata.obs.stim.isin(['control','beta'])].copy()
adata.shape<jupyter_output><empty_output><jupyter_text>Perform 1D test between control and 3 hours for gamma<jupyter_code>tp_list = ['3']
beta_dicts = {}
for tp in tp_list:
    
    adata_stim = adata[
        adata.obs.cell_type.isin(['basal/club', 'club', 'basal', 'goblet']) & \
        adata.obs.stim.isin(['control','beta']) & \
        adata.obs.time.isin(['0',tp])
        ].copy()
    time_converter={0:0, int(tp):1}
    adata_stim.obs['time_step'] = adata_stim.obs['time'].astype(int).apply(lambda x: time_converter[x])
    
    memento.create_groups(adata_stim, label_columns=['time_step', 'donor'], inplace=True, q=0.25*0.2)
    memento.compute_size_factors(adata_stim, trim_percent=0.05)
    memento.compute_1d_moments(adata_stim, inplace=True, filter_mean_thresh=0.07, min_perc_group=.9)

    memento.ht_1d_moments(
        adata_stim, 
        formula_like='1 + time_step',
        cov_column='time_step', 
        num_boot=10000, 
        verbose=1,
        num_cpus=13)
    
    beta_dicts[tp] = adata_stim.copy()

gamma_dicts = {}
for tp in tp_list:
    
    adata_stim = adata[
        adata.obs.cell_type.isin(['basal/club', 'club', 'basal', 'goblet']) & \
        adata.obs.stim.isin(['control','gamma']) & \
        adata.obs.time.isin(['0',tp])
        ].copy()
    time_converter={0:0, int(tp):1}
    adata_stim.obs['time_step'] = adata_stim.obs['time'].astype(int).apply(lambda x: time_converter[x])
    
    memento.create_groups(adata_stim, label_columns=['time_step', 'donor'], inplace=True, q=0.25*0.2)
    memento.compute_size_factors(adata_stim, trim_percent=0.05)
    memento.compute_1d_moments(adata_stim, inplace=True, filter_mean_thresh=0.07, min_perc_group=.9)

    memento.ht_1d_moments(
        adata_stim, 
        formula_like='1 + time_step',
        cov_column='time_step', 
        num_boot=10000, 
        verbose=1,
        num_cpus=13)
    
    gamma_dicts[tp] = adata_stim.copy()

alpha_dicts = {}
for tp in tp_list:
    
    adata_stim = adata[
        adata.obs.cell_type.isin(['basal/club', 'club', 'basal', 'goblet']) & \
        adata.obs.stim.isin(['control','alpha']) & \
        adata.obs.time.isin(['0',tp])
        ].copy()
    time_converter={0:0, int(tp):1}
    adata_stim.obs['time_step'] = adata_stim.obs['time'].astype(int).apply(lambda x: time_converter[x])
    
    memento.create_groups(adata_stim, label_columns=['time_step', 'donor'], inplace=True, q=0.25*0.2)
    memento.compute_size_factors(adata_stim, trim_percent=0.05)
    memento.compute_1d_moments(adata_stim, inplace=True, filter_mean_thresh=0.07, min_perc_group=.9)

    memento.ht_1d_moments(
        adata_stim, 
        formula_like='1 + time_step',
        cov_column='time_step', 
        num_boot=10000, 
        verbose=1,
        num_cpus=13)
    
    alpha_dicts[tp] = adata_stim.copy()

def classify_gene(x):
    
    isgf3 = x in isgf3_genes
    gaf = x in gaf_genes
    
    if gaf & isgf3:
        return ['isgf3', 'gaf']
    elif gaf:
        return ['gaf']
    elif isgf3:
        return ['isgf3']
    else:
        return ['rest']

beta_dfs = []
for tp in tp_list:
    
    result_1d_stim = memento.get_1d_ht_result(beta_dicts[tp])
    result_1d_stim['isgf3'] = result_1d_stim.gene.apply(lambda x: x in isgf3_genes)
    result_1d_stim['gaf'] = result_1d_stim.gene.apply(lambda x: x in gaf_genes)
    result_1d_stim['cat'] = result_1d_stim.gene.apply(classify_gene)
    result_1d_stim['de_fdr'] = memento.util._fdrcorrect(result_1d_stim['de_pval'])

    upreg_1d_stim = result_1d_stim.query('de_fdr < 0.1 & de_coef > 0 & ~gene.str.contains("MT-").values').copy()
    upreg_1d_stim['tp'] = int(tp)
    beta_dfs.append(upreg_1d_stim)
beta_df = pd.concat(beta_dfs)
beta_df = beta_df.explode('cat').sort_values('cat', ascending=False)

gamma_dfs = []
for tp in tp_list:
    
    result_1d_stim = memento.get_1d_ht_result(gamma_dicts[tp])
    result_1d_stim['isgf3'] = result_1d_stim.gene.apply(lambda x: x in isgf3_genes)
    result_1d_stim['gaf'] = result_1d_stim.gene.apply(lambda x: x in gaf_genes)
    result_1d_stim['cat'] = result_1d_stim.gene.apply(classify_gene)
    result_1d_stim['de_fdr'] = memento.util._fdrcorrect(result_1d_stim['de_pval'])

    upreg_1d_stim = result_1d_stim.query('de_fdr < 0.1 & de_coef > 0 & ~gene.str.contains("MT-").values').copy()
    upreg_1d_stim['tp'] = int(tp)
    gamma_dfs.append(upreg_1d_stim)
gamma_df = pd.concat(gamma_dfs)
gamma_df = gamma_df.explode('cat').sort_values('cat', ascending=False)

alpha_dfs = []
for tp in tp_list:
    
    result_1d_stim = memento.get_1d_ht_result(alpha_dicts[tp])
    result_1d_stim['isgf3'] = result_1d_stim.gene.apply(lambda x: x in isgf3_genes)
    result_1d_stim['gaf'] = result_1d_stim.gene.apply(lambda x: x in gaf_genes)
    result_1d_stim['cat'] = result_1d_stim.gene.apply(classify_gene)
    result_1d_stim['de_fdr'] = memento.util._fdrcorrect(result_1d_stim['de_pval'])

    upreg_1d_stim = result_1d_stim.query('de_fdr < 0.1 & de_coef > 0 & ~gene.str.contains("MT-").values').copy()
    upreg_1d_stim['tp'] = int(tp)
    alpha_dfs.append(upreg_1d_stim)
alpha_df = pd.concat(alpha_dfs)
alpha_df= alpha_df.explode('cat').sort_values('cat', ascending=False)

# plt.figure(figsize=(3, 5))
plt.figure(figsize=(5, 3))
plt.subplots_adjust(wspace=0.5)
plt.subplot(1, 2, 1)
sns.boxplot(x='tp', y='de_coef', hue='cat', data=gamma_df, fliersize=1)
plt.legend(loc='upper center', ncol=3, bbox_to_anchor=(1, 1.2))
plt.subplot(1, 2, 2)
sns.boxplot(x='tp', y='dv_coef', hue='cat', data=gamma_df, fliersize=1)
plt.legend([])
plt.ylim(-2.5, 2);
# plt.figure(figsize=(3, 5))
plt.figure(figsize=(5, 3))
plt.subplots_adjust(wspace=0.5)
plt.subplot(1, 2, 1)
sns.boxplot(x='tp', y='de_coef', hue='cat', data=beta_df, fliersize=1)
plt.legend(loc='upper center', ncol=3, bbox_to_anchor=(1, 1.2))
plt.subplot(1, 2, 2)
sns.boxplot(x='tp', y='dv_coef', hue='cat', data=beta_df, fliersize=1)
plt.legend([])
plt.ylim(-2.5, 2);
# plt.figure(figsize=(3, 5))
plt.figure(figsize=(5, 3))
plt.subplots_adjust(wspace=0.5)
plt.subplot(1, 2, 1)
sns.boxplot(x='tp', y='de_coef', hue='cat', data=alpha_df, fliersize=1)
plt.legend(loc='upper center', ncol=3, bbox_to_anchor=(1, 1.2))
plt.subplot(1, 2, 2)
sns.boxplot(x='tp', y='dv_coef', hue='cat', data=alpha_df, fliersize=1)
plt.legend([])
plt.ylim(-2.5, 2);
gamma_df['stim'] = 'gamma'
beta_df['stim'] = 'beta'
alpha_df['stim'] = 'alpha'
sns.boxplot(x='cat', y='dv_coef', hue='stim', data=pd.concat([gamma_df, alpha_df, beta_df]), fliersize=2)
plt.legend(ncol=3, loc='upper center', bbox_to_anchor=(0.5,1.2))
plt.ylim(-2.5, 2)
sns.boxplot(x='cat', y='de_coef', hue='stim', data=pd.concat([gamma_df, alpha_df, beta_df]), fliersize=2)
plt.legend(ncol=3, loc='upper center', bbox_to_anchor=(0.5,1.2))
# plt.ylim(-2.5, 2)
stats.mannwhitneyu(gamma_df.query('cat == "gaf"').dv_coef.values, gamma_df.query('cat == "rest"').dv_coef.values)<jupyter_output><empty_output><jupyter_text>Predicting the network topology<jupyter_code>df = beta_df.query('de_coef > ').copy()
df['prod'] = df['de_coef'] * df['dv_coef'].clip()
df.groupby('isgf3').de_coef.mean()
from sklearn import datasets, svm, linear_model
from sklearn.model_selection import train_test_split
import numpy as np

iris = datasets.load_iris()
X = df[['de_coef', 'prod']].values
y = df['isgf3'].values

# Limit to the two first classes, and split into training and test
X_train, X_test, y_train, y_test = train_test_split(X[y < 2], y[y < 2],
                                                    test_size=.3,
                                                    random_state=0)

# Create a simple classifier
classifier = linear_model.LogisticRegression(random_state=0)
classifier.fit(X_train, y_train)
y_score = classifier.predict_proba(X_test)[:, 1]

from sklearn.metrics import average_precision_score, roc_auc_score
average_precision = average_precision_score(y_test, y_score)
roc_score = roc_auc_score(y_test, y_score)


print('Average precision-recall score: {0:0.2f}'.format(
      average_precision))
print('ROC AUC score: {0:0.2f}'.format(
      roc_score))

from sklearn.metrics import precision_recall_curve, roc_curve

fpr, tpr, t_roc = roc_curve(y_test, y_score)
p, r, t = precision_recall_curve(y_test, y_score)

plt.plot(fpr, tpr)
plt.plot(r, p)
result_1d_stim = memento.get_1d_ht_result(adata_beta_stim)
result_1d_stim['isgf3'] = result_1d_stim.gene.apply(lambda x: x in isgf3_genes)
result_1d_stim['gaf'] = result_1d_stim.gene.apply(lambda x: x in gaf_genes)
result_1d_stim['de_fdr'] = memento.util._fdrcorrect(result_1d_stim['de_pval'])

upreg_1d_stim = result_1d_stim.query('de_fdr < 0.1 & de_coef > 0 & ~gene.str.contains("MT-").values')
# de_genes = result_1d_overall.query('de_fdr < 0.05 & de_coef > 0 & ~gene.str.contains("MT-").values').gene.tolist()
isgf3_coef_beta_stim = upreg_1d_beta_stim.query('isgf3').dv_coef
gaf_coef_beta_stim = upreg_1d_beta_stim.query('gaf').dv_coef
rest_beta_stim = upreg_1d_beta_stim.query('~gaf & ~isgf3').dv_coef
upreg_1d_beta_stim.query('gaf').shape
print(upreg_1d_beta_stim.query('isgf3').de_coef.mean())
print(upreg_1d_beta_stim.query('gaf').de_coef.mean())
print(upreg_1d_beta_stim.query('~gaf & ~isgf3').de_coef.mean())
print(upreg_1d_beta_stim.query('isgf3').dv_coef.mean())
print(upreg_1d_beta_stim.query('gaf').dv_coef.mean())
print(upreg_1d_beta_stim.query('~gaf & ~isgf3').dv_coef.mean())
upreg_1d_beta_stim.query('gaf').de_coef.mean()
result_1d_overall.query('gene == "HES4"')<jupyter_output><empty_output><jupyter_text>Plot the progression for these DE genes<jupyter_code>adata_processed.obs.cell_type.drop_duplicates()
def plot_progression(ct, gene, stim):
    
    adata = sc.read(data_path + 'HBEC_type_I_filtered_counts.h5ad')
    
    adata = adata[
        adata.obs.cell_type.isin(ct) & \
        adata.obs.stim.isin(['control',stim])].copy()
    time_converter={0:0, 3:1, 6:2, 9:3, 24:4, 48:5}
    adata.obs['time_step'] = adata.obs['time'].astype(int).apply(lambda x: time_converter[x])
    memento.create_groups(adata, label_columns=['time_step'], inplace=True, q=0.25*0.2)
    memento.compute_size_factors(adata, trim_percent=0.05)
    memento.compute_1d_moments(adata, inplace=True, filter_mean_thresh=0.2, min_perc_group=.9, filter_genes=False)
    moments_mean, moments_var, _ = memento.get_1d_moments(adata)
    plot_time_expr(gene, moments_mean, moments_var, ct, stim)
    
def plot_time_expr(gene, moments_mean, moments_var, ct, stim):
    
    plt.figure(figsize=(10,3));
    plt.subplots_adjust(wspace=0.3)
    plt.subplot(1, 2, 1);
    plt.title('{} mean'.format(gene))
    plt.scatter(np.arange(6), moments_mean[moments_mean.columns.sort_values()].query('gene == "{}"'.format(gene)).iloc[0, 1:].values)
    plt.subplot(1, 2, 2);
    plt.title('{} variability {} {}'.format(gene, ct, stim))
    plt.scatter(np.arange(6), moments_var[moments_var.columns.sort_values()].query('gene == "{}"'.format(gene)).iloc[0, 1:].values)
    plt.xlabel('tp'); plt.ylabel('variability')
for ct in adata_processed.obs.cell_type.drop_duplicates():
    print(ct)
    plot_progression([ct], 'ACE2', 'gamma')
adata = sc.read(data_path + 'HBEC_type_I_filtered_counts.h5ad')
adata.obs.cell_type.value_counts()
adata = adata[
#     adata.obs.cell_type.isin(['basal/club', 'club', 'basal', 'goblet']) & \
    adata.obs.stim.isin(['control','gamma'])].copy()
time_converter={0:0, 3:1, 6:2, 9:3, 24:4, 48:5}
adata.obs['time_step'] = adata.obs['time'].astype(int).apply(lambda x: time_converter[x])
memento.create_groups(adata, label_columns=['time_step'], inplace=True, q=0.25*0.2)
memento.compute_size_factors(adata, trim_percent=0.05)
memento.compute_1d_moments(adata, inplace=True, filter_mean_thresh=0.2, min_perc_group=.9, filter_genes=False)
moments_mean, moments_var, _ = memento.get_1d_moments(adata)
def plot_time_expr(gene):
    
    plt.figure(figsize=(10,3));
    plt.subplots_adjust(wspace=0.3)
    plt.subplot(1, 2, 1);
    plt.title('{} mean'.format(gene))
    plt.scatter(np.arange(6), moments_mean[moments_mean.columns.sort_values()].query('gene == "{}"'.format(gene)).iloc[0, 1:].values)
    plt.subplot(1, 2, 2);
    plt.title('{} variability'.format(gene))
    plt.scatter(np.arange(6), moments_var[moments_var.columns.sort_values()].query('gene == "{}"'.format(gene)).iloc[0, 1:].values)
    plt.xlabel('tp'); plt.ylabel('variability')<jupyter_output><empty_output><jupyter_text>for gene in de_genes:        plot_time_expr(gene)    plt.savefig('figures/{}.pdf'.format(gene))    plt.close()<jupyter_code>plot_time_expr('ACE2');
plot_time_expr('APOBEC3A');
plt.figure(figsize=(6, 3));
plot_time_expr('IRF1');
moments_mean[moments_mean.columns.sort_values()].query('gene == "ISG15"')
moments_var[moments_var.columns.sort_values()].query('gene == "ISG15"')
moments_mean[moments_mean.columns.sort_values()].query('gene == "SAT1"')
moments_var[moments_var.columns.sort_values()].query('gene == "SAT1"')
adata.shape
memento.ht_1d_moments(
    adata, 
    formula_like='1 + time_step',
    cov_column='time_step', 
    num_boot=10000, 
    verbose=1,
    num_cpus=13)
moments_1d = memento.get_1d_moments(adata)
moments_1d[1][moments_1d[1].columns.sort_values()].query('gene == "TXN"')
result_1d = memento.get_1d_ht_result(adata)

moments_1d[1].query('gene == "SAT1"')
result_1d.query('gene == "SAT1"')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16734939759036144, 'fraction_numerical': 0.029578313253012047, 'mean_word_length': 5.215275177836016, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 42, 'https://': 8, 'lorem ipsum': 0, 'www.': 8, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '284052', 'license': 'MIT', 'n_tokens_mistral': 7257, 'n_tokens_neox': 6717, 'n_words': 907, 'path': 'analysis/ifn_hbec/chipseq/chip_variability.ipynb', 'repo_name': 'yelabucsf/scrna-parameter-estimation', 'size': '16600'}