<jupyter_start><jupyter_text>MHC II Analysis - PBMCs<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
import itertools
%matplotlib inline
import matplotlib
matplotlib.rcParams['pdf.fonttype'] = 42
matplotlib.rcParams['ps.fonttype'] = 42
import sys
sys.path.append('/home/ssm-user/Github/scrna-parameter-estimation/dist/memento-0.0.5-py3.8.egg')
sys.path.append('/home/ssm-user/Github/misc-seq/miscseq/')
import encode
import memento
data_path = '/data_volume/parameter_estimation/'
fig_path = '/home/ssm-user/Github/scrna-parameter-estimation/figures/fig5/supp/'<jupyter_output><empty_output><jupyter_text>Read PBMC data<jupyter_code>adata = sc.read(data_path + 'interferon_filtered.h5ad')
adata.obs['q'] = 0.07
memento.setup_memento(adata, q_column='q', trim_percent=0.15)<jupyter_output><empty_output><jupyter_text>Check correlations in ctrl<jupyter_code>adata_cm = adata.copy()[(adata.obs['stim']=='ctrl')].copy()
memento.create_groups(adata_cm, label_columns=['cell', 'ind'])
memento.compute_1d_moments(adata_cm, min_perc_group=.1)

# adata_ctrl = adata.copy()[(adata.obs['cell']=='B cells') & (adata.obs['stim']=='ctrl')].copy()
canonical_isgs = pd.read_csv('alpha_beta_reactome.txt', header=None)[0].tolist()
m, v, count = memento.get_1d_moments(adata_cm)
mhc_genes = m.query('gene.str.contains("HLA")', engine='python').gene.tolist()
tf_list = ['STAT1']

gene_list = adata_cm.var.index.tolist()
isgs = list(set(canonical_isgs) & set(gene_list))
target_list = list(set(isgs) | set(mhc_genes))

memento.compute_2d_moments(adata_cm, list(itertools.product(tf_list,target_list)))
corr_df = memento.get_2d_moments(adata_cm, groupby='cell')
corr_df.query('gene_2.str.contains("HLA")',engine='python')<jupyter_output><empty_output><jupyter_text>Check correlations in stim<jupyter_code>adata_cm = adata.copy()[(adata.obs['stim']=='stim') & (adata.obs['cell'] != 'Megakaryocytes')].copy()
memento.create_groups(adata_cm, label_columns=['cell', 'ind'])
memento.compute_1d_moments(adata_cm, min_perc_group=.7)

# adata_ctrl = adata.copy()[(adata.obs['cell']=='B cells') & (adata.obs['stim']=='ctrl')].copy()
canonical_isgs = pd.read_csv('alpha_beta_reactome.txt', header=None)[0].tolist()
m, v, count = memento.get_1d_moments(adata_cm)
mhc_genes = m.query('gene.str.contains("HLA")', engine='python').gene.tolist()
tf_list = ['STAT1']

gene_list = adata_cm.var.index.tolist()
isgs = list(set(canonical_isgs) & set(gene_list))
target_list = list(set(isgs) | set(mhc_genes))

memento.compute_2d_moments(adata_cm, list(itertools.product(tf_list,target_list)))
corr_df = memento.get_2d_moments(adata_cm, groupby='cell')
m,v = memento.get_1d_moments(adata_cm, groupby='cell')
m.query('gene.str.contains("HLA")',engine='python')
corr_df.query('gene_2.str.contains("HLA")',engine='python')<jupyter_output><empty_output><jupyter_text>Make boxplots<jupyter_code>corr_df, counts = memento.get_2d_moments(adata_cm)

def classify(g):
    
    if g in ['HLA-A', 'HLA-B', 'HLA-C']:
        return 'type1'
    elif len(g) > 5 and 'HLA' in g:
        return 'type2'
    elif 'HLA' in g:
        return 'Misc HLA'
    else:
        return 'ISG'
