<jupyter_start><jupyter_text>05/14/2020 1. reproduce max's P-P, E-P, and E-E 's Maxs paper"Promoter/Enhancer Annotations of HiChIP Loops Promoters were defined as 1 kb regions centered at the TSS, and enhancers were identified as chromHMM enhancers not overlapping with promoters in any cell type. We annotated loop anchors as ‘other’ if the anchor did not contain a promoter or enhancer as defined above."\<jupyter_code>import pandas as pd
import numpy as np
import subprocess
import seaborn as sns
import pybedtools<jupyter_output><empty_output><jupyter_text>1A. make promoter bed<jupyter_code>tss_file = '../data/external/TSS_annon_hg19.csv'
tss_df = pd.read_csv(tss_file,index_col=0)
tss_df['start'] = tss_df['loc'].apply(lambda x: max(0,x-1000))
tss_df['stop'] = tss_df['end'].apply(lambda x: x+1000)
tss_df = tss_df[['chr','start','stop','TSS']]
tss_bed = pybedtools.BedTool.from_dataframe(tss_df)
tss_bed.to_dataframe()[:5]<jupyter_output><empty_output><jupyter_text>1B. get enhancer bedE116 - GM12878E123 - K561<jupyter_code>gm12878_chromHMM_bed = pybedtools.BedTool('../data/external/chromHMM/E116_15_coreMarks_dense.bed')
K562_chromHMM_bed = pybedtools.BedTool('../data/external/chromHMM/E123_15_coreMarks_dense.bed')<jupyter_output><empty_output><jupyter_text>K562_chromHMM_bed.to_dataframe()  1C. get loops/anchors<jupyter_code>gm_loops = pd.read_excel('../data/external/mumbach2017/supp_table_2.xlsx',sheet_name='GM12878 H3K27ac Loops')
k562_loops = pd.read_excel('../data/external/mumbach2017/supp_table_2.xlsx',sheet_name='K562 H3K27ac Loops')
print(gm_loops.shape, k562_loops.shape)
gm_loops[:5]
gm_loops['source'] =  'chr' + gm_loops['chr_x'].map(str) + '_' + gm_loops['start_x'].map(str)+ '_' + gm_loops['stop_x'].map(str)
gm_loops['target'] =  'chr' + gm_loops['chr_y'].map(str) + '_' + gm_loops['start_y'].map(str)+ '_' + gm_loops['stop_y'].map(str)
k562_loops['source'] =  'chr' + k562_loops['chr_x'].map(str) + '_' + k562_loops['start_x'].map(str)+ '_' + k562_loops['stop_x'].map(str)
k562_loops['target'] =  'chr' + k562_loops['chr_y'].map(str) + '_' + k562_loops['start_y'].map(str)+ '_' + k562_loops['stop_y'].map(str)
anchor_1 = gm_loops[['chr_x','start_x','stop_x']]
anchor_2 = gm_loops[['chr_y','start_y','stop_y']]
anchor_1.columns = ['chr','start','stop']
anchor_2.columns = ['chr','start','stop']
gm_anchors = pd.concat([anchor_1,anchor_2],axis=0).drop_duplicates()
gm_anchors['chr'] = 'chr'+gm_anchors['chr'].map(str)
gm_anchors['anchor_name'] = gm_anchors['chr']+ '_' + gm_anchors['start'].map(str)+ '_' + gm_anchors['stop'].map(str)

anchor_1 = k562_loops[['chr_x','start_x','stop_x']]
anchor_2 = k562_loops[['chr_y','start_y','stop_y']]
anchor_1.columns = ['chr','start','stop']
anchor_2.columns = ['chr','start','stop']
k562_anchors = pd.concat([anchor_1,anchor_2],axis=0).drop_duplicates()
k562_anchors['chr'] = 'chr'+k562_anchors['chr'].map(str)
k562_anchors['anchor_name'] = k562_anchors['chr']+ '_' + k562_anchors['start'].map(str)+ '_' + k562_anchors['stop'].map(str)

