# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # 8F_controls
#
# 08/25/2020
# - positive controls = common loops, common ATAC, and expressed genes across all types
# - negative controls = no loop, no ATAC, no expression across all cell types 
#
#
#
#
#
#

from collections import Counter, defaultdict
import pandas as pd
import numpy as np; np.random.seed(0)
import os, glob
import seaborn as sns; sns.set()
import matplotlib
import matplotlib.pyplot as plt
import pybedtools
from scipy.stats import ttest_ind

# %load_ext autoreload
# %autoreload 2

save_dir = '../data/processed/fig4_modelling/vocab_sum_final/'
if not os.path.exists(save_dir):
    os.makedirs(save_dir)

data_all = pd.read_csv('/Users/mguo123/Google Drive/1_khavari/omics_project-LD/pan_omics/data/processed/tissue_crms/all_count_sep_overall.csv',index_col=0,header=0)


print(len(data_all.tissue.unique()))
data_all.tissue.unique()

# +
# # 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)

# genes_all = data_all.index.values

# -

# # negative controls
# no loop, no ATAC, no expression across all cell types for negative controls
#
# for the negative controls, narrow down the list by loci that don’t have an ATAC peak within a certain distance: Like 100 KB or some such? I
#
# make sure actually in a gene desert or region of inaccessible chromatin so really should be “off” / negative
#

NEG_DIST_MIN = 25000 # for negative control, gene must by at least 10kb away from closest atac peak

# %%time
neg_crms = data_all[data_all.iloc[:,1:].sum(axis=1)<0.001]

print(data_all.shape)
print(neg_crms.shape)

gene_tissue_counts = pd.Series(neg_crms.index).value_counts()
neg_genes = gene_tissue_counts[gene_tissue_counts==len(data_all.tissue.unique())].index
len(neg_genes) # number of genes with no loop, ATAC or expression (looped or proximal)

gene_regions = pybedtools.BedTool('../data/external/gencode.v19.gene.bed').to_dataframe()
gene_regions_neg = gene_regions[gene_regions.name.isin(neg_genes)]
print(gene_regions_neg.shape) # # genes with genomic location
gene_regions_neg_bed = pybedtools.BedTool.from_dataframe(gene_regions_neg).sort()

print(len(glob.glob('../data/interim/merged/atac/*bed')))
sorted(glob.glob('../data/interim/merged/atac/*bed'))

# %%time
all_far_genes = []
all_far_genes_dict = {}
for atac_bed in sorted(glob.glob('../data/interim/merged/atac/*bed')):
    tissue = os.path.basename(atac_bed).split('_merged')[0]
    print(tissue)
    gene_regions_neg_bed_tissue = gene_regions_neg_bed.closest(atac_bed,d=True)

    gene_dist_df = gene_regions_neg_bed_tissue.to_dataframe().groupby('name').agg({'thickEnd':min})
    gene_dist_df.columns = ['dist']
    far_genes_tissue = list(gene_dist_df.index[gene_dist_df.dist > NEG_DIST_MIN])
    print(len(far_genes_tissue))
    all_far_genes+=far_genes_tissue
    all_far_genes_dict[tissue] = far_genes_tissue

far_genes_count = pd.Series(Counter(all_far_genes))
neg_genes_far = sorted(far_genes_count.index[far_genes_count==len(glob.glob('../data/interim/merged/atac/*bed'))])
print(len(neg_genes_far))

print(neg_genes_far)

pd.Series(neg_genes_far).to_csv(os.path.join(save_dir, 'neg_genes_far.csv'),index=None,header=None)

# # 2. positive controls
#
# positive controls = common loops, common ATAC, and expressed genes across all types

# +
THRES=1 # for RNA
normal_tissues_pos = [ 'Airway', 'Astrocytes', 'Bladder', 'Colon', 
         'Esophageal', 'GDSD0', 'GDSD3', 'GDSD6', 'GM12878', 'HMEC', 'Melanocytes', 'Ovarian', 'Pancreas',
       'Prostate', 'Renal', 'Thyroid',       'Uterine']
