<jupyter_start><jupyter_text>06/29/2020KRT (GDSD6 checks)<jupyter_code># basic packages
import os, glob
import pandas as pd
import numpy as np; np.random.seed(0)
import itertools
from collections import Counter, defaultdict
import time

# machine learning packages from sklearn
from sklearn.preprocessing import MinMaxScaler #StandardScaler 
from sklearn import preprocessing, metrics
from sklearn.feature_selection import VarianceThreshold
from sklearn.model_selection import train_test_split, KFold, cross_validate, cross_val_score, StratifiedKFold
from sklearn.linear_model import LogisticRegression, Lasso, LassoCV
from sklearn.naive_bayes import GaussianNB
from sklearn.svm import LinearSVC
from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor, VotingClassifier, AdaBoostClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.neural_network import MLPClassifier
from sklearn.metrics import roc_auc_score, auc, roc_curve, plot_roc_curve, confusion_matrix, accuracy_score
from sklearn.metrics import explained_variance_score
from scipy import interp
import scipy.stats as stats
from subprocess import call
from IPython.display import Image


# for IRF
from functools import reduce
# Needed for the scikit-learn wrapper function
import irf
from irf import (irf_utils, utils,
                 irf_jupyter_utils)
from irf.ensemble.wrf import RandomForestClassifierWithWeights
from math import ceil

# Import our custom utilities
from imp import reload


# Import tools needed for visualization
import seaborn as sns; sns.set()
import matplotlib
import matplotlib.pyplot as plt
from sklearn.tree import export_graphviz
import pydot
%load_ext autoreload
%autoreload 2
save_dir = '../data/processed/fig4_modelling/KRT_ex'
if not os.path.exists(save_dir):
    os.makedirs(save_dir)
THRES=1
normal_tissues = ['Airway','Astrocytes','Bladder','Colon','Esophageal','GDSD6','GM12878','HMEC','Melanocytes','Ovarian',
'Pancreas','Prostate','Renal','Thyroid','Uterine']
normal_tissues_dict = dict(zip(normal_tissues,range(len(normal_tissues))))
rna_df = pd.read_csv('../data/interim/rna/tissue_tpm_sym.csv',index_col=0)
rna_KRT_dict = pd.Series(rna_df.GDSD6.values, index=rna_df.index.values).to_dict()<jupyter_output><empty_output><jupyter_text>0. Data Wrangling- import- preprocess<jupyter_code># import
data_all = pd.read_csv('/Users/mguo123/Google Drive/1_khavari/omics_project-LD/pan_omics/data/processed/tissue_crms/all_count_comb_overall.csv',index_col=0,header=0)
data_all = data_all[data_all.tissue.isin(normal_tissues)]
data_all = data_all[data_all.iloc[:,2:].sum(axis=1)>1e-1]

# expression labels
exp_label = list(np.log10(data_all.exp.values+1e-2))
labels_all  = np.array(np.array(exp_label)>THRES)

tissues_label  = data_all.tissue.values#np.array((data_all.exp>THRES).values)
tissue_num_labels =  data_all.tissue.map(normal_tissues_dict).values

genes_all = data_all.index.values
gene_to_num_dict = dict(zip(np.unique(genes_all),range(len(np.unique(genes_all)))))
genes_num_all = np.vectorize(gene_to_num_dict.get)(genes_all)


print('files_loaded', data_all.shape)

data_all[:5]
## only tfs

data_all.drop(['tissue','exp','num_loop_counts','num_loops','num_atac_regions_pro','num_atac_regions_loop'],axis=1,inplace=True)

data_all.shape



