<jupyter_start><jupyter_text>Coexpresion 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.4-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>Read and select TFs<jupyter_code>tf_df = pd.read_csv('DatabaseExtract_v_1.01.csv', index_col=0)
tf_df = tf_df[tf_df['TF assessment'] == 'Known motif']
tfs = tf_df['HGNC symbol'].tolist()
condensed_tfs = []
for tf in tfs:
    if 'STAT' in tf:
        condensed_tfs.append(tf)
    elif 'IRF' in tf:
        condensed_tfs.append(tf)
    elif 'NFKB' in tf:
        condensed_tfs.append(tf)
    elif 'JUN' in tf:
        condensed_tfs.append(tf)
    elif 'MYC' in tf:
        condensed_tfs.append(tf)
    elif 'FOS' in tf:
        condensed_tfs.append(tf)
    else:
        continue<jupyter_output><empty_output><jupyter_text>Read 2d HT result comparing to control<jupyter_code>cts = [['ciliated'], ['bc', 'club'], ['basal']]
tps = ['3', '6', '9', '24', '48']
stims = ['alpha', 'beta', 'gamma', 'lambda']
def read_result(tf, ct):
    
    ht_dict = {}
    moments_dict = {}
    
    for stim in stims:
        ht_dict[stim] = {}
        moments_dict[stim] = {}
        for tp in tps:
            try:
                adata = sc.read('/data_volume/ifn_hbec/coex_test/{}_{}_{}_{}.h5ad'.format('-'.join(ct), stim, tp, tf))
            except:
                continue
            ht_dict[stim][tp] = memento.get_2d_ht_result(adata)
            ht_dict[stim][tp]['corr_fdr'] = memento.util._fdrcorrect(ht_dict[stim][tp]['corr_pval'])
            moments_dict[stim][tp] = memento.get_2d_moments(adata, groupby='time_step')
    return ht_dict,moments_dict
all_result = {}

for tf in condensed_tfs:
    all_result[tf] = {}
    for ct in cts:

        key = '-'.join(ct)
        all_result[tf][key] = {}
        all_result[tf][key]['ht'], all_result[tf][key]['moments'] = read_result(tf, ct)<jupyter_output>/data/home/Github/scrna-parameter-estimation/dist/memento-0.0.4-py3.7.egg/memento/memento.py:539: RuntimeWarning: invalid value encountered in true_divide<jupyter_text>Number of DC genes over time<jupyter_code>all_result['STAT1']['bc-club']['moments']['alpha']['6']#.sort_values('corr_fdr')#.head(100)#['corr_pval'].hist()
def plot_num_de_dv(ct, tf):
    
    dc_count = np.zeros((len(stims), len(tps)))
    
    for i, stim in enumerate(stims):
        for j, tp in enumerate(tps):
            
            if ct == 'all':
                
                combined_df = pd.concat([all_result[tf]['-'.join(c)]['ht'][stim][tp] for c in cts])
                dc_count[i, j] = combined_df.query('corr_fdr < 0.1'.format(tf)).gene_2.drop_duplicates().shape[0]
            else:
                df = all_result[tf][ct]['ht'][stim][tp]
                dc_count[i, j] = df.query('corr_fdr < 0.1').shape[0]
                
    dc_count = np.hstack([np.zeros((len(stims), 1)), dc_count])
    plt.figure(figsize=(5, 2));
    for i in range(len(stims)):
        plt.plot([0.0]+[float(int(tp)) for tp in tps], dc_count[i, :], '--o')
    plt.ylabel('DC count'); plt.legend(stims, loc='upper center', ncol=4, bbox_to_anchor=(1.3, 1.25), frameon=False); plt.xlabel('time')
    
plot_num_de_dv('bc-club', 'STAT1')
# plot_num_de_dv('ciliated', 'STAT1')
def plot_num_de_dv(ct, tf):
    
    dc_count = np.zeros((len(stims), len(tps)))
    
    for i, stim in enumerate(stims):
        for j, tp in enumerate(tps):
            
            if ct == 'all':
                
                df_list = []
                for c in cts:
                    df = all_result['-'.join(c)]['ht'][stim][tp].query('gene_1 == "{}"'.format(tf)).copy()
                    df['corr_fdr'] = memento.util._fdrcorrect(df['corr_pval'])
                    df_list.append(df)
                
                combined_df = pd.concat(df_list)
                dc_count[i, j] = combined_df.query('corr_pval < 0.1'.format(tf)).gene_2.drop_duplicates().shape[0]
            else:
                df = all_result[ct]['ht'][stim][tp].query('gene_1 == "{}"'.format(tf)).copy()
                df['corr_fdr'] = memento.util._fdrcorrect(df['corr_pval'])
                dc_count[i, j] = df.query('corr_pval < 0.1').shape[0]
                
    dc_count = np.hstack([np.zeros((len(stims), 1)), dc_count])
    plt.figure(figsize=(5, 2));
    for i in range(len(stims)):
        plt.plot([0.0]+[float(int(tp)) for tp in tps], dc_count[i, :], '--o')
    plt.ylabel('DC count'); plt.legend(stims, loc='upper center', ncol=4, bbox_to_anchor=(1.3, 1.25), frameon=False); plt.xlabel('time')
    