cancer_tissues_pos = ['A431-CTRLi', 'CAL27-CTRLi', 'SCC13-CTRLi','COLO_SCR_DMSO','WM_SCR_DMSO',
                     'A431-p63i','CAL27-p63i', 'SCC13-p63i', 'COLO_SCR_PLX','COLO_shMITF_DMSO', 'COLO_shMITF_PLX']

all_tissue_pos = normal_tissues_pos + ['A431-CTRLi', 'CAL27-CTRLi', 'SCC13-CTRLi','COLO_SCR_DMSO','WM_SCR_DMSO']
#                'D0-CTRLi','D3-CTRLi', 'SCC13-p63i', 'D0-p63i', 'D3-p63i',
#                'A431-CTRLi', 'CAL27-CTRLi',  'SCC13-CTRLi','COLO_SCR_DMSO','WM_SCR_DMSO']

# +
# expressed_genes

expr_df = data_all[data_all.tissue.isin(all_tissue_pos)]['exp'].reset_index()
expr_df.columns = ['gene','exp']
expr_df_min = expr_df.groupby('gene').agg({'exp':'min'}) # basically find lowest expressed gene 
expr_genes = expr_df_min.index[expr_df_min.exp>THRES]
print(len(expr_genes)) # number of commonly expressed genes

#common ATAC
atac_pro_df = data_all[data_all.tissue.isin(all_tissue_pos)].num_atac_regions_pro.reset_index()
atac_pro_df.columns = ['gene','atac']
atac_pro_df_min = atac_pro_df.groupby('gene').agg({'atac':'min'}) # basically find lowest # of atac peaks in the promoter region gene 
atac_genes = atac_pro_df_min.index[atac_pro_df_min.atac>0]
print(len(atac_genes)) # number of atac genes that have at least 1 atac peak in promoter region (2kb upstream, 500 bp downstream of tss) in all tissues

# common loop
hichip_pro_df = data_all[data_all.tissue.isin(all_tissue_pos)].num_loops.reset_index()
hichip_pro_df.columns = ['gene','hichip']
hichip_pro_df_min = hichip_pro_df.groupby('gene').agg({'hichip':'min'}) # basically find lowest # of atac peaks in the promoter region gene 
hichip_genes = hichip_pro_df_min.index[hichip_pro_df_min.hichip>0]
print(len(hichip_genes)) # number of atac genes that have at least 1 loop to promoter region (2kb upstream, 500 bp downstream of tss) in all tissues

all_pos_genes = sorted(set(expr_genes).intersection(set(atac_genes)).intersection(set(hichip_genes)))
len(all_pos_genes)
# -

expr_genes

all_pos_genes

pd.Series(all_pos_genes).to_csv(os.path.join(save_dir, 'all_pos_genes.csv'),index=None,header=None)

# +
# expressed_genes

expr_df = data_all[data_all.tissue.isin(normal_tissues_pos)]['exp'].reset_index()
expr_df.columns = ['gene','exp']
expr_df_min = expr_df.groupby('gene').agg({'exp':'min'}) # basically find lowest expressed gene 
expr_genes = expr_df_min.index[expr_df_min.exp>THRES]
print(len(expr_genes)) # number of commonly expressed genes

#common ATAC
atac_pro_df = data_all[data_all.tissue.isin(normal_tissues_pos)].num_atac_regions_pro.reset_index()
atac_pro_df.columns = ['gene','atac']
atac_pro_df_min = atac_pro_df.groupby('gene').agg({'atac':'min'}) # basically find lowest # of atac peaks in the promoter region gene 
atac_genes = atac_pro_df_min.index[atac_pro_df_min.atac>0]
print(len(atac_genes)) # number of atac genes that have at least 1 atac peak in promoter region (2kb upstream, 500 bp downstream of tss) in all tissues

# common loop
hichip_pro_df = data_all[data_all.tissue.isin(normal_tissues_pos)].num_loops.reset_index()
hichip_pro_df.columns = ['gene','hichip']
hichip_pro_df_min = hichip_pro_df.groupby('gene').agg({'hichip':'min'}) # basically find lowest # of atac peaks in the promoter region gene 
hichip_genes = hichip_pro_df_min.index[hichip_pro_df_min.hichip>0]
print(len(hichip_genes)) # number of atac genes that have at least 1 loop to promoter region (2kb upstream, 500 bp downstream of tss) in all tissues