selector = VarianceThreshold()
data_all_varfilt = selector.fit_transform(data_all)
data_all_varfilt_cols = data_all.columns[selector.get_support()]
print(data_all.shape, data_all_varfilt.shape, len(data_all_varfilt_cols))
scaler = MinMaxScaler()
data_all_norm = scaler.fit_transform(data_all_varfilt)
data_all_norm = pd.DataFrame(data_all_norm, columns = data_all_varfilt_cols)
data_all_norm[:5]<jupyter_output><empty_output><jupyter_text>1. See if there is even a correlation in expression between GDSD6 genes and some of the famourus TFs (AP2B, JUN, FOS, etc)<jupyter_code>KRT_genes = pd.read_csv('../../rnaseq/unique_gene_lists/'+'GDSD6'+'_genes.txt',header=None).loc[:,0]
len(KRT_genes)
KRT_tfs = sorted(set(list(pd.read_csv('../data/external/krt_tfs_063020.csv')['tfs'])))
print(KRT_tfs,len(KRT_tfs))
# data_all['tissue'] = tissues_label
tfs_feat_dict = defaultdict(list)
for feat in data_all.columns:
    tfs_feat_dict[feat.split('_')[0]].append(feat)
KRT_TF_feat = []
KRT_TF_dict = defaultdict(list)
for tf in KRT_tfs:
    if tf+'_pro' in data_all.columns:
        KRT_TF_feat.append(tf+'_pro')
        KRT_TF_dict[tf].append(tf+'_pro')
    if tf+'_loop' in data_all.columns:
        KRT_TF_feat.append(tf+'_loop')
        KRT_TF_dict[tf].append(tf+'_loop')
KRT_TF_dict['all'] = KRT_TF_feat  
KRT_crm = data_all[tissues_label=='GDSD6']
KRT_exp_arr = np.array(exp_label)[tissues_label=='GDSD6']
KRT_exp_genes_arr = set(list(KRT_crm.index[np.array(KRT_exp_arr)>THRES]))
KRT_exp_genes_num = len(KRT_exp_genes_arr)#np.where(np.array(KRT_exp_arr)>THRES)[0])
KRT_crm_KRT_genes = KRT_crm[KRT_crm.index.isin(KRT_genes)]
# num KRT expressed genes, 
KRT_exp_genes_num, KRT_crm.shape, KRT_exp_genes_num/KRT_crm.shape[0]
KRT_crm#.loc['TP63',]
KRT_tf_expr = {}
not_expr_tf=[]
for tf in KRT_tfs:
    exp = rna_KRT_dict.get(tf)
    if exp is None:
        exp=0
    KRT_tf_expr[tf]=exp
    if exp<THRES:
        not_expr_tf.append(tf)
        print(tf)
print(len(not_expr_tf))<jupyter_output>CEBPA
E2F1
FOXF2
KER2
LDB2
OTX1
PBX1
PBX2
POU1F1
POU3F1
POU3F2
POU3F3
POU3F4
POU4F1
POU4F2
POU4F3
POU5F1
POU5F1B
POU6F1
POU6F2
PRRX1
SOX11
TFAP2B
TFAP2D
TWIST2
ZEB1
26<jupyter_text>get correlation between expression and feature value<jupyter_code>corr_dict = {}
expr_genes_set = set()
for feat in data_all.columns:
    expr_genes = set(list(KRT_crm.index[KRT_crm[feat]>0])).intersection(KRT_exp_genes_arr)
    expr_genes_set = expr_genes_set.union(expr_genes)
    num_expr_genes = len(expr_genes)
    corr_dict[feat] = [np.corrcoef(KRT_exp_arr, KRT_crm[feat])[0][1], feat in KRT_TF_feat, num_expr_genes/KRT_exp_genes]