print(gm_anchors.shape, k562_anchors.shape)
gm_anchors_bed = pybedtools.BedTool.from_dataframe(gm_anchors)
k562_anchors_bed = pybedtools.BedTool.from_dataframe(k562_anchors)
gm_anchors_bed.to_dataframe()<jupyter_output><empty_output><jupyter_text>1Ci for gms<jupyter_code>anchor_bed_prom = gm_anchors_bed.intersect(tss_bed,wa=True,wb=True)
promoter_anchors = anchor_bed_prom.to_dataframe()['name'].unique()
print(promoter_anchors.shape)
anchor_bed_chromHMM = gm_anchors_bed.intersect(gm12878_chromHMM_bed,wa=True,wb=True).to_dataframe(header=None)
anchor_bed_chromHMM_enh = anchor_bed_chromHMM[anchor_bed_chromHMM.loc[:,7].isin(['7_Enh','6_EnhG'])]
enhancer_anchors = anchor_bed_chromHMM_enh.loc[:,3].unique()#.shape 
print(enhancer_anchors.shape) # num enhancers before filter promoters out
enhancer_anchors = np.array(list(set(enhancer_anchors) - set(promoter_anchors)))
print(enhancer_anchors.shape) # num enhancers after filter promoters out

print('num anchors', gm_anchors.shape[0])
print('% anchors captured by P or E: ',(enhancer_anchors.shape[0] +promoter_anchors.shape[0] )/gm_anchors.shape[0] )
print('% anchors captured by P: ',(promoter_anchors.shape[0] )/gm_anchors.shape[0] )
print('% anchors captured by E: ',(enhancer_anchors.shape[0] )/gm_anchors.shape[0] )

loop_df = gm_loops
loop_df['source_P'] = loop_df.source.isin(promoter_anchors)
loop_df['target_P'] = loop_df.target.isin(promoter_anchors)
loop_df['source_E'] = loop_df.source.isin(enhancer_anchors)
loop_df['target_E'] = loop_df.target.isin(enhancer_anchors)
    
loop_df['is_P_P'] = loop_df['source_P'] & loop_df['target_P']
loop_df['is_P_E'] = (loop_df['source_P'] & loop_df['target_E']) |  (loop_df['source_E'] & loop_df['target_P'])
loop_df['is_E_E'] =loop_df['source_E'] & loop_df['target_E']
print('num loops',loop_df.shape[0])
print('num P-P: ',loop_df['is_P_P'].sum())#,loop_df['is_P_P'].sum()/loop_df.shape[0])
print('num P-E: ',loop_df['is_P_E'].sum())#,loop_df['is_P_P'].sum()/loop_df.shape[0])
print('num E-E: ',loop_df['is_E_E'].sum())#,loop_df['is_P_P'].sum()/loop_df.shape[0])
sum_loops = loop_df['is_P_P'].sum() + loop_df['is_P_E'].sum() + loop_df['is_E_E'].sum()
print('num loops to focus on', sum_loops)
print('% P-P: ',loop_df['is_P_P'].sum()/sum_loops)#,loop_df['is_P_P'].sum()/loop_df.shape[0])
print('% P-E: ',loop_df['is_P_E'].sum()/sum_loops)#,loop_df['is_P_P'].sum()/loop_df.shape[0])
print('% E-E: ',loop_df['is_E_E'].sum()/sum_loops)#,loop_df['is_P_P'].sum()/loop_df.shape[0])<jupyter_output>num loops 6396
num P-P:  28
num P-E:  68
num E-E:  43
num loops to focus on 139
% P-P:  0.2014388489208633
% P-E:  0.4892086330935252
% E-E:  0.30935251798561153<jupyter_text>1Cii for K562<jupyter_code>anchor_bed_prom = k562_anchors_bed.intersect(tss_bed,wa=True,wb=True)
promoter_anchors = anchor_bed_prom.to_dataframe()['name'].unique()
print(promoter_anchors.shape)
anchor_bed_chromHMM = k562_anchors_bed.intersect(K562_chromHMM_bed,wa=True,wb=True).to_dataframe(header=None)
anchor_bed_chromHMM_enh = anchor_bed_chromHMM[anchor_bed_chromHMM.loc[:,7].isin(['7_Enh','6_EnhG'])]
enhancer_anchors = anchor_bed_chromHMM_enh.loc[:,3].unique()#.shape 
print(enhancer_anchors.shape) # num enhancers before filter promoters out
enhancer_anchors = np.array(list(set(enhancer_anchors) - set(promoter_anchors)))
print(enhancer_anchors.shape) # num enhancers after filter promoters out