plot_num_de_dv('bc-club', 'STAT1')
plot_num_de_dv('ciliated', 'STAT1')<jupyter_output><empty_output><jupyter_text>Get HT result and moments<jupyter_code>moment_df_list = []
ht_df_list = []
for ct in cts:
    
    key = key = '-'.join(ct)
    for stim in stims:
        for tp in tps:
            
            df = all_result[key]['ht'][stim][tp]
            df['stim'] = stim
            df['ct'] = key
            df['tp'] = int(tp)
            
            ht_df_list.append(df.copy())
ht_df = pd.concat(ht_df_list)
for ct in cts:
    
    key = key = '-'.join(ct)
    for stim in stims:
        for tp in tps:
            
            df_ctrl = pd.DataFrame()
            df_ctrl['gene'] = all_result[key]['moments'][stim][tp][0]['gene'].tolist()
            df_ctrl['mean'] = all_result[key]['moments'][stim][tp][0]['0'].tolist()
            df_ctrl['var'] = all_result[key]['moments'][stim][tp][1]['0'].tolist()
            df_ctrl['stim'] = 'ctrl'
            df_ctrl['ct'] = key
            df_ctrl['tp'] = 0
            
            df_stim = pd.DataFrame()
            df_stim['gene'] = all_result[key]['moments'][stim][tp][0]['gene'].tolist()
            df_stim['mean'] = all_result[key]['moments'][stim][tp][0]['1'].tolist()
            df_stim['var'] = all_result[key]['moments'][stim][tp][1]['1'].tolist()
            df_stim['stim'] = stim
            df_stim['ct'] = key
            df_stim['tp'] = int(tp)
            
            moment_df_list.append(df_ctrl)
            moment_df_list.append(df_stim)
moment_df = pd.concat(moment_df_list).drop_duplicates()<jupyter_output><empty_output><jupyter_text>Plot effect sizes<jupyter_code>def plot_es(ct):
    plt.figure(figsize=(5, 2));
    plt.subplots_adjust(wspace=0.5)
    plt.subplot(1, 2, 1)
    sns.pointplot(x='tp', y='de_coef', hue='stim', data=ht_df.query('ct == "{}"'.format(ct)), ci=66, dodge=True, linestyles='--', scale=.5); plt.legend(loc='upper center', bbox_to_anchor=(0.9, 1.25), ncol=4, frameon=False)
    plt.subplot(1, 2, 2)
    sns.pointplot(x='tp', y='dv_coef', hue='stim', data=ht_df.query('ct == "{}"'.format(ct)), ci=66, dodge=True, linestyles='--', scale=.5); plt.legend([], frameon=False)
plot_es('basal')
plot_es('bc-club')
plot_es('ciliated')<jupyter_output><empty_output><jupyter_text>Plot moments<jupyter_code>def plot_moments(ct):
    plt.figure(figsize=(5, 2));
    plt.subplots_adjust(wspace=0.5)
    plt.subplot(1, 2, 1)
    sns.pointplot(x='tp', y='mean', hue='stim', data=moment_df.query('ct == "{}"'.format(ct)), ci=66, dodge=True, linestyles='--', scale=.7); plt.legend(loc='upper center', bbox_to_anchor=(0.9, 1.25), ncol=4, frameon=False)
    plt.subplot(1, 2, 2)
    sns.pointplot(x='tp', y='var', hue='stim', data=moment_df.query('ct == "{}"'.format(ct)), ci=66, dodge=True, linestyles='--', scale=.7); plt.legend([], frameon=False)
plot_moments('basal')
plot_moments('bc-club')
plot_moments('ciliated')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18411322198138466, 'fraction_numerical': 0.016575290067576184, 'mean_word_length': 3.0287621982537236, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 28, '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': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '140076', 'license': 'MIT', 'n_tokens_mistral': 3109, 'n_tokens_neox': 2959, 'n_words': 483, 'path': 'analysis/ifn_hbec/version1/all_stim_2d_analyze.ipynb', 'repo_name': 'yelabucsf/scrna-parameter-estimation', 'size': '7843'}