corr_df = pd.DataFrame.from_dict(corr_dict,orient='index')
corr_df.columns = ['corr', 'KRT_TF','frac_expr_genes']
corr_df['corr_sq'] = corr_df['corr'].apply(lambda x:x**2)
corr_df = corr_df.sort_values('corr_sq',ascending=False)
corr_df.reset_index(inplace=True)
corr_df.columns = ['feat', 'corr', 'KRT_TF','frac_expr_genes', 'corr_sq']
corr_df.reset_index(inplace=True)
corr_df.fillna(0,inplace=True)
print('all', len(expr_genes_set)/KRT_exp_genes)
corr_df.to_csv(os.path.join(save_dir, 'stats_corr_KRT.csv'))
corr_df.sort_values('corr_sq',ascending=False)[:50]
sns.scatterplot(x='index', y='corr_sq',alpha=1, data=corr_df[corr_df.KRT_TF])<jupyter_output><empty_output><jupyter_text>Results: we can see that the correlation between expression of KRT genes and TFs footprinting features is quite low, furthermore there seems to not be an association beween the r^2 rank and whether or not the KRT feature is contains a KRT tf (manually annotated)<jupyter_code># KRT_crm['SOX13_loop'].describe()
# corr_df#[corr_df.KRT_TF]<jupyter_output><empty_output><jupyter_text>next we see if the KRT tfs are significantly enriched in KRT genes<jupyter_code>KRT_crm.shape
results = {}
count_all = data_all.sum().sum()
count_KRTgene = KRT_crm.sum().sum() #mat_counts.sum(axis=1)[0], sum first row
for tf, feat_list in KRT_TF_dict.items():
    if len(feat_list)>0:
#         print(tf)
        KRT_crm_KRT_genes = KRT_crm[feat_list]
        count_KRTtf_KRTgene = KRT_crm_KRT_genes.sum().sum() # A
        count_KRTtf = data_all[feat_list].sum().sum() #mat_counts.sum(axis=0)[0], sum down first col
        count_KRTtf_neg = count_KRTtf - count_KRTtf_KRTgene # B
        count_neg_KRTgene = count_KRTgene - count_KRTtf_KRTgene #C
        count_neg_neg = count_all - count_KRTgene - count_KRTtf_neg #D
        mat_counts = np.array([[count_KRTtf_KRTgene,count_neg_KRTgene],
                       [count_KRTtf_neg, count_neg_neg]]).reshape((2,2))
        pseudo = 1
        mat_counts_pseudo = mat_counts+pseudo
        num_in_1 = mat_counts.sum(axis=1)[0] #count_KRTgene
        num_in_2 = mat_counts.sum(axis=0)[0] #count_KRTtf
        in_1_and_in_2 = count_KRTtf_KRTgene
        in_1_or_in_2 = count_KRTgene +count_KRTtf_neg
        in_1 = count_KRTgene
        in_2 = count_KRTtf
        observed_num = mat_counts[0][0] #count_KRTtf_KRTgene
        expected_num = num_in_1*num_in_2/sum(sum(mat_counts))
        oddsratio_pseudo, pvalue_pseudo = stats.fisher_exact(mat_counts_pseudo,alternative='greater')
        jaccard = in_1_and_in_2/in_1_or_in_2
        intersect_over_min = in_1_and_in_2/min(in_1,in_2)

        results[tf] = { 'jaccard':jaccard,'intersect_over_min':intersect_over_min,
                                'intersection':in_1_and_in_2, 
                               'union':in_1_or_in_2, 
                               'num_in_1':num_in_1,'num_in_2':num_in_2,
                               'observed':observed_num, 'expected':expected_num, 
                               'oddsratio':oddsratio_pseudo, 'pval':pvalue_pseudo}
        
result_df = pd.DataFrame.from_dict(results,orient='index')
result_df['pval_bonf'] = result_df.pval.apply(lambda x: min(1, x* sum(sum(mat_counts))))#result_df.shape[0]))
result_df['log_pval_bonf'] = result_df.pval_bonf.apply(lambda x: min(100,-np.log10(x+1e-100)))
result_df.to_csv(os.path.join(save_dir, 'stats_fisher_KRTtfs.csv'))
display(result_df[:5])
result_df.shape
result_df_filt
# # for all
# KRT_crm = data_all[tissues_label=='GDSD6']
# KRT_crm_KRT_genes = KRT_crm[KRT_TF_feat]
# count_KRTtf_KRTgene = KRT_crm_KRT_genes.sum().sum() # A
# count_all = data_all.sum().sum()
# count_KRTgene = KRT_crm.sum().sum() #mat_counts.sum(axis=1)[0], sum first row
# count_KRTtf = data_all[KRT_TF_feat].sum().sum() #mat_counts.sum(axis=0)[0], sum down first col
# count_KRTtf_neg = count_KRTtf - count_KRTtf_KRTgene # B
# count_neg_KRTgene = count_KRTgene - count_KRTtf_KRTgene #C
# count_neg_neg = count_all - count_KRTgene - count_KRTtf_neg #D