print('num anchors', k562_anchors.shape[0])
print('% anchors captured by P or E: ',(enhancer_anchors.shape[0] +promoter_anchors.shape[0] )/k562_anchors.shape[0] )
print('% anchors captured by P: ',(promoter_anchors.shape[0] )/k562_anchors.shape[0] )
print('% anchors captured by E: ',(enhancer_anchors.shape[0] )/k562_anchors.shape[0] )

loop_df = k562_loops
loop_df['source_P'] = loop_df.source.isin(promoter_anchors)
loop_df['target_P'] = loop_df.target.isin(promoter_anchors)
loop_df['source_E'] = loop_df.source.isin(enhancer_anchors)
loop_df['target_E'] = loop_df.target.isin(enhancer_anchors)
    
loop_df['is_P_P'] = loop_df['source_P'] & loop_df['target_P']
loop_df['is_P_E'] = (loop_df['source_P'] & loop_df['target_E']) |  (loop_df['source_E'] & loop_df['target_P'])
loop_df['is_E_E'] =loop_df['source_E'] & loop_df['target_E']
print('num loops',loop_df.shape[0])
print('num P-P: ',loop_df['is_P_P'].sum())#,loop_df['is_P_P'].sum()/loop_df.shape[0])
print('num P-E: ',loop_df['is_P_E'].sum())#,loop_df['is_P_P'].sum()/loop_df.shape[0])
print('num E-E: ',loop_df['is_E_E'].sum())#,loop_df['is_P_P'].sum()/loop_df.shape[0])
sum_loops = loop_df['is_P_P'].sum() + loop_df['is_P_E'].sum() + loop_df['is_E_E'].sum()
print('num loops to focus on', sum_loops)
print('% P-P: ',loop_df['is_P_P'].sum()/sum_loops)#,loop_df['is_P_P'].sum()/loop_df.shape[0])
print('% P-E: ',loop_df['is_P_E'].sum()/sum_loops)#,loop_df['is_P_P'].sum()/loop_df.shape[0])
print('% E-E: ',loop_df['is_E_E'].sum()/sum_loops)#,loop_df['is_P_P'].sum()/loop_df.shape[0])<jupyter_output>num loops 1553
num P-P:  118
num P-E:  459
num E-E:  394
num loops to focus on 971
% P-P:  0.121524201853759
% P-E:  0.4727085478887745
% E-E:  0.4057672502574665<jupyter_text>```kk i compared to max's results9:21this is max's GM12878 data percentage:num anchors 10714% anchors captured by P or E:  0.21849915997759942% anchors captured by P:  0.08418891170431211% anchors captured by E:  0.13431024827328727num loops 6396num P-P:  338num P-E:  1263num E-E:  2136num loops to focus on 3737% P-P:  0.09044688252609044% P-E:  0.337971635001338% E-E:  0.57158148247257169:21i modified our stuff9:22P-P (don't care about accessibility here9:22P-E = P-PIRacc (so a promoter region connecting to an accessible region that is not a promoter)9:22E-E = PIRacc-PIRacc9:22and then you get these percentages9:22perc_is_P_P_adj 0.08977712805519222perc_is_P_PIRacc_adj 0.3681875304557462perc_is_PIRacc_PIRacc_adj 0.54203534148906159:27this is max's K562 datanum anchors 2955% anchors captured by P or E:  0.7922165820642978% anchors captured by P:  0.30524534686971233% anchors captured by E:  0.48697123519458546num loops 1553num P-P:  118num P-E:  459num E-E:  394num loops to focus on 971% P-P:  0.121524201853759% P-E:  0.4727085478887745% E-E:  0.4057672502574665```  2. compare our GM loops/anchors to max's loops/anchors 2A anchors<jupyter_code>tissue='GM12878'
loop_file = '../data/interim/merged/loops/'+tissue+'.loops.csv'
loop_df = pd.read_csv(loop_file, index_col=0)
anchor_file = '../data/interim/merged/anchors/'+tissue+'.anchors.csv'
anchor_df = pd.read_csv(anchor_file, index_col=0)
anchor_df = anchor_df[['chr','start','end','anchors','sample']]#[:5]
anchor_bed = pybedtools.BedTool.from_dataframe(anchor_df)