normal_pos_genes = list(set(expr_genes).intersection(set(atac_genes)).intersection(set(hichip_genes)))
len(normal_pos_genes)

# +
# expressed_genes

expr_df = data_all[data_all.tissue.isin(cancer_tissues_pos)]['exp'].reset_index()
expr_df.columns = ['gene','exp']
expr_df_min = expr_df.groupby('gene').agg({'exp':'min'}) # basically find lowest expressed gene 
expr_genes = expr_df_min.index[expr_df_min.exp>THRES]
print(len(expr_genes)) # number of commonly expressed genes

#common ATAC
atac_pro_df = data_all[data_all.tissue.isin(cancer_tissues_pos)].num_atac_regions_pro.reset_index()
atac_pro_df.columns = ['gene','atac']
atac_pro_df_min = atac_pro_df.groupby('gene').agg({'atac':'min'}) # basically find lowest # of atac peaks in the promoter region gene 
atac_genes = atac_pro_df_min.index[atac_pro_df_min.atac>0]
print(len(atac_genes)) # number of atac genes that have at least 1 atac peak in promoter region (2kb upstream, 500 bp downstream of tss) in all tissues

# common loop
hichip_pro_df = data_all[data_all.tissue.isin(cancer_tissues_pos)].num_loops.reset_index()
hichip_pro_df.columns = ['gene','hichip']
hichip_pro_df_min = hichip_pro_df.groupby('gene').agg({'hichip':'min'}) # basically find lowest # of atac peaks in the promoter region gene 
hichip_genes = hichip_pro_df_min.index[hichip_pro_df_min.hichip>0]
print(len(hichip_genes)) # number of atac genes that have at least 1 loop to promoter region (2kb upstream, 500 bp downstream of tss) in all tissues


cancer_pos_genes = list(set(expr_genes).intersection(set(atac_genes)).intersection(set(hichip_genes)))
len(cancer_pos_genes)
# -

len(set(normal_pos_genes).union(set(cancer_pos_genes)))

pos_genes_df = pd.concat([pd.DataFrame({'gene':normal_pos_genes,'type':'normal'}),
           pd.DataFrame({'gene':cancer_pos_genes,'type':'cancer'})],axis=0)
pos_genes_df = pos_genes_df.groupby('gene').agg({'type':'|'.join }).reset_index()
pos_genes_df.to_csv(os.path.join(save_dir, 'pos_genes_df.csv'))

pos_genes_df



# # 3. Check Controls for Motifs (MOODS
#
# 09/09/2020
#
# after 8F1_make_mpra
#

from scipy.stats import ttest_ind
import pandas as pd
import numpy as np; np.random.seed(0)
import os, glob
import seaborn as sns; sns.set()
import matplotlib
import matplotlib.pyplot as plt
from collections import Counter

glob.glob('../data/processed/fig4_modelling/vocab_sum_final/mpra*')

# read rna
rna_df = pd.read_csv('../data/interim/rna/tissue_tpm_sym.csv',index_col=0)
print(rna_df.columns)

# read mpra info 
mpra_df = pd.read_csv( '../data/processed/fig4_modelling/vocab_sum_final/mpra_oligo_df_final.txt', sep='\t')
name_to_seq_type = pd.Series(mpra_df.seq_type.values, index=mpra_df.name.values).to_dict()
# want seq_type 'pos' or 'neg'
len(name_to_seq_type)