# mat_counts = np.array([[count_KRTtf_KRTgene,count_neg_KRTgene],
#                        [count_KRTtf_neg, count_neg_neg]]).reshape((2,2))
# pseudo = 1
# mat_counts_pseudo = mat_counts+pseudo
# num_in_1 = mat_counts.sum(axis=1)[0] #count_KRTgene
# num_in_2 = mat_counts.sum(axis=0)[0] #count_KRTtf
# in_1_and_in_2 = count_KRTtf_KRTgene
# in_1_or_in_2 = count_KRTgene +count_KRTtf_neg
# in_1 = count_KRTgene
# in_2 = count_KRTtf
# observed_num = mat_counts[0][0] #count_KRTtf_KRTgene
# expected_num = num_in_1*num_in_2/sum(sum(mat_counts))
# oddsratio_pseudo, pvalue_pseudo = stats.fisher_exact(mat_counts_pseudo,alternative='greater')
# jaccard = in_1_and_in_2/in_1_or_in_2
# intersect_over_min = in_1_and_in_2/min(in_1,in_2)

# results = { 'jaccard':jaccard,'intersect_over_min':intersect_over_min,
#                         'intersection':in_1_and_in_2, 
#                        'union':in_1_or_in_2, 
#                        'num_in_1':num_in_1,'num_in_2':num_in_2,
#                        'observed':observed_num, 'expected':expected_num, 
#                        'oddsratio':oddsratio_pseudo, 'pval':pvalue_pseudo}
# # results<jupyter_output><empty_output><jupyter_text>do some excel work:see pptx  next we check which tf footprints are significantly associated with having KRT genes be expressed<jupyter_code># expression
KRT_exp_crm = KRT_crm.iloc[np.where(np.array(KRT_exp_arr)>THRES)[0],:]

results = {}
count_all = KRT_crm.sum().sum()
count_KRTgene = KRT_exp_crm.sum().sum() #mat_counts.sum(axis=1)[0], sum first row
for tf, feat_list in tfs_feat_dict.items():
    if len(feat_list)>0:
#         print(tf)
        KRT_crm_KRT_genes = KRT_exp_crm[feat_list]
        count_KRTtf_KRTgene = KRT_crm_KRT_genes.sum().sum() # A
        if count_KRTtf_KRTgene==0:
            continue
        count_KRTtf = KRT_crm[feat_list].sum().sum() #mat_counts.sum(axis=0)[0], sum down first col
        count_KRTtf_neg = count_KRTtf - count_KRTtf_KRTgene # B
        count_neg_KRTgene = count_KRTgene - count_KRTtf_KRTgene #C
        count_neg_neg = count_all - count_KRTgene - count_KRTtf_neg #D
        mat_counts = np.array([[count_KRTtf_KRTgene,count_neg_KRTgene],
                       [count_KRTtf_neg, count_neg_neg]]).reshape((2,2))
        pseudo = 1
        mat_counts_pseudo = mat_counts+pseudo
        num_in_1 = mat_counts.sum(axis=1)[0] #count_KRTgene
        num_in_2 = mat_counts.sum(axis=0)[0] #count_KRTtf
        in_1_and_in_2 = count_KRTtf_KRTgene
        in_1_or_in_2 = count_KRTgene +count_KRTtf_neg
        in_1 = count_KRTgene
        in_2 = count_KRTtf
        observed_num = mat_counts[0][0] #count_KRTtf_KRTgene
        expected_num = num_in_1*num_in_2/sum(sum(mat_counts))
        oddsratio_pseudo, pvalue_pseudo = stats.fisher_exact(mat_counts_pseudo,alternative='greater')
        jaccard = in_1_and_in_2/in_1_or_in_2
        intersect_over_min = in_1_and_in_2/min(in_1,in_2)

        results[tf] = { 'jaccard':jaccard,'intersect_over_min':intersect_over_min,
                                'intersection':in_1_and_in_2, 
                               'union':in_1_or_in_2, 
                               'num_in_1':num_in_1,'num_in_2':num_in_2,
                               'observed':observed_num, 'expected':expected_num, 
                               'oddsratio':oddsratio_pseudo, 'pval':pvalue_pseudo}
        