# anchor_annon_file = os.path.join(anchor_annon_dir, tissue+'_annon.bed')
# anchor_annon_df = pd.read_csv(anchor_annon_file, names=['chr', 'start', 'stop', 'name', 'atac_chr', 'atac_start', 'atac_stop'], sep='\t')
# anchor_to_count_dict = anchor_annon_df.groupby('name').name.count().to_dict()
# anchors_acc = anchor_to_count_dict.keys()

# print(anchor_df.shape)

# promoter_anchor_file = os.path.join(promoter_annon_dir, 'promoter_'+'GM12878'+'_annon.bed')
# promoter_anchor_df = pd.read_csv(promoter_anchor_file, names=['chr', 'start', 'stop', 'tss', 'anchor_chr', 'anchor_start', 
#                                                          'anchor_stop', 'anchor', 'overlap'], sep='\t')
# promoter_anchors = promoter_anchor_df.anchor.unique()
# print(promoter_anchors.shape)
# union_anchors_df = pd.concat([gm_anchors_bed.to_dataframe(),anchor_bed.to_dataframe()])
# union_anchors_df = union_anchors_df[['chrom','start','end']]
shared_anchors = (gm_anchors_bed+anchor_bed)#.sort().merge()
# union_bed = pybedtools.BedTool.from_dataframe(union_anchors_df).sort().merge()
print('# anchors we have: ',anchor_bed.count()) 
print('# anchors Max has: ', gm_anchors_bed.count())
print('# shared anchors: ', shared_anchors.count())
print('# anchors Max had we didnt: ',(gm_anchors_bed-anchor_bed).count())#.merge().count())
print('# anchors we had Max didnt: ',(anchor_bed-gm_anchors_bed).count())#.merge().count())
print('# total union set of anchors: ', union_bed.count())

shared_anchors.count()/gm_anchors_bed.count()
anchor_bed.count()/gm_anchors_bed.count()<jupyter_output><empty_output><jupyter_text>we capture 81.5% of the anchors (sensitivitity) that Max predicted, we also have almost 8x as many anchor points  2B loops<jupyter_code>def name_loop_row(row):
    
    if row.x_before_y:
        return '::'.join([row.source,row.target])
    else:
        return '::'.join([row.target,row.source])

gm_loops = gm_loops[gm_loops.chr_x ==gm_loops.chr_y]
gm_loops['x_before_y'] = gm_loops.start_y > gm_loops.start_x
gm_loops['loop_name'] = gm_loops.apply(name_loop_row,axis=1)

loop_df[['chr_x','start_x','stop_x']] = loop_df.source.str.split('_',expand=True)
loop_df[['chr_y','start_y','stop_y']] = loop_df.target.str.split('_',expand=True)
loop_df = loop_df[loop_df.chr_x ==loop_df.chr_y]
loop_df['x_before_y'] = loop_df.start_y.map(int) > loop_df.start_x.map(int)
loop_df['loop_name'] = loop_df.apply(name_loop_row,axis=1)
our_loops_set = set(loop_df.loop_name)
max_loops_set = set(gm_loops.loop_name)

print('num our loops: ', len(our_loops_set))
print('num max loops: ', len(max_loops_set))
print('shared loops: ', len(our_loops_set.intersection(max_loops_set)))
len(our_loops_set.intersection(max_loops_set))/len(max_loops_set)
len(our_loops_set)/len(max_loops_set)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17427952155414625, 'fraction_numerical': 0.05836076101790158, 'mean_word_length': 7.627423822714682, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '59978', 'license': 'MIT', 'n_tokens_mistral': 5680, 'n_tokens_neox': 5165, 'n_words': 726, 'path': 'notebooks/4B_fig2_hichip_stats_mumbach.ipynb', 'repo_name': 'mguo123/pan_omics', 'size': '12457'}