# check rna tpm values of the pos controls 
mpra_df_controls = mpra_df[mpra_df.seq_type.isin(['pos','neg'])]
mpra_df_controls_instance = mpra_df_controls[['name','seq_type']].drop_duplicates().reset_index(drop=True)
print(mpra_df_controls_instance.shape)
mpra_df_controls_instance['Colon_exp'] = mpra_df_controls_instance.name.apply(lambda x: rna_df.loc[x, 'Colon'])
mpra_df_controls_instance['GDSD6_exp'] = mpra_df_controls_instance.name.apply(lambda x: rna_df.loc[x, 'GDSD6'])
mpra_df_controls_instance['GM12878_exp'] = mpra_df_controls_instance.name.apply(lambda x: rna_df.loc[x, 'GM12878'])
mpra_df_controls_instance['Melanocytes_exp'] = mpra_df_controls_instance.name.apply(lambda x: rna_df.loc[x, 'Melanocytes'])
mpra_df_controls_instance['A431-CTRLi_exp'] = mpra_df_controls_instance.name.apply(lambda x: rna_df.loc[x, 'A431-CTRLi'])
mpra_df_controls_instance['CAL27-CTRLi_exp'] = mpra_df_controls_instance.name.apply(lambda x: rna_df.loc[x, 'CAL27-CTRLi'])
mpra_df_controls_instance['SCC13-CTRLi_exp'] = mpra_df_controls_instance.name.apply(lambda x: rna_df.loc[x, 'SCC13-CTRLi'])
mpra_df_controls_instance['COLO_SCR_DMSO_exp'] = mpra_df_controls_instance.name.apply(lambda x: rna_df.loc[x, 'COLO_SCR_DMSO'])
mpra_df_controls_instance['WM_SCR_DMSO_exp'] = mpra_df_controls_instance.name.apply(lambda x: rna_df.loc[x, 'WM_SCR_DMSO'])
mpra_df_controls_instance.sort_values(['seq_type','name'],inplace=True)
mpra_df_pos_controls_instance = mpra_df_controls_instance[mpra_df_controls_instance.seq_type=='pos'].reset_index(drop=True)
display(mpra_df_pos_controls_instance.describe())

name_to_seq_type_count = Counter(name_to_seq_type.values())
name_to_seq_type_count

#read tf info
tf_annon_df = pd.read_csv('../data/external/HOCOMOCOv11_annotation.csv',index_col=0)
tf_annon_df['id_trim'] = tf_annon_df['id'] + '.pwm.trim'
tf_name_to_id_dict = pd.Series(tf_annon_df.id_trim.values, index=tf_annon_df.tf.values).to_dict()
tf_id_to_name_dict = pd.Series(tf_annon_df.tf.values, index=tf_annon_df.id_trim.values).to_dict()
tf_name_to_id_abbr_dict = pd.Series(tf_annon_df.id.values, index=tf_annon_df.tf.values).to_dict()
tf_id_abbr_to_name_dict = pd.Series(tf_annon_df.tf.values, index=tf_annon_df.id.values).to_dict()


# motif_scan_control_dir = '../data/processed/fig4_modelling/vocab_sum_final/motif_scan_control/'
motif_scan_control_dir = '../data/processed/fig4_modelling/vocab_sum_final/motif_scan_control_pval001/'
motif_scan_control_files = glob.glob(os.path.join(motif_scan_control_dir, '*.csv'))
print(len(motif_scan_control_files))


# +
def get_str_arr(l):
    return '|'.join([str(x) for x in l])

def ttest_motif(row):
    list_neg = [float(x) for x in row['score_get_str_arr_neg'].split('|')]
    list_pos = [float(x) for x in row['score_get_str_arr_pos'].split('|')]
    t, p = ttest_ind(list_pos, list_neg)
    return t,p


# +
# %%time
pos_neg_scan_info_all=pd.DataFrame()
pos_neg_instance_scan_info_all=pd.DataFrame()