result_df_exp = pd.DataFrame.from_dict(results,orient='index')
result_df_exp['pval_bonf'] = result_df_exp.pval.apply(lambda x: min(1, x* result_df.shape[0]))
result_df_exp['log_pval_bonf'] = result_df_exp.pval_bonf.apply(lambda x: min(100,-np.log10(x+1e-100)))
result_df_exp['is_KRT_TF'] = result_df_exp.index.isin(KRT_tfs)
result_df_exp.to_csv(os.path.join(save_dir, 'stats_fisher_KRTexp.csv'))
result_df_exp.sort_values('oddsratio',ascending=False)
result_df_exp.sort_values('oddsratio',ascending=False)[:10]


# results[result_idx] = {'color_row':color_1, 'color_col':color_2,
#                        'jaccard':jaccard,'intersect_over_min':intersect_over_min,
#                         'intersection':len(in_1_and_in_2), 
#                        'union':len(in_1_or_in_2), 
#                        'num_in_1':num_in_1,'num_in_2':num_in_2,
#                        'observed':observed_num, 'expected':expected_num, 
#                        'oddsratio':oddsratio_pseudo, 'pval':pvalue_pseudo}
# sns.distplot(KRT_crm_KRT_genes['JUND_loop'])
# def comp_two_gene_dicts(gene_dict1, gene_dict2, pseudo=0):
#     """
#     1: rows
#     2: columns
#     """
#     bg = set(gene_dict1['background']) | set(gene_dict2['background'])
    
#     results = {}
#     result_idx = 0
#     for color_1 in ['purple','green', 'grey', 'blue']:
#         for color_2 in ['purple','green', 'grey', 'blue']:
#             geneset_1 = set(gene_dict1[color_1])
#             geneset_2 = set(gene_dict2[color_2])
            
            
#             in_1_and_in_2 = geneset_1 & geneset_2
#             in_1_not_2 = geneset_1 - geneset_2
#             not_1_in_2 = geneset_2 - geneset_1
#             in_1_or_in_2 = geneset_1 | geneset_2
#             not_1_not_2 = bg - in_1_or_in_2
            
#             mat_counts = np.array([[len(in_1_and_in_2), len(in_1_not_2)],
#                                   [len(not_1_in_2), len(not_1_not_2)]]).reshape((2,2))
#             mat_counts_pseudo = mat_counts+pseudo
# #             if (color_1=='grey') & (color_2=='grey'):
# #                 print(mat_counts)
    
#             num_in_1 = mat_counts.sum(axis=1)[0]
#             num_in_2 = mat_counts.sum(axis=0)[0]
#             observed_num = mat_counts[0][0]
#             expected_num = num_in_1*num_in_2/sum(sum(mat_counts))
#             oddsratio_pseudo, pvalue_pseudo = stats.fisher_exact(mat_counts_pseudo,alternative='greater')
#             jaccard = len(in_1_and_in_2)/len(in_1_or_in_2)
#             intersect_over_min = len(in_1_and_in_2)/min(num_in_1,num_in_2)
            
#             results[result_idx] = {'color_row':color_1, 'color_col':color_2,
#                                    'jaccard':jaccard,'intersect_over_min':intersect_over_min,
#                                     'intersection':len(in_1_and_in_2), 
#                                    'union':len(in_1_or_in_2), 
#                                    'num_in_1':num_in_1,'num_in_2':num_in_2,
#                                    'observed':observed_num, 'expected':expected_num, 
#                                    'oddsratio':oddsratio_pseudo, 'pval':pvalue_pseudo}
    
