<jupyter_start><jupyter_text>Variability and coexpression of ISGs at baseline<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 itertools
import sys
sys.path.append('/home/ssm-user/Github/scrna-parameter-estimation/dist/memento-0.0.6-py3.8.egg')
sys.path.append('/home/ssm-user/Github/misc-seq/miscseq')
import encode
import memento
data_path = '/data_volume/memento/hbec/'
fig_path = '/data/home/Github/scrna-parameter-estimation/figures/fig6/'
import matplotlib
matplotlib.rcParams['pdf.fonttype'] = 42
matplotlib.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)<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_deep.h5ad')
adata = sc.read(data_path + 'HBEC_type_I_filtered_counts_deep.h5ad')
adata = adata[:, ~adata.var.index.str.startswith('MT-')].copy()
# adata.obs['cell_type'] = adata.obs['cell_type'].apply(lambda x: x if x != 'basal/club' else 'bc')
# adata.obs['cell_type'] = adata.obs['cell_type'].apply(lambda x: x if x != 'ionocyte/tuft' else 'ion-tuft')<jupyter_output>/home/ssm-user/anaconda3/envs/single_cell/lib/python3.8/site-packages/pandas/core/arrays/categorical.py:2487: FutureWarning: The `inplace` parameter in pandas.Categorical.remove_unused_categories is deprecated and will be removed in a future version.
  res = method(*args, **kwargs)<jupyter_text>sc.pl.umap(adata_processed, color=['cell_type', 'time', 'stim'])<jupyter_code>adata.shape
converter = {'basal/club':'BC', 'basal':'B', 'ciliated':'C', 'goblet':'G', 'ionocyte/tuft':'IT', 'neuroendo':'N'}
adata.obs['ct'] = adata.obs['cell_type'].apply(lambda x: converter[x])
adata_processed.obs['ct'] = adata_processed.obs['cell_type'].apply(lambda x: converter[x])<jupyter_output><empty_output><jupyter_text>Define ISG classes<jupyter_code>isg_classes = pd.read_csv('../isg_classes.csv')
t1_isg = isg_classes.query('overall_type =="type1"').gene.tolist()
t2_isg = isg_classes.query('overall_type =="type2"').gene.tolist()
shared_isg = isg_classes.query('overall_type == "shared"').gene.tolist()
all_isg = isg_classes[isg_classes.sum(axis=1) > 0]
cytokines = pd.read_csv('mart_export.txt')['Gene name'].tolist()
# cytokines = v_df.query('cytokine').gene.tolist()<jupyter_output><empty_output><jupyter_text>Setup memento<jupyter_code>def assign_q(batch):
    
    if batch == 0:
        return 0.387*0.25
    elif batch == 1:
        return 0.392*0.25
    elif batch == 2:
        return 0.436*0.25
    else:
        return 0.417*0.25
adata.obs['q'] = adata.obs['batch'].apply(assign_q)
memento.setup_memento(adata, q_column='q', trim_percent=0.05)<jupyter_output>Version 0.0.6<jupyter_text>Variability of DEGs in control cells<jupyter_code>m_list, v_list = [], []
for ct in ['C', 'B', 'BC']:
    ctrl = adata[(adata.obs.stim=='control') & adata.obs.ct.isin([ct])].copy().copy()
    ctrl.obs['group'] = 1
    memento.create_groups(ctrl, label_columns=['group', 'ct'])
    memento.compute_1d_moments(ctrl, min_perc_group=.9) 
    
    m_df, v_df = memento.get_1d_moments(ctrl, groupby='group')
    
    v_df['cytokine'] = v_df['gene'].apply(lambda x: 'cytokine' if x in cytokines else 'rest')
    v_df = v_df.merge(isg_classes[['gene', 'overall_type']], on='gene', how='left').fillna('rest')
    m_df['ct'] = ct
    v_df['ct'] = ct
    
    m_list.append(m_df.copy())
    v_list.append(v_df.copy())
m_df = pd.concat(m_list)
v_df = pd.concat(v_list)
v_df['is_isg'] = v_df['overall_type'].apply(lambda x: 'ISG' if x != 'rest' else 'rest')
plt.scatter(m_df['group_1'], v_df['group_1'], s=1, color='gray')
# plt.plot([-3, 5], [0, 0])
# plt.scatter(m_df.query('gene in @t1_isg')['group_1'], v_df.query('gene in @t1_isg')['group_1'], s=5, linewidth=1, color='k')

overlap = v_df.query('is_isg !="rest" & cytokine != "rest"').copy()
overlap['gene_class'] = 'cytokine'
def classify(row):
    
    if row['is_isg'] == 'ISG':
        return 'ISG'
    elif row['cytokine'] == 'cytokine':
        return 'cytokine'
    else:
        return 'rest'
v_df['gene_class'] = v_df.apply(classify, axis=1)
barplot_df = pd.concat([v_df, overlap])
plt.figure(figsize=(3,2))
sns.boxplot(x='ct', y='group_1', hue='gene_class', data=v_df, fliersize=0.0, palette='Set2')
# sns.stripplot(x='ct', y='group_1', hue='overall_type', data=v_df)

plt.legend(bbox_to_anchor=(0.6,1.2), frameon=False, title=None, ncol=3, loc='upper center', handletextpad=0.5, columnspacing=1); 
plt.ylabel('Variability'); plt.xlabel(None)
plt.ylim(-2,5)
plt.savefig('isg_cytokine_var_baseline.pdf', bbox_inches='tight')
plt.figure(figsize=(2,2))
sns.boxplot(x='ct', y='group_1', hue='cytokine', data=v_df, fliersize=0.0, palette='Set2')
# sns.stripplot(x='ct', y='group_1', hue='overall_type', data=v_df)

plt.legend(bbox_to_anchor=(0.6,1.2), frameon=False, title=None, ncol=2, loc='upper center', handletextpad=0.5, columnspacing=1); 
plt.ylabel('Variability'); plt.xlabel(None)
plt.ylim(-2,4)
plt.savefig('cytokine_var_baseline.pdf', bbox_inches='tight')
plt.figure(figsize=(2,2))
sns.boxplot(x='ct', y='group_1', hue='is_isg', data=v_df, fliersize=0.0, palette='Set2')
# sns.stripplot(x='ct', y='group_1', hue='overall_type', data=v_df)