# for motif_scan_control_file in motif_scan_control_files[:5]: ## DEBUG
for motif_scan_control_file in motif_scan_control_files:

    motif_scan_results = pd.read_csv(motif_scan_control_file,names=['name_config','pwm_file','pos','strand','score','seq','seq_actual'])
    motif = tf_id_abbr_to_name_dict[os.path.basename(motif_scan_control_file).split('_scan')[0]]
    motif_scan_results['motif'] = motif_scan_results.pwm_file.str.split('.pwm',expand=True).iloc[:,0].map(tf_id_abbr_to_name_dict)
    motif_scan_results['seq_type']  = motif_scan_results.name_config.map(name_to_seq_type)
    
    # create  df, name_config,seq_type motif, avg score
    instance_scan_info = motif_scan_results.groupby(['name_config','seq_type','motif']).agg({'score':['count','mean','std','min','max']})
    instance_scan_info.columns = ['_'.join(x) for x in instance_scan_info.columns ]
    instance_scan_info.reset_index(inplace=True)
    pos_neg_instance_scan_info_all = pd.concat([pos_neg_instance_scan_info_all, instance_scan_info],axis=0)

    # # create  df, seq_type motif, avg score
    pos_neg_scan_info = motif_scan_results.groupby(['seq_type','motif']).agg({'score':['count','mean','std','min','max', get_str_arr ]})
    pos_neg_scan_info.columns = ['_'.join(x) for x in pos_neg_scan_info.columns ]
    pos_neg_scan_info.reset_index(inplace=True)
    pos_neg_scan_info_all = pd.concat([pos_neg_scan_info_all, pos_neg_scan_info],axis=0)
    

# -

# ### check to see which motifs are enriched in controls, and do a t test to compare pos vs neg (MOODS)

# there are duplicate motif values
print(pos_neg_scan_info_all.shape)
pos_neg_scan_info_all = pos_neg_scan_info_all.groupby(['seq_type','motif']).agg({'score_count':sum,
                                                        'score_mean': 'mean',
                                                        'score_std': 'mean',
                                                         'score_min': min,
                                                         'score_max':max,
                                                         'score_get_str_arr': lambda x: '|'.join(x)
                                                        }).reset_index()
print(pos_neg_scan_info_all.shape)


pos_neg_scan_info_wide = pos_neg_scan_info_all.pivot(index='motif',columns='seq_type',values=['score_count','score_mean','score_std','score_min','score_max','score_get_str_arr'])
pos_neg_scan_info_wide.columns = ['_'.join(x) for x in pos_neg_scan_info_wide.columns ]
pos_neg_scan_info_wide['mean_pos_neg_range'] = pos_neg_scan_info_wide.score_mean_pos - pos_neg_scan_info_wide.score_mean_neg
pos_neg_scan_info_wide = pd.concat([pos_neg_scan_info_wide, pos_neg_scan_info_wide.apply(ttest_motif,axis=1).apply(lambda s: pd.Series({'t_val':s[0], 'p_val':s[1]}))],axis=1)
pos_neg_scan_info_wide['p_val_bonf'] = pos_neg_scan_info_wide.p_val.apply(lambda x: min (1, x*pos_neg_scan_info_wide.shape[0]))

print(pos_neg_scan_info_wide.shape)
pos_neg_scan_info_wide_filt = pos_neg_scan_info_wide[pos_neg_scan_info_wide.p_val_bonf<0.05]
print(pos_neg_scan_info_wide_filt.shape)

motif_count_pos_neg_ratio = pd.Series(pos_neg_scan_info_wide.score_count_pos/pos_neg_scan_info_wide.score_count_neg*name_to_seq_type_count['neg']/name_to_seq_type_count['pos']).sort_values()
print(motif_count_pos_neg_ratio.shape)
motif_count_pos_neg_ratio[motif_count_pos_neg_ratio>1].shape, motif_count_pos_neg_ratio[motif_count_pos_neg_ratio<1].shape


pos_neg_scan_info_wide_filt[pos_neg_scan_info_wide_filt.t_val<0].index.values

pos_neg_scan_info_wide_filt[pos_neg_scan_info_wide_filt.t_val>0].index.values

pos_neg_scan_info_wide_filt[pos_neg_scan_info_wide_filt.t_val<0].index.values

# ## try to find vocab pairs within a single instance

vocab_pairs = mpra_df[mpra_df.seq_type=='vocab'].name.apply(lambda x:x.split('+')[0]).unique()
len(vocab_pairs)

pos_neg_instance_scan_info_all[:5]

neg_instance_scan_info_all = pos_neg_instance_scan_info_all[pos_neg_instance_scan_info_all.seq_type=='neg']