#             result_idx+=1
#     result_df = pd.DataFrame.from_dict(results,orient='index')
#     result_df['pval_bonf'] = result_df.pval.apply(lambda x: min(1, x* sum(sum(mat_counts))))#result_df.shape[0]))
#     result_df['log_pval_bonf'] = result_df.pval_bonf.apply(lambda x: min(100,-np.log10(x+1e-100)))
#     return result_df<jupyter_output><empty_output><jupyter_text>2. check relationship between TFAP2C and KLF4<jupyter_code>cols_to_check=['TFAP2C_pro','TFAP2C_loop','KLF4_pro','KLF4_loop']
f = sns.pairplot(KRT_crm[cols_to_check],kind='reg')
plt.savefig(os.path.join(save_dir, 'TFAP2C_KLF4_corr_check.png'))<jupyter_output><empty_output><jupyter_text>5. Next check would be for pairs of KRT tfs to see which are significantly enriched in KRT genes<jupyter_code>tfs_w_feats = set()
for col in data_all.columns:
    tfs_w_feats.add(col.split('_')[0])
len(tfs_w_feats)  
KRT_TF_pair_dict = defaultdict(list)
for tf1 in KRT_tfs:
    for tf2 in KRT_tfs:
        if (tf1<tf2):
            if (tf1 in tfs_w_feats) and (tf2 in tfs_w_feats):
                possible_feats = [tf1+'_pro',tf1+'_loop',tf2+'_pro',tf2+'_loop']
                for feat in possible_feats:
                    if feat in data_all.columns:
                        KRT_TF_pair_dict[tf1+'::'+tf2].append(feat)
    #                 KRT_TF_pair_dict[tf1+'::'+tf2]=[]
len(KRT_TF_pair_dict)
# for key, list_feat in KRT_TF_pair_dict.items():
#     if len(list_feat)<4:
#         print(key, list_feat)
%%time
counter = 0
results = {}
count_all = data_all.sum().sum()
count_KRTgene = KRT_crm.sum().sum() #mat_counts.sum(axis=1)[0], sum first row

for tf, feat_list in KRT_TF_pair_dict.items():
    if len(feat_list)>0:
        counter+=1
        if (counter%100)==0:
            print(counter, tf)
#         print(tf)
        KRT_crm_KRT_genes = KRT_crm[feat_list]
        count_KRTtf_KRTgene = KRT_crm_KRT_genes.sum().sum() # A
        count_KRTtf = data_all[feat_list].sum().sum() #mat_counts.sum(axis=0)[0], sum down first col
        count_KRTtf_neg = count_KRTtf - count_KRTtf_KRTgene # B
        count_neg_KRTgene = count_KRTgene - count_KRTtf_KRTgene #C
        count_neg_neg = count_all - count_KRTgene - count_KRTtf_neg #D
        mat_counts = np.array([[count_KRTtf_KRTgene,count_neg_KRTgene],
                       [count_KRTtf_neg, count_neg_neg]]).reshape((2,2))
        pseudo = 1
        mat_counts_pseudo = mat_counts+pseudo
        num_in_1 = mat_counts.sum(axis=1)[0] #count_KRTgene
        num_in_2 = mat_counts.sum(axis=0)[0] #count_KRTtf
        in_1_and_in_2 = count_KRTtf_KRTgene
        in_1_or_in_2 = count_KRTgene +count_KRTtf_neg
        in_1 = count_KRTgene
        in_2 = count_KRTtf
        observed_num = mat_counts[0][0] #count_KRTtf_KRTgene
        expected_num = num_in_1*num_in_2/sum(sum(mat_counts))
        oddsratio_pseudo, pvalue_pseudo = stats.fisher_exact(mat_counts_pseudo,alternative='greater')
        jaccard = in_1_and_in_2/in_1_or_in_2
        intersect_over_min = in_1_and_in_2/min(in_1,in_2)

        results[tf] = { 'jaccard':jaccard,'intersect_over_min':intersect_over_min,
                                'intersection':in_1_and_in_2, 
                               'union':in_1_or_in_2, 
                               'num_in_1':num_in_1,'num_in_2':num_in_2,
                               'observed':observed_num, 'expected':expected_num, 
                               'oddsratio':oddsratio_pseudo, 'pval':pvalue_pseudo}
        