plt.legend(bbox_to_anchor=(0.5,1.2), frameon=False, title=None, ncol=2, loc='upper center', handletextpad=0.5, columnspacing=1);  
plt.ylabel('Variability'); plt.xlabel(None)
plt.ylim(-2,3)
plt.savefig('isg_var_baseline.pdf', bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>Find a UMAP example<jupyter_code>a = v_df.merge(m_df[['gene', 'ct', 'group_1']], on=['gene', 'ct'])
a.head(2)
adata_processed.obs['is_control'] = (adata_processed.obs['stim']=='control').astype(int)
sc.pl.umap(adata_processed, color='is_control', s=20)
adata_ctrl = adata_processed[adata_processed.obs.stim=='control'].copy()
sc.tl.umap(adata_ctrl)<jupyter_output>/home/ssm-user/anaconda3/envs/single_cell/lib/python3.8/site-packages/pandas/core/arrays/categorical.py:2487: FutureWarning: The `inplace` parameter in pandas.Categorical.remove_unused_categories is deprecated and will be removed in a future version.
  res = method(*args, **kwargs)
/home/ssm-user/anaconda3/envs/single_cell/lib/python3.8/site-packages/sklearn/manifold/_spectral_embedding.py:245: UserWarning: Graph is not fully connected, spectral embedding may not work as expected.
  warnings.warn("Graph is not fully connected, spectral embedding"<jupyter_text>Heatmap for correlations at baseline<jupyter_code>ct = 'C'
ctrl = adata[(adata.obs.stim=='control') & adata.obs.ct.isin([ct])].copy().copy()
ctrl.obs['group'] = 1
memento.create_groups(ctrl, label_columns=['group', 'ct'])
memento.compute_1d_moments(ctrl, min_perc_group=.9) 
included_de_genes = list(set(t1_isg+t2_isg+shared_isg) & set(ctrl.var.index))
memento.compute_2d_moments(ctrl, gene_pairs=list(itertools.combinations(included_de_genes,2)))
c_df = memento.get_2d_moments(ctrl, groupby='ct')
ct_c_df = c_df[['gene_1', 'gene_2', 'ct_C']]
ct_heatmap = pd.DataFrame(np.zeros((len(included_de_genes), len(included_de_genes))), columns=included_de_genes, index=included_de_genes)
for idx, row in ct_c_df.iterrows():
    
    ct_heatmap.loc[row['gene_1'], row['gene_2']] = row['ct_C']
    ct_heatmap.loc[row['gene_2'], row['gene_1']] = row['ct_C']
print(ct_heatmap.shape)
cg = sns.clustermap(ct_heatmap.fillna(0.0), center=0, vmax=0.3, metric=''); plt.close()
row_inds = cg.dendrogram_row.reordered_ind
col_inds = cg.dendrogram_col.reordered_ind
print(ct_heatmap.shape)

def classify_color(x):
    
    if x in shared_isg:
        return -1
    if x in t1_isg:
        return 0
    if x in t2_isg:
        return 1
    return 'impossible'
# class_list = np.array([classify(g) for g in included_de_genes])[:, np.newaxis]
from matplotlib.colors import LinearSegmentedColormap


myColors = sns.color_palette('husl')[:3]
class_list = np.array([classify_color(g) for g in included_de_genes])[np.newaxis, :]
cmap = LinearSegmentedColormap.from_list('Custom', myColors, len(myColors))
'ISG15' in included_de_genes
ct_heatmap.index.get_loc(g)
row_inds.indexof(ct_heatmap.index.get_loc(g))
g = 'ISG15'
np.where(np.array(row_inds)==ct_heatmap.index.get_loc(g))[0][0]
for i in selected_genes[5]:
    print(i)
colors = sns.color_palette('husl')
fig = plt.figure(figsize=(20,20))
ax1 = plt.subplot2grid((21,20), (1,0), colspan=19, rowspan=19)
ax2 = plt.subplot2grid((21,20), (0,0), colspan=19, rowspan=1)

sns.heatmap(ct_heatmap.iloc[row_inds, col_inds], cmap='viridis', vmax=0.8, vmin=0.1, ax=ax1, yticklabels=False, xticklabels=False, cbar_kws={"shrink": 0.5})
sns.heatmap(class_list[:, row_inds], cmap=cmap, ax=ax2, yticklabels=False, xticklabels=False,)
colorbar = ax2.collections[0].colorbar
colorbar.set_ticks([-0.66, 0,0.66,])
colorbar.set_ticklabels(['shared', 'type1', 'type2'])

# some_genes = list(set(['ISG15', 'MX1', 'IFIT3', 'OAS1', 'OAS2', 'OASL', 'MX2', 'IFI6', 'STAT1'])&set(included_de_genes))
# for g in some_genes:
#     idx = np.where(np.array(row_inds)==ct_heatmap.index.get_loc(g))[0][0]
#     ax1.plot([0, len(included_de_genes)], [idx, idx], color='magenta')
#     ax1.plot([idx, idx],[0, len(included_de_genes)], color='magenta',)
    
selected_genes = {}
locations = [
    (235, 245),
    (215, 224),
    (174, 183),
    (80, 90),
    (18, 29),
    (322, 350),
    (100, 100)
]
colors = sns.color_palette('husl', len(locations))
i = 0

for x1, x2 in locations:
    
    ax1.plot([x1, x1],[0, len(included_de_genes)], color=colors[i],)
    ax1.plot([x2, x2],[0, len(included_de_genes)], color=colors[i])

    selected_genes[i]= ct_heatmap.index[np.array(row_inds)[x1:x2]]
    
    i+=1
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)
fig = plt.figure(figsize=(4,3.5))
ax1 = plt.subplot2grid((21,20), (1,0), colspan=19, rowspan=19)
ax2 = plt.subplot2grid((21,20), (0,0), colspan=19, rowspan=1)

sns.heatmap(ct_heatmap.iloc[row_inds[100:355], col_inds[100:355]], cmap='viridis', vmax=0.8, vmin=0.1, ax=ax1, yticklabels=False, xticklabels=False, cbar_kws={"shrink": 0.5})
sns.heatmap(class_list[:, row_inds[100:355]], cmap=cmap, ax=ax2, yticklabels=False, xticklabels=False, cbar_kws={"shrink": 3})
colorbar = ax2.collections[0].colorbar
colorbar.set_ticks([-0.66, 0,0.66,])
colorbar.set_ticklabels(['shared', 'type1', 'type2'])
plt.savefig('baseline_isg_corr.png', dpi=800, bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>Compare correlation to random genes<jupyter_code>ct = 'C'
ctrl = adata[(adata.obs.stim=='control') & adata.obs.ct.isin([ct])].copy().copy()
ctrl.obs['group'] = 1
memento.create_groups(ctrl, label_columns=['group', 'ct'])
memento.compute_1d_moments(ctrl, min_perc_group=.9) 
included_de_genes = list(set(t1_isg+t2_isg+shared_isg) & set(ctrl.var.index))
memento.compute_2d_moments(ctrl, gene_pairs=list(itertools.combinations(included_de_genes,2)))
c_df = memento.get_2d_moments(ctrl, groupby='ct')
ctrl.var.index.shape
list(zip(np.random.choice(ctrl.var.index, 1000), np.random.choice(ctrl.var.indexm, 1000)))
# memento.compute_2d_moments(ctrl, gene_pairs=list(itertools.combinations(np.random.choice(ctrl.var.index, size=len(included_de_genes), replace=False),2)))
memento.compute_2d_moments(ctrl, gene_pairs=list(zip(np.random.choice(ctrl.var.index, 1000), np.random.choice(ctrl.var.index, 1000))))

random_c_df = memento.get_2d_moments(ctrl, groupby='ct')
c_df['which'] = 'between ISG'
random_c_df['which'] = 'random'
comb = pd.concat([c_df, random_c_df])
plt.figure(figsize=(2,2))

sns.kdeplot(c_df['ct_C'])
sns.kdeplot(random_c_df['ct_C'])
plt.xlabel('correlation')
plt.ylabel(None); 
plt.legend(['between ISGs', 'random'], frameon=False, loc='upper center', bbox_to_anchor=(0.5,1.2), ncol=2, columnspacing=0.7)<jupyter_output><empty_output><jupyter_text>Correlation to TFs<jupyter_code>tf_df = pd.read_csv('human_tf.txt', sep='\t')
tf_list = tf_df['Symbol'].tolist()
tf_list += ['CIITA']
adata_stim
adata
stim = 'control'
ct = 'C'
adata_stim = adata.copy()[

    adata.obs.ct.isin([ct]) & \
    adata.obs.stim.isin([stim])].copy()
adata_stim.obs['group'] = 1

memento.create_groups(adata_stim, label_columns=['donor', 'group'])
memento.compute_1d_moments(adata_stim, min_perc_group=.9)

available_tfs = list(set(tf_list) & set(adata_stim.var.index))
available_isgs= list(set(t1_isg+t2_isg+shared_isg) & set(adata_stim.var.index))

# moment_1d_df = memento.get_1d_moments(adata_stim, groupby='is_stim')

memento.compute_2d_moments(adata_stim, list(itertools.product(available_tfs, available_isgs)))
corr_df = memento.get_2d_moments(adata_stim, groupby='group')
avg_corr = corr_df.groupby('gene_1')['group_1'].mean().sort_values(ascending=False).reset_index(name='corr')
avg_corr = corr_df.groupby('gene_1')['group_1'].mean().sort_values(ascending=False).reset_index(name='corr')
for idx, row in avg_corr.iterrows():
    plt.plot([idx, idx], [0, row['corr']], color='k', lw=1)
# sns.barplot(x='gene_1', y='corr', data=corr_df.groupby('gene_1')['group_1'].mean().reset_index(name='corr'));
plt.xticks([], []); plt.xlabel('Transcriptional regulator'); plt.ylabel('Avg. corr with ISGs at baseline')
plt.figure(figsize=(2,0.8))

avg_corr = corr_df.groupby('gene_1')['group_1'].mean().sort_values(ascending=False).reset_index(name='corr').head(12)
sns.barplot(x='gene_1', y='corr', data=avg_corr);
plt.xticks(rotation=90)
plt.xlabel(None); plt.ylabel(None)
plt.savefig('baseline_tf_isg_corr_zoom.pdf', bbox_inches='tight')
avg_corr = corr_df.groupby('gene_1')['group_1'].mean().sort_values(ascending=False).reset_index(name='corr')

plt.figure(figsize=(4,2))
for idx, row in avg_corr.iterrows():
    plt.plot([idx, idx], [0, row['corr']], color='k', lw=0.5)
# sns.barplot(x='gene_1', y='corr', data=corr_df.groupby('gene_1')['group_1'].mean().reset_index(name='corr'));
plt.xticks([], []); plt.xlabel('Transcriptional regulator'); plt.ylabel('Avg. corr with ISGs')
plt.savefig('baseline_tf_isg_corr.png', dpi=800, bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>STAT1 binding site enrichment<jupyter_code>enc = encode.Encode('/home/ssm-user/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()
# Having a binding site makes you more variable at baseline
v_df['tss_binding'] = v_df['gene'].isin(isgf3_genes)
v_df['deg_type'] = v_df['class'].astype(str).apply(lambda x: {
    'uniq_gamma':'uniq\ntype II', 'uniq_alpha/beta':'uniq\ntype I', 'uniq_lambda':'uniq\ntype III', 'partially_shared':'partial', 'shared':'shared', 'none-DE':'non-DE'}[x])
sns.barplot(x='deg_type', y='tss_binding', data=v_df[v_df['deg_type']!='partial'].sort_values('deg_type'), capsize=0.2, palette='Set2');
plt.ylabel('Proportion with ISGF3 \nbinding site in TSS')
# plt.xticks([0, 1, 2, 3, 4], ['none\nDE', 'shared','type III', 'type II', 'type I']);
v_df.head(2)
stats.mannwhitneyu(v_df.query('tss_binding').ct_C, v_df.query('~tss_binding').ct_C)
plt.figure(figsize=(2,4))
sns.boxplot(x='tss_binding', y='ct_BC', data=v_df, fliersize=0.1); plt.ylabel('Variability of genes')
plt.ylim(-2, 2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1863357843137255, 'fraction_numerical': 0.029963235294117648, 'mean_word_length': 6.6195144724556485, 'pattern_counts': {'":': 3, '<': 42, '<?xml version=': 0, '>': 43, 'https://': 4, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '77751', 'license': 'MIT', 'n_tokens_mistral': 7091, 'n_tokens_neox': 6624, 'n_words': 950, 'path': 'analysis/ifn_hbec/version3/baseline/deg_in_baseline.ipynb', 'repo_name': 'yelabucsf/scrna-parameter-estimation', 'size': '16320'}