neg_instance_motifs = neg_instance_scan_info_all.groupby('name_config').agg({'motif':'|'.join}).reset_index()


def check_vocab_pair(motif_str, vocab_pairs=vocab_pairs):
    motif_list = motif_str.split('|')
    vocab_pairs_found = []
    for vocab_pair in vocab_pairs:
        v1, v2 = vocab_pair.split('::')
        if (v1 in motif_list) & (v2 in motif_list):
            vocab_pairs_found.append(vocab_pair)
    return '|'.join(sorted(vocab_pairs_found)),len(vocab_pairs_found)


neg_instance_motifs = pd.concat([neg_instance_motifs,
                                 neg_instance_motifs.motif.apply(check_vocab_pair).apply(lambda s: pd.Series({'vocab_pairs_found':s[0], 'num_vocab_pairs_found':s[1]}))],axis=1)
neg_instance_motifs

neg_instance_motifs_counts = Counter('|'.join(list(neg_instance_motifs.vocab_pairs_found)).split('|'))

pos_instance_scan_info_all = pos_neg_instance_scan_info_all[pos_neg_instance_scan_info_all.seq_type=='pos']
pos_instance_motifs = pos_instance_scan_info_all.groupby('name_config').agg({'motif':'|'.join}).reset_index()
pos_instance_motifs = pd.concat([pos_instance_motifs,
                                 pos_instance_motifs.motif.apply(check_vocab_pair).apply(lambda s: pd.Series({'vocab_pairs_found':s[0], 'num_vocab_pairs_found':s[1]}))],axis=1)
pos_instance_motifs_counts = Counter('|'.join(list(pos_instance_motifs.vocab_pairs_found)).split('|'))

control_vocab_pair_dict = {}
for vocab in list(neg_instance_motifs_counts.keys()) + list(pos_instance_motifs_counts.keys()):
    control_vocab_pair_dict[vocab] = {'in_pos':vocab in pos_instance_motifs_counts,
                                      'in_neg':vocab in neg_instance_motifs_counts,
                                      'count_pos':pos_instance_motifs_counts.get(vocab, 0),
                                      'count_neg':neg_instance_motifs_counts.get(vocab, 0),
        
    }

control_vocab_pair_df = pd.DataFrame.from_dict(control_vocab_pair_dict,orient='index')
control_vocab_pair_df['pos_neg_ratio'] = (control_vocab_pair_df.count_pos+1)/(control_vocab_pair_df.count_neg+1)
control_vocab_pair_df['pos_neg_ratio_corr'] = control_vocab_pair_df.pos_neg_ratio*neg_instance_motifs.shape[0]/pos_instance_motifs.shape[0]
control_vocab_pair_df.sort_values('pos_neg_ratio_corr',inplace=True)
control_vocab_pair_df.describe()







# # 4. AME enrichment

# ## AME motif enrichment for negative controls
#
# `ame --verbose 1 --oc . --scoring avg --method fisher --hit-lo-fraction 0.25 --evalue-report-threshold 10.0 --control --shuffle-- --kmer 2 mpra_neg_control_seq.fasta motifs.meme`
#
# AME motif enrichment neg control
#
# motifs: '../data/external/hocomoco_human_trim_jaspar_format.txt'

save_dir

# num actual tfs used 
vocab_pairs = mpra_df[mpra_df.seq_type=='vocab'].name.apply(lambda x:x.split('+')[0]).unique()
list_tfs = []
for vocab_pair in vocab_pairs:
    vocab1, vocab2 = vocab_pair.split('::')
    list_tfs.append(vocab1)
    list_tfs.append(vocab2)
list_tfs = sorted(set(list_tfs))
len(list_tfs)

#get motif annotation file
tf_annon_df = pd.read_csv('../data/external/HOCOMOCOv11_annotation.csv',index_col=0)
tf_annon_df['id_trim'] = tf_annon_df['id'] + '.pwm.trim'
tf_name_to_id_dict = pd.Series(tf_annon_df.id_trim.values, index=tf_annon_df.tf.values).to_dict()
tf_id_to_name_dict = pd.Series(tf_annon_df.tf.values, index=tf_annon_df.id_trim.values).to_dict()
tf_name_to_id_abbr_dict = pd.Series(tf_annon_df.id.values, index=tf_annon_df.tf.values).to_dict()
tf_id_abbr_to_name_dict = pd.Series(tf_annon_df.tf.values, index=tf_annon_df.id.values).to_dict()


#get motif lengths
motif_to_len = {}
for file in sorted(glob.glob('../data/external/hocomoco_human_trim/*')):
    motif_id = os.path.basename(file).split('.pwm')[0]
#     print(tf_id_abbr_to_name_dict[motif_id])
    with open(file, 'r') as g:
        len_motif = len(g.readlines()[0].strip().split())
        motif_to_len[tf_id_abbr_to_name_dict[motif_id]] = len_motif

len(set(tf_id_abbr_to_name_dict.values())), len(motif_to_len)

#read in motif enrichment results
ame_result_file = os.path.join(save_dir, 'motif_scan_neg_control_ame.tsv')
ame_result = pd.read_csv(ame_result_file, sep='\t')
ame_result['motif_DB'] = 'hocomoco_human_trim'
ame_result['motif'] = ame_result.motif_ID.map(tf_id_abbr_to_name_dict)
ame_result = ame_result[~ame_result.motif.isna()]
ame_result['in_vocab'] = ame_result.motif.isin(list_tfs)

ame_result

# #### these are the motifs that are enriched in the negative control sequences

enriched_motifs_in_neg = sorted(ame_result.motif.unique())
print(len(enriched_motifs_in_neg))
print(enriched_motifs_in_neg)

# #### these are the motifs that are enriched in the negative control sequences AND are part of vocab words

ame_result_filt = ame_result[ame_result.in_vocab]
print(ame_result.shape, ame_result_filt.shape)
enriched_motifs_in_neg_vocab = sorted(ame_result_filt.motif.unique())
print(len(enriched_motifs_in_neg_vocab))
print(enriched_motifs_in_neg_vocab)

# ## AME motif enrichment negative controls (background vocab words)
# `ame --verbose 1 --oc . --scoring avg --method fisher --hit-lo-fraction 0.25 --evalue-report-threshold 10.0 --control mpra_vocab_seq2.fasta mpra_neg_control_seq.fasta motifs.meme`

# +
#read in motif enrichment results
ame_result_file = os.path.join(save_dir, 'motif_scan_neg_control_bckgdvocab_ame.tsv')
ame_result = pd.read_csv(ame_result_file, sep='\t')
ame_result['motif_DB'] = 'hocomoco_human_trim'
ame_result['motif'] = ame_result.motif_ID.map(tf_id_abbr_to_name_dict)
ame_result = ame_result[~ame_result.motif.isna()]
ame_result['in_vocab'] = ame_result.motif.isin(list_tfs)

print('enriched motifs')
enriched_motifs = sorted(ame_result.motif.unique())
print(len(enriched_motifs))
print(enriched_motifs)

print('enriched motifs with tested vocab pair')
ame_result_filt = ame_result[ame_result.in_vocab]
print(ame_result.shape, ame_result_filt.shape)
enriched_motifs_vocab = sorted(ame_result_filt.motif.unique())
print(len(enriched_motifs_vocab))
print(enriched_motifs_vocab)
# -





# ## AME motif enrichment for positive controls
#
# `ame --verbose 1 --oc . --scoring avg --method fisher --hit-lo-fraction 0.25 --evalue-report-threshold 10.0 --control --shuffle-- --kmer 2 mpra_pos_control_seq.fasta motifs.meme`
#
# AME motif enrichment neg control
#
# motifs: '../data/external/hocomoco_human_trim_jaspar_format.txt'

# +
#read in motif enrichment results
ame_result_file = os.path.join(save_dir, 'motif_scan_pos_control_ame.tsv')
ame_result = pd.read_csv(ame_result_file, sep='\t')
ame_result['motif_DB'] = 'hocomoco_human_trim'
ame_result['motif'] = ame_result.motif_ID.map(tf_id_abbr_to_name_dict)
ame_result = ame_result[~ame_result.motif.isna()]
ame_result['in_vocab'] = ame_result.motif.isin(list_tfs)