result_df = pd.DataFrame.from_dict(results,orient='index')
result_df['pval_bonf'] = result_df.pval.apply(lambda x: min(1, x* sum(sum(mat_counts))))#result_df.shape[0]))
result_df['log_pval_bonf'] = result_df.pval_bonf.apply(lambda x: min(100,-np.log10(x+1e-100)))
result_df.to_csv(os.path.join(save_dir, 'stats_fisher_KRTtfs_pairs.csv'))
result_df
result_df_filt = result_df[(result_df.pval_bonf<0.05) & (result_df.oddsratio>1)& (result_df.oddsratio<1000)]
print(result_df_filt.shape, result_df.shape)
result_df_filt.sort_values('oddsratio',ascending=False)[:10]
result_mat_df = pd.DataFrame(index=KRT_tfs, columns = KRT_tfs).fillna(0)
for pair in result_df_filt.index.values:
    tf1,tf2 = pair.split('::')
    result_mat_df.at[tf1,tf2]=1
    result_mat_df.at[tf2,tf1]=1
result_mat_df.to_csv(os.path.join(save_dir, 'stats_fisher_KRTtfs_pairs_mat.csv'))<jupyter_output><empty_output><jupyter_text>KRT gene pairs with expression<jupyter_code>counter = 0
results = {}
count_all = KRT_crm.sum().sum()
count_KRTgene = KRT_exp_crm.sum().sum() #mat_counts.sum(axis=1)[0], sum first row
for tf, feat_list in KRT_TF_pair_dict.items():
    if len(feat_list)>0:
        counter+=1
        if (counter%100)==0:
            print(counter, tf)
#         print(tf)
        KRT_crm_KRT_genes = KRT_exp_crm[feat_list]
        count_KRTtf_KRTgene = KRT_crm_KRT_genes.sum().sum() # A
        if count_KRTtf_KRTgene==0:
            continue
        count_KRTtf = KRT_crm[feat_list].sum().sum() #mat_counts.sum(axis=0)[0], sum down first col
        count_KRTtf_neg = count_KRTtf - count_KRTtf_KRTgene # B
        count_neg_KRTgene = count_KRTgene - count_KRTtf_KRTgene #C
        count_neg_neg = count_all - count_KRTgene - count_KRTtf_neg #D
        mat_counts = np.array([[count_KRTtf_KRTgene,count_neg_KRTgene],
                       [count_KRTtf_neg, count_neg_neg]]).reshape((2,2))
        pseudo = 1
        mat_counts_pseudo = mat_counts+pseudo
        num_in_1 = mat_counts.sum(axis=1)[0] #count_KRTgene
        num_in_2 = mat_counts.sum(axis=0)[0] #count_KRTtf
        in_1_and_in_2 = count_KRTtf_KRTgene
        in_1_or_in_2 = count_KRTgene +count_KRTtf_neg
        in_1 = count_KRTgene
        in_2 = count_KRTtf
        observed_num = mat_counts[0][0] #count_KRTtf_KRTgene
        expected_num = num_in_1*num_in_2/sum(sum(mat_counts))
        oddsratio_pseudo, pvalue_pseudo = stats.fisher_exact(mat_counts_pseudo,alternative='greater')
        jaccard = in_1_and_in_2/in_1_or_in_2
        intersect_over_min = in_1_and_in_2/min(in_1,in_2)

        results[tf] = { 'jaccard':jaccard,'intersect_over_min':intersect_over_min,
                                'intersection':in_1_and_in_2, 
                               'union':in_1_or_in_2, 
                               'num_in_1':num_in_1,'num_in_2':num_in_2,
                               'observed':observed_num, 'expected':expected_num, 
                               'oddsratio':oddsratio_pseudo, 'pval':pvalue_pseudo}
        