order = dict(zip(['cM', 'ncM', 'DC', 'B','Tc', 'Th', 'NK','Mega'], range(8)))
cts = ['CD4 T cells',  'CD14+ Monocytes', 'FCGR3A+ Monocytes', 'NK cells','CD8 T cells', 'B cells', 'Dendritic cells', 'Megakaryocytes']
label_converter = dict(zip(cts, ['Th', 'cM', 'ncM', 'NK', 'Tc', 'B', 'DC', 'Mega']))

long_df = []
for col in corr_df.columns[2:]:
    
    df = corr_df[['gene_1', 'gene_2']].copy()
    df['corr'] = corr_df[col]
    df['ct'] = label_converter[col.split('^')[1]]
    df['ind'] = col.split('^')[2]
    long_df.append(df)
    
long_df = pd.concat(long_df)

long_df =long_df.sort_values('ct', key=lambda x: x.map(order))
long_df['type'] = [classify(g) for g in long_df.gene_2]
long_df['gene'] = long_df['gene_2']
plt.figure(figsize=(10, 3))
plt.subplot(1, 2, 1)
dat = long_df.query('type=="type1"', engine='python')
plt.plot([-0.5, 7.5], [0, 0], '--', lw=3, c='k')
sns.boxplot(x='ct', y='corr', hue='gene',data=dat)
plt.subplot(1, 2, 2)
dat = long_df.query('type=="type2"', engine='python')
plt.plot([-0.5, 7.5], [0, 0], '--', lw=3, c='k')
sns.boxplot(x='ct', y='corr', hue='gene',data=dat)

plt.savefig(fig_path + 'mhc_mean_corr_pbmc.pdf', bbox_inches='tight')
dat = long_df_mhc.query('type=="type2"', engine='python')
plt.plot([-0.5, 7.5], [0, 0], '--', lw=3, c='k')
sns.boxplot(x='ct', y='corr', hue='gene_2',data=dat)
# sns.stripplot(x='ct', y='corr', data=dat)

dat = long_df_mhc.query('type=="type2"', engine='python')
plt.plot([-0.5, 7.5], [0, 0], '--', lw=3, c='k')
sns.boxplot(x='ct', y='corr',data=dat)
# sns.stripplot(x='ct', y='corr', data=dat)

sns.boxplot(x='ct', y='corr', data=long_df.query('gene_2 == "HLA-A"').sort_values('ct', key=lambda x: x.map(order)))
long_df.head(10)<jupyter_output><empty_output><jupyter_text>Get difference in STAT1 correlations<jupyter_code>adata_cm = adata.copy()[(adata.obs['cell']=='CD14+ Monocytes')].copy()
memento.create_groups(adata_cm, label_columns=['stim'])
memento.compute_1d_moments(adata_cm, min_perc_group=.7)
m, v, count = memento.get_1d_moments(adata_cm)
canonical_isgs = pd.read_csv('alpha_beta_reactome.txt', header=None)[0].tolist()
m, v, count = memento.get_1d_moments(adata_cm)
mhc_genes = m.query('gene.str.contains("HLA")', engine='python').gene.tolist()
tf_list = ['STAT1']

gene_list = adata_cm.var.index.tolist()
isgs = list(set(canonical_isgs) & set(gene_list))
target_list = list(set(isgs) | set(mhc_genes))

memento.compute_2d_moments(adata_cm, list(itertools.product(tf_list,target_list)))
memento.ht_2d_moments( 
    adata_cm, 
    formula_like='1 + stim',
    cov_column='stim', 
    num_boot=20000, 
    verbose=1,
    num_cpus=13)
corr_df = memento.get_2d_moments(adata_cm, groupby='stim')
a = memento.get_2d_ht_result(adata_cm)
a
cts = [['ciliated']]
tps = ['3', '6', '9', '24', '48']
stims = ['beta']


for ct in cts:
    for tp in tps:
        for stim in stims:
                            
            print('starting', ct, tp, stim)

            adata_stim = adata.copy()[
                adata.obs.cell_type.isin(ct) & \
                adata.obs.stim.isin(['control', stim]) & \
                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'])
            memento.compute_1d_moments(adata_stim, min_perc_group=.9)
            
            tf_list = ['STAT1']
            gene_list = adata_stim.var.index.tolist()
            
            isgs = list(set(canonical_isgs) & set(gene_list))
            
            target_list = list(set(isgs) | set(mhc_genes))
            
            memento.compute_2d_moments(adata_stim, list(itertools.product(tf_list,target_list)))

            moment_df = memento.get_2d_moments(adata_stim, groupby='time_step')
            moment_df.to_csv('/data_volume/ifn_hbec/tf_test/moments_{}_{}_{}_mhc1_20210321.csv'.format('-'.join(ct), stim, tp), index=False)
stim = 'beta'
ct = ['ciliated']
subset = adata[adata.obs['cell_type'].isin(ct) & adata.obs.stim.isin(['control', stim])].copy().copy()
memento.create_groups(subset, label_columns=['time', 'donor'])
memento.compute_1d_moments(subset, min_perc_group=.9)
mean_df = memento.get_1d_moments(subset, groupby='time')[0][['gene'] + ['time_' + str(i) for i in [0, 3, 6, 9, 24, 48]]]
var_df = memento.get_1d_moments(subset, groupby='time')[1][['gene'] + ['time_' + str(i) for i in [0, 3, 6, 9, 24, 48]]]
d = mean_df.query('gene.isin(["REL", "NFKB2"])', engine='python')
plt.figure(figsize=(5,2))
plt.subplots_adjust(wspace=0.4)
locs = [1, 2]
cs = [0, 2]
num=0
for idx, row in d.iterrows():
    plt.subplot(1,2,locs[num])
    plt.plot(row.values[1:], '-o')
    plt.title(row['gene'])
    plt.xticks([0, 1, 2, 3, 4, 5], ['0', '3', '6', '9', '24', '48'])
    if num == 0:
        plt.ylabel('Mean expr')
    plt.xlabel('time (hr)')
    num +=1

plt.savefig(fig_path + 'mhc_aux_mean.pdf', bbox_inches='tight')

d.values[2:].T
d
plt.figure(figsize=(2,2))
d = mean_df.query('gene.isin(["HLA-A", "HLA-B", "HLA-C", "B2M"])', engine='python').values[:,1:]
d = d - d[:, 0].reshape(-1,1)
plt.plot(d.T, '-o'); plt.ylabel('Relative mean expr')
plt.xticks([0, 1, 2, 3, 4, 5], ['0', '3', '6', '9', '24', '48']); plt.xlabel('time (hr)')<jupyter_output><empty_output><jupyter_text>Correlation of MHC I transcripts to STAT1<jupyter_code>canonical_isgs = pd.read_csv('alpha_beta_reactome.txt', header=None)[0].tolist()
mhc_genes = ["HLA-A", "HLA-B", "HLA-C", "B2M"]
cts = [['ciliated']]
tps = ['3', '6', '9', '24', '48']
stims = ['beta']


for ct in cts:
    for tp in tps:
        for stim in stims:
                            
            print('starting', ct, tp, stim)

            adata_stim = adata.copy()[
                adata.obs.cell_type.isin(ct) & \
                adata.obs.stim.isin(['control', stim]) & \
                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'])
            memento.compute_1d_moments(adata_stim, min_perc_group=.9)
            
            tf_list = ['STAT1']
            gene_list = adata_stim.var.index.tolist()
            
            isgs = list(set(canonical_isgs) & set(gene_list))
            
            target_list = list(set(isgs) | set(mhc_genes))
            
            memento.compute_2d_moments(adata_stim, list(itertools.product(tf_list,target_list)))

            moment_df = memento.get_2d_moments(adata_stim, groupby='time_step')
            moment_df.to_csv('/data_volume/ifn_hbec/tf_test/moments_{}_{}_{}_mhc1_20210321.csv'.format('-'.join(ct), stim, tp), index=False)
stim = 'beta'
ct = ['ciliated']
for tp in ['3', '6', '9', '24', '48']:
    
    moments = pd.read_csv('/data_volume/ifn_hbec/tf_test/moments_{}_{}_{}_mhc1_20210321.csv'.format('-'.join(ct), stim, tp)).rename(columns={'time_step_0':'time_0', 'time_step_1':'time_'+tp})
    
    if tp == '3':
        moment_df = moments[['gene_1', 'gene_2', 'time_0', 'time_3']].copy()
        
    else:
        moment_df = moment_df.merge(moments[['gene_1', 'gene_2', 'time_{}'.format(tp)]].copy(), on=['gene_1', 'gene_2'], how='inner')
plt.rcParams['axes.prop_cycle'].by_key()['color']
plt.rcParams['axes.prop_cycle'].by_key()['color'][1:5]
tf = 'STAT1'
plt.figure(figsize=(6,2))
plt.subplots_adjust(wspace=0.5)
plt.subplot(1, 3, 1);

d = mean_df.query('gene.isin(["HLA-A", "HLA-B", "HLA-C", "B2M"])', engine='python').values[:,1:]
d = d - d[:, 0].reshape(-1,1)
plt.plot(d.T, '-o'); plt.ylabel('Relative mean expr')
plt.xticks([0, 1, 2, 3, 4, 5], ['0', '3', '6', '9', '24', '48']); plt.xlabel('time (hr)')

plt.subplot(1, 3, 2)
plt.plot(moment_df.query('gene_1 == "{}" & gene_2 in @mhc_genes'.format(tf)).values[:, 2:].T, '--o')
plt.xticks([0, 1, 2, 3, 4, 5], ['0', '3', '6', '9', '24', '48'])
plt.xlabel('Time (hr)'); plt.ylabel('Corr {}'.format(tf))
plt.ylim(0, 0.55)
plt.legend(moment_df.query('gene_1 == "{}" & gene_2 in @mhc_genes'.format(tf)).gene_2, frameon=False, ncol=4, bbox_to_anchor=(1.3,1.2))


select_isgs = ['MX1','MX2', 'OAS3', 'IFIT3', 'IFI2']
plt.subplot(1, 3, 3)
for i in range(4):
    plt.plot(moment_df.query('gene_1 == "{}" & gene_2 in @select_isgs'.format(tf)).values[i, 2:].T, '--o', color=plt.rcParams['axes.prop_cycle'].by_key()['color'][4+i])
plt.xticks([0, 1, 2, 3, 4, 5], ['0', '3', '6', '9', '24', '48'])
plt.xlabel('Time (hr)'); plt.ylabel('Corr {}'.format(tf))
plt.ylim(0.4, 1);
plt.legend(moment_df.query('gene_1 == "{}" & gene_2 in @select_isgs'.format(tf)).gene_2, frameon=False, ncol=2, bbox_to_anchor=(1,1))

plt.savefig(fig_path + 'mhc_mean_corr.pdf', bbox_inches='tight')
moment_df
select_isgs = ['MX1', 'OAS3', 'IFIT3', 'IFI2']
moment_df.query('gene_1 == "{}" & gene_2 in @select_isgs'.format(tf))
tf = 'STAT1'
select_isgs = ['MX1','MX2', 'OAS3', 'IFIT3', 'IFI2']
plt.figure(figsize=(2,2))
plt.plot(moment_df.query('gene_1 == "{}" & gene_2 in @select_isgs'.format(tf)).values[:, 2:].T, '--o')
plt.xticks([0, 1, 2, 3, 4, 5], ['0', '3', '6', '9', '24', '48'])
plt.xlabel('Time (hr)'); plt.ylabel('Corr {}'.format(tf))
plt.ylim(0.2, 1);
plt.legend(moment_df.query('gene_1 == "{}" & gene_2 in @select_isgs'.format(tf)).gene_2, frameon=False, ncol=2, bbox_to_anchor=(1,1))


# plt.savefig(fig_path + 'mhc_stat1_corr.pdf', bbox_inches='tight')
tf = 'NFKB2'
plt.figure(figsize=(2,2))
plt.plot(moment_df.query('gene_1 == "{}"'.format(tf)).values[:, 2:].T, '--o')
plt.xticks([0, 1, 2, 3, 4, 5], ['0', '3', '6', '9', '24', '48'])
plt.xlabel('Time (hr)'); plt.ylabel('Corr {}'.format(tf))
plt.savefig(fig_path + 'mhc_corr_NFKB2.pdf', bbox_inches='tight')

tf = 'REL'
plt.figure(figsize=(2,2))
plt.plot(moment_df.query('gene_1 == "{}"'.format(tf)).values[:, 2:].T, '--o')
plt.xticks([0, 1, 2, 3, 4, 5], ['0', '3', '6', '9', '24', '48'])
plt.xlabel('Time (hr)'); plt.ylabel('Corr {}'.format(tf))
plt.savefig(fig_path + 'mhc_corr_REL.pdf', bbox_inches='tight')

tf = 'RFX5'
plt.figure(figsize=(2,2))
plt.plot(moment_df.query('gene_1 == "{}"'.format(tf)).values[:, 2:].T, '--o')
plt.xticks([0, 1, 2, 3, 4, 5], ['0', '3', '6', '9', '24', '48'])
plt.xlabel('Time (hr)'); plt.ylabel('Corr {}'.format(tf))
tf = 'CREB3'
plt.figure(figsize=(2,2))
plt.plot(moment_df.query('gene_1 == "{}"'.format(tf)).values[:, 2:].T, '--o')
plt.xticks([0, 1, 2, 3, 4, 5], ['0', '3', '6', '9', '24', '48'])
plt.xlabel('Time (hr)'); plt.ylabel('Corr {}'.format(tf))
plt.savefig(fig_path + 'mhc_creb3_corr.pdf', bbox_inches='tight')
tf = 'STAT1'
plt.figure(figsize=(2,2))
plt.plot(moment_df.query('gene_1 == "{}"'.format(tf)).values[:, 2:].T, '--o')
plt.xticks([0, 1, 2, 3, 4, 5], ['0', '3', '6', '9', '24', '48'])
plt.xlabel('Time (hr)'); plt.ylabel('Corr {}'.format(tf))
plt.savefig(fig_path + 'mhc_stat1_corr.pdf', bbox_inches='tight')
tf = 'CREB3'
moment_df.query('gene_1 == "{}"'.format(tf))
tf = 'CREB3'
plt.figure(figsize=(2,2))
plt.plot(moment_df.query('gene_1 == "{}"'.format(tf)).values[:, 2:].T, '-o')
plt.xticks([0, 1, 2, 3, 4, 5], ['0', '3', '6', '9', '24', '48'])
plt.xlabel('Time (hr)'); plt.ylabel('Corr {}'.format(tf))
plt.legend(['HLA-A', 'HLA-B', 'HLA-C', 'B2M'], ncol=4, frameon=False, bbox_to_anchor=(1,1))
tf = 'STAT6'
plt.figure(figsize=(2,2))
plt.plot(moment_df.query('gene_1 == "{}"'.format(tf)).values[:, 2:].T, '-o')
plt.xticks([0, 1, 2, 3, 4, 5], ['0', '3', '6', '9', '24', '48'])
plt.xlabel('Time (hr)'); plt.ylabel('Corr {}'.format(tf))
tf = 'JUND'
plt.figure(figsize=(2,2))
plt.plot(moment_df.query('gene_1 == "{}"'.format(tf)).values[:, 2:].T, '-o')
plt.xticks([0, 1, 2, 3, 4, 5], ['0', '3', '6', '9', '24', '48'])
plt.xlabel('Time (hr)'); plt.ylabel('Corr {}'.format(tf))
tf = 'ETS2'
plt.figure(figsize=(2,2))
plt.plot(moment_df.query('gene_1 == "{}"'.format(tf)).values[:, 2:].T, '-o')
plt.xticks([0, 1, 2, 3, 4, 5], ['0', '3', '6', '9', '24', '48'])
plt.xlabel('Time (hr)'); plt.ylabel('Corr {}'.format(tf))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.23470580594882454, 'fraction_numerical': 0.03818677750357189, 'mean_word_length': 4.870758673274876, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 30, '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': '578159', 'license': 'MIT', 'n_tokens_mistral': 7058, 'n_tokens_neox': 6540, 'n_words': 1061, 'path': 'analysis/ifn_hbec/version2/figure_5/mhcII_analysis_pbmc.ipynb', 'repo_name': 'yelabucsf/scrna-parameter-estimation', 'size': '15398'}