print('enriched motifs')
enriched_motifs = sorted(ame_result.motif.unique())
print(len(enriched_motifs))
print(enriched_motifs)

print('enriched motifs with tested vocab pair')
ame_result_filt = ame_result[ame_result.in_vocab]
print(ame_result.shape, ame_result_filt.shape)
enriched_motifs_vocab = sorted(ame_result_filt.motif.unique())
print(len(enriched_motifs_vocab))
print(enriched_motifs_vocab)
# -

# ## AME motif enrichment for positive controls background vocab words
#
# `ame --verbose 1 --oc . --scoring avg --method fisher --hit-lo-fraction 0.25 --evalue-report-threshold 10.0 --control mpra_vocab_seq2.fasta mpra_pos_control_seq.fasta motifs.meme`
#
# AME motif enrichment neg control
#
# motifs: '../data/external/hocomoco_human_trim_jaspar_format.txt'

# +
#read in motif enrichment results
ame_result_file = os.path.join(save_dir, 'motif_scan_pos_control_bckgdvocab_ame.tsv')
ame_result = pd.read_csv(ame_result_file, sep='\t')
ame_result['motif_DB'] = 'hocomoco_human_trim'
ame_result['motif'] = ame_result.motif_ID.map(tf_id_abbr_to_name_dict)
ame_result = ame_result[~ame_result.motif.isna()]
ame_result['in_vocab'] = ame_result.motif.isin(list_tfs)

print('enriched motifs')
enriched_motifs = sorted(ame_result.motif.unique())
print(len(enriched_motifs))
print(enriched_motifs)

print('enriched motifs with tested vocab pair')
ame_result_filt = ame_result[ame_result.in_vocab]
print(ame_result.shape, ame_result_filt.shape)
enriched_motifs_vocab = sorted(ame_result_filt.motif.unique())
print(len(enriched_motifs_vocab))
print(enriched_motifs_vocab)
# -

# ## AME motif enrichment for vocab words
#
# `ame --verbose 1 --oc . --scoring avg --method fisher --hit-lo-fraction 0.25 --evalue-report-threshold 10.0 --control --shuffle-- --kmer 2 mpra_vocab_seq2.fasta motifs.meme
# `
#
# AME motif enrichment neg control
#
# motifs: '../data/external/hocomoco_human_trim_jaspar_format.txt'

# +
#read in motif enrichment results
ame_result_file = os.path.join(save_dir, 'motif_scan_vocab_ame.tsv')
ame_result = pd.read_csv(ame_result_file, sep='\t')
ame_result['motif_DB'] = 'hocomoco_human_trim'
ame_result['motif'] = ame_result.motif_ID.map(tf_id_abbr_to_name_dict)
ame_result = ame_result[~ame_result.motif.isna()]
ame_result['in_vocab'] = ame_result.motif.isin(list_tfs)

print('enriched motifs')
enriched_motifs = sorted(ame_result.motif.unique())
print(len(enriched_motifs))
print(enriched_motifs)

print('enriched motifs with tested vocab pair')
ame_result_filt = ame_result[ame_result.in_vocab]
print(ame_result.shape, ame_result_filt.shape)
enriched_motifs_vocab = sorted(ame_result_filt.motif.unique())
print(len(enriched_motifs_vocab))
print(enriched_motifs_vocab)
# -

motif_df = pd.DataFrame({'motif':list(tf_id_abbr_to_name_dict.values())})
motif_df['len_motif'] = motif_df.motif.map(motif_to_len)
motif_df['in_vocab'] = motif_df.motif.isin(list_tfs)
motif_df['enriched_vocab'] = motif_df.motif.isin(enriched_motifs)
motif_df = motif_df[~motif_df.len_motif.isna()]

motif_df[motif_df.in_vocab & (~motif_df.enriched_vocab)]

pd.read_csv(os.path.join(save_dir, 'mpra_oligo_df_final_SUBMIT_091120.txt'),sep='\t').shape#oligo.str.len().value_counts()



StarcoderdataJupyterScriptsDedupFiltered