result_df = pd.DataFrame.from_dict(results,orient='index')
result_df['pval_bonf'] = result_df.pval.apply(lambda x: min(1, x* sum(sum(mat_counts))))#result_df.shape[0]))
result_df['log_pval_bonf'] = result_df.pval_bonf.apply(lambda x: min(100,-np.log10(x+1e-100)))
result_df.to_csv(os.path.join(save_dir, 'stats_fisher_KRTtfs_exp_pairs.csv'))
result_df
result_df_filt = result_df[(result_df.pval<0.05) & (result_df.oddsratio>1)& (result_df.oddsratio<1000)]
print(result_df_filt.shape, result_df.shape)
result_df_filt.sort_values('oddsratio',ascending=False)[:50]
result_df_filt.loc['KLF4::TFAP2C',:]<jupyter_output><empty_output><jupyter_text>5. Next check stability score found interactions and annotate  of KRT tfs each interaction has<jupyter_code>score_thres=.2
def check_interaction(interaction_str):
    feat_arr = interaction_str.split('::')
    num_KRT = 0
    num_tot = len(feat_arr)
    for feat in feat_arr:
        if feat in KRT_TF_feat:
            num_KRT+=1
    return [num_KRT,num_tot]
# GDSD6_stability_df = pd.read_csv('../data/processed/fig4_modelling/irf_manual/test_GDSD6_GDSD6_boosted_stability_score.csv')
GDSD6_stability_df = pd.read_csv('../data/processed/fig4_modelling/irf_manual/test_GDSD6_purple_boosted_stability_score.csv')
GDSD6_stability_df
GDSD6_stability_df[['num_KRT','num_in_interaction']]=GDSD6_stability_df['index'].apply(func=check_interaction).apply(pd.Series)#,result_type='expand')

GDSD6_stability_df['frac_KRT'] = GDSD6_stability_df['num_KRT']/GDSD6_stability_df['num_in_interaction']
sns.pairplot(GDSD6_stability_df[['score','frac_KRT']],kind='reg')
GDSD6_stability_df[(GDSD6_stability_df.score >score_thres)]

GDSD6_stability_df[(GDSD6_stability_df.frac_KRT>0) & (GDSD6_stability_df.score >score_thres)]<jupyter_output><empty_output><jupyter_text>SKIN trajectory analysis<jupyter_code># rna
rna_df = pd.read_csv('../data/interim/rna/tissue_tpm_sym.csv',index_col=0)
rna_df_norm = rna_df[normal_tissues]
rna_D0_dict = pd.Series(rna_df.GDSD0.values, index=rna_df.index.values).to_dict()
rna_D3_dict = pd.Series(rna_df.GDSD3.values, index=rna_df.index.values).to_dict()
rna_D6_dict = pd.Series(rna_df.GDSD6.values, index=rna_df.index.values).to_dict()
%%time
D0_crm = pd.read_csv('/Users/mguo123/Google Drive/1_khavari/omics_project-LD/pan_omics/data/processed/tissue_crms/pro_loop_tissue/GDSD0_crm.csv',index_col=0,header=0)
D0_crm.drop(['tissue','exp','num_loop_counts','num_loops','num_atac_regions_pro','num_atac_regions_loop'],axis=1,inplace=True)

D3_crm = pd.read_csv('/Users/mguo123/Google Drive/1_khavari/omics_project-LD/pan_omics/data/processed/tissue_crms/pro_loop_tissue/GDSD3_crm.csv',index_col=0,header=0)
D3_crm.drop(['tissue','exp','num_loop_counts','num_loops','num_atac_regions_pro','num_atac_regions_loop'],axis=1,inplace=True)

D6_crm = pd.read_csv('/Users/mguo123/Google Drive/1_khavari/omics_project-LD/pan_omics/data/processed/tissue_crms/pro_loop_tissue/GDSD6_crm.csv',index_col=0,header=0)
D6_crm.drop(['tissue','exp','num_loop_counts','num_loops','num_atac_regions_pro','num_atac_regions_loop'],axis=1,inplace=True)

what genes

D0_crm<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup