<jupyter_start><jupyter_text>HiC-analysis - use only cell type-specific pairsTodo:- assign 900 clusters into 8 cell types? (maybe not)- or how do we figure out the cluster assignment for each linked pair?<jupyter_code># enhancers lifted
min_match = 0.1

f = '/cndd2/fangming/projects/scf_enhancers/enhancer_metacells_r100_n935/liftover/enhs_list_reformat_{:.2f}_lifted.bed'.format(min_match)
enh_list = pd.read_csv(f, sep='\t', header=None, names=['chr', 'start', 'end', 'enh_id'])
enh_list['center'] = (enh_list['start'] + enh_list['end'])/2
print(enh_list.shape)
enh_list.head()

# tss lifted
# f = '/cndd2/fangming/projects/scf_enhancers/enhancer_metacells_r100_n935/liftover/gencode.vM16.annotation_genes_lifted.bed'
f = '/cndd2/fangming/projects/scf_enhancers/enhancer_metacells_r100_n935/liftover/gencode.vM16.annotation_genes_tss_lifted.bed'
tss_list = pd.read_csv(f, sep='\t', header=None, names=['chr', 'start', 'end', 'gid', 'gene_name'])
tss_list['gid'] = tss_list['gid'].apply(lambda x: x.split('.')[0])
print(tss_list.shape)
tss_list.head()

# enh - tss linked list
# orig significant pairs (enh_index and gid)
f = '/cndd2/fangming/projects/scf_enhancers/enhancer_metacells_r100_n935/results_new/fig2_metacell_201213_linked_pairs_mc.tsv'
linked_pairs_mc = pd.read_csv(f, sep="\t")
print(linked_pairs_mc.shape)

f = '/cndd2/fangming/projects/scf_enhancers/enhancer_metacells_r100_n935/results_new/fig2_metacell_201213_linked_pairs_atac.tsv'
linked_pairs_atac = pd.read_csv(f, sep="\t")
print(linked_pairs_atac.shape)

f = '/cndd2/fangming/projects/scf_enhancers/enhancer_metacells_r100_n935/results_new/fig2_metacell_201213_correlated_pairs_mc.tsv'
correlated_pairs_mc = pd.read_csv(f, sep="\t")
print(correlated_pairs_mc.shape)

f = '/cndd2/fangming/projects/scf_enhancers/enhancer_metacells_r100_n935/results_new/fig2_metacell_201213_correlated_pairs_atac.tsv'
correlated_pairs_atac = pd.read_csv(f, sep="\t")
print(correlated_pairs_atac.shape)

linked_pairs_mc.head()
linked_pairs_both = pd.merge(linked_pairs_mc, linked_pairs_atac, on=['gene', 'enh', 'dist'])
correlated_pairs_both = pd.merge(correlated_pairs_mc, correlated_pairs_atac, on=['gene', 'enh', 'dist'])
print(linked_pairs_both.shape, 
      correlated_pairs_both.shape, 
     )
# gene expression profiles
f = '/sphere/fangming/enhancers/scripts/data_organized/gene_rna_clsts.tsv'
gene_rna_clsts = pd.read_csv(f, sep='\t', index_col=0) #header=True, index=True)
print(gene_rna_clsts.shape)
gene_rna_clsts.head()

# cluster annotation
f = '/cndd2/fangming/projects/scf_enhancers/enhancer_ethan38_200520/ref/annotations_order.tsv'
clst_annot = pd.read_csv(f, sep='\t').set_index('cluster')['annotation']
print(clst_annot.shape)
clst_annot.head()
# merge cluster name
annot_rename_lookup = {
    'L2/3': "L2_3",
    'L4/5': "L4",
    'L6b': "L6",
    'Sncg': "Ndnf",
    'Lamp5': "Ndnf",
}

def rename_annot(x, annot_rename_lookup):
    """
    """
    if x in annot_rename_lookup.keys():
        return annot_rename_lookup[x]
    else:
        return x
    
clst_annot_broad = pd.Series([annot.split(' ')[0] for annot in clst_annot.values], 
                             index=clst_annot.index)
clst_annot_broad = clst_annot_broad.apply(lambda x: rename_annot(x, annot_rename_lookup))

print(clst_annot_broad.shape, 
      len(np.unique(clst_annot_broad.values)),
      np.unique(clst_annot_broad.values),
     )
# merge RNA profiles
_tmp = gene_rna_clsts.T
_tmp['annot'] = clst_annot_broad.loc[_tmp.index]
gene_rna_clsts_broad = _tmp.groupby('annot').median().T
gene_rna_clsts_broad_rank = gene_rna_clsts_broad.rank(pct=True, axis=1)
print(gene_rna_clsts_broad_rank.shape)
gene_rna_clsts_broad_rank.head()
active_th = 0.5
# distances
distance_cutoff = 1e5
distance_idx = np.arange(1, int(distance_cutoff/resolution)+1)
distances = distance_idx * resolution

for celltype in celltypes:
    print(celltype)
    
    contacts = [[] for i in np.arange(len(distances))]
    contacts_sig_pairs = collections.OrderedDict({
            'linked_mc': [[] for i in np.arange(len(distances))],
            'linked_atac': [[] for i in np.arange(len(distances))],
            'linked_both': [[] for i in np.arange(len(distances))],
            'correlated_mc': [[] for i in np.arange(len(distances))],
            'correlated_atac': [[] for i in np.arange(len(distances))],
            'correlated_both': [[] for i in np.arange(len(distances))],
        })

    for chromosome in chromosomes:
        print(celltype, resolution, chromosome)

        # read in the matrix 
        resource_string = (
             '/sphere/fangming/enhancers/data_m3c/{}_all_brain.txt_1kb_contacts.mcool'
             '::/resolutions/{}'.format(celltype, resolution)
             )
        c = cooler.Cooler(resource_string)
        mat_all = c.matrix(balance=True).fetch(chromosome)

        # all enhancers and tss on that chromosomes
        enhs = enh_list[enh_list['chr']==chromosome]
        tsss = tss_list[tss_list['chr']==chromosome]
        
        # only include the gene if it is actively expressed (>50% across cell types)
        active_genes = gene_rna_clsts_broad_rank[gene_rna_clsts_broad_rank[celltype]>active_th].index.values
        tsss = tsss[tsss['gid'].isin(active_genes)]
        
        # get paired bins
        paired_bin_tables = collections.OrderedDict({
            'linked_mc': enh_gene_id_to_binidx(linked_pairs_mc, tsss, enhs, resolution),
            'linked_atac': enh_gene_id_to_binidx(linked_pairs_atac, tsss, enhs, resolution),
            'linked_both': enh_gene_id_to_binidx(linked_pairs_both, tsss, enhs, resolution),
            'correlated_mc': enh_gene_id_to_binidx(correlated_pairs_mc, tsss, enhs, resolution),
            'correlated_atac': enh_gene_id_to_binidx(correlated_pairs_atac, tsss, enhs, resolution),
            'correlated_both': enh_gene_id_to_binidx(correlated_pairs_both, tsss, enhs, resolution),
        })
        track_names = list(paired_bin_tables.keys())

        # contact stats: distance bin by corr values in that bin
        (_contacts, _contacts_sig_pairs
        ) = get_contacts(mat_all, paired_bin_tables, resolution, distance_idx)

        # combine all chromosomes 
        for idx, dist in enumerate(distances):
            contacts[idx] += list(_contacts[idx])
            for key in _contacts_sig_pairs.keys():
                contacts_sig_pairs[key][idx] += list(_contacts_sig_pairs[key][idx])
#         break

    # stats
    contacts_mean = np.array([np.nanmean(val) for val in contacts])
    contacts_std = np.array([np.nanstd(val) for val in contacts])
    contacts_n = np.array([len(val) for val in contacts])

    contacts_sig_pairs_mean = collections.OrderedDict({
        key: np.zeros(len(distances)) 
        for key in contacts_sig_pairs.keys()
    })
    contacts_sig_pairs_std = contacts_sig_pairs_mean.copy()
    contacts_sig_pairs_n = contacts_sig_pairs_mean.copy()

    for key in contacts_sig_pairs.keys():
        contacts_sig_pairs_mean[key] = np.array([np.nanmean(val) 
                                                 for val in contacts_sig_pairs[key]])

        contacts_sig_pairs_std[key] = np.array([np.nanstd(val) 
                                                for val in contacts_sig_pairs[key]])
        
        contacts_sig_pairs_n[key] = np.array([len(val) 
                                              for val in contacts_sig_pairs[key]])

    
    # run and save results
    output_res = output_res_format.format(int(100*active_th), celltype, resolution, today)
    print(output_res)

    # save results
    # distances, contacts_mean, contacts_std ## these are numpy arrays 
    # contacts_sig_pairs_mean, contacts_sig_pairs_std ## these
    with open(output_res, 'wb') as fwh:
        pickle.dump((distances, 
                     contacts_mean, contacts_std, contacts_n,
                     contacts_sig_pairs_mean, 
                     contacts_sig_pairs_std,
                     contacts_sig_pairs_n,
                    ), fwh)
#     break<jupyter_output>L2_3
L2_3 10000 chr1
{'fraction_non_alphanumeric': 0.11816600397614314, 'fraction_numerical': 0.01677435387673956, 'mean_word_length': 3.6850989522700814, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '40316', 'license': 'MIT', 'n_tokens_mistral': 3038, 'n_tokens_neox': 2907, 'n_words': 457, 'path': 'archives/hic/contact_freq_enrichment-metacell-Dec13-bck.ipynb', 'repo_name': 'mukamel-lab/SingleCellFusion_EnhancerPaper', 'size': '8048'}