<jupyter_start><jupyter_text>LB1 in KDDs<jupyter_code>density_dat = pd.read_csv('tables_out/all_dat_tog.tsv', sep='\t').replace('EarlySomite','Early somite')

density_dat['score_lb1'] = density_dat[['score0_lb1','score1_lb1']].mean(axis=1)
density_dat['score_k9'] = density_dat[['score0_k9','score1_k9']].mean(axis=1)

f, axB = plt.subplots(figsize=(6,8))

b = sns.violinplot(x='score_lb1', y='cell_type', data=density_dat, kind='violin',
           hue='category_k9',linewidth=0.5, inner=None, ax=axB,
                   hue_order=['T1-KDD','T2-KDD','nonKDD'],
          order=['ESCs','Mid-hindgut','Liver','Anterior primitive streak',
                'Paraxial mesoderm','Early somite','Artery progenitors',d
                'Endothelial progenitors','Cardiac myocytes','Epicardium',
                'Definitive ectoderm','Border ectoderm','Midbrain progenitors'],
                  palette = mypal)
plt.legend(frameon=False, title='KDD category', bbox_to_anchor=(1, 1.15),
          ncol=1)
# plt.xticks(rotation=45, ha='right')
plt.gca().spines['top'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
axB.set(ylabel='cell type', xlabel='LB1 binding [log2(read counts)]')
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Proportion expressed genes in LADs A and B are in .ai file already (track view and violins)C: gene expressionD: gene density / ATAC-seq / B compartment enrichments<jupyter_code>gtf_path = 'hg38/kallisto_indices/homo_sapiens/'
lads_path = 'LAD_BEDs/'
lad_dat = 'LADs/dat'
path_to_gene_results = 'LADs/dat/genes'
path_to_atac = 'LADs/dat/ATACseq'

## gene expression

rnaseq_dat = pd.read_table('RNAseq_quantile_normalized_dat.tsv',
                          header=0, names=['name','tpm_CardiacMyocytes','tpm_H9ESC',
                                           'tpm_EarlySomite','tpm_ParaxMesoderm',
                                           'tpm_MidHindgut','tpm_EndoProgenitor',
                                           'tpm_DefEctoderm'])




# get locations of genes

# get LAD annotations per gene

# Kallisto analysis used Ensembl gene annotations v96, should do that here as well for continuity

ensembl = pd.read_csv(os.path.join(gtf_path,'Homo_sapiens.GRCh38.96.gtf'),
                     sep='\t', skiprows=5, header=None, 
                      names=['chrom','source','desc','start','end','score','strand','score2','desc_further'],
                     low_memory=False)
ensembl['chr'] = 'chr' + ensembl['chrom'].astype(str)
ensembl['name'] = ensembl['desc_further'].str.split('gene_name ').str[1].str.split(';').str[0].str.strip('""')

# filter for genes only

ensembl_genes = ensembl.query('desc == "gene"').copy()
n_ens_genes = pd.DataFrame(ensembl_genes['name'].value_counts())
n_ens_genes['gene'] = n_ens_genes.index
ensembl_genes_single = n_ens_genes.query('name == 1')['gene'].tolist()
ensembl_genes_single_df = ensembl_genes.query('(name in @ensembl_genes_single)').copy()

# make BED file

ens_single_bed = BedTool.from_dataframe(ensembl_genes_single_df[['chr','start','end','name']].copy()).sort()



def get_lad_state(row):
    if row['T1-LAD']:
        return('T1-LAD')
    elif row['T2-LAD']:
        return('T2-LAD')
    else:
        return('nonLAD')

dfs = []

for val in rnaseq_dat.columns[1:]:
    ct = val.split('_')[1]
    t1lads = BedTool(os.path.join(lads_path, f'{ct}_T1-LADs.bed'))
    genes_in_t1lads = ens_single_bed.intersect(t1lads, f=1.0).to_dataframe()['name'].drop_duplicates().tolist()
    t2lads = BedTool(os.path.join(lads_path, f'{ct}_T2-LADs.bed'))
    genes_in_t2lads = ens_single_bed.intersect(t2lads, f=1.0).to_dataframe()['name'].drop_duplicates().tolist()
    df = rnaseq_dat[['name', f'tpm_{ct}']].copy()
    df.columns = ['name','expression']
    df['T1-LAD'] = df['name'].isin(genes_in_t1lads)
    df['T2-LAD'] = df['name'].isin(genes_in_t2lads)
    df['LAD_state'] = df.apply(lambda row: get_lad_state(row), axis=1)
    df['celltype'] = ct
    dfs.append(df)
    
tog = pd.concat(dfs, sort=False)

# look at proportion of expressed genes

prop_exp = []
cats = []
cts = []

for ct in tog['celltype'].unique():
    for cat in tog['LAD_state'].unique():
        prop_exp.append(len(tog.query('celltype == @ct and LAD_state == @cat and expression > 5')) / len(tog.query('celltype == @ct and LAD_state == @cat')))
        cats.append(cat)
        cts.append(ct)
        

df = pd.DataFrame({
    'proportion_expressed':prop_exp,
    'LAD_cat':cats,
    'celltype':cts
})

mypal = {'T1-LAD':'indigo','T2-LAD':'mediumpurple','nonLAD':'lightgrey'}

f, ax = plt.subplots(figsize=(4,3))
sns.barplot(x='LAD_cat', y='proportion_expressed', data=df,
           order=['T1-LAD','T2-LAD','nonLAD'], ci='sd',
           palette = mypal)<jupyter_output><empty_output><jupyter_text>Gene, ATAC-peaks, B compartment enrichments<jupyter_code>## feature enrichments

def get_enrichment_df(enr_results_path, enr_feature):
    """
    Collects enrichment results (pre-computed) into a dataframe for plotting
    """
    zscores = []
    pvalues = []
    lad_cts = []
    enr_elements = []
    ladtypes = []

    for ct in cts:

        enr_dat_t1lad = pd.read_csv(f'{enr_results_path}/{enr_feature}_{ct}_T1LADs_enrichment_results.tsv',
                            sep='\t').replace(0.0, 0.00001)
        enr_dat_t2lad = pd.read_csv(f'{enr_results_path}/{enr_feature}_{ct}_T2LADs_enrichment_results.tsv',
                            sep='\t').replace(0.0, 0.00001)

        with open(f'{enr_results_path}/{enr_feature}_{ct}_T1LADs_actual_intersections.json', 'r') as f:
            actual_ints_t1lads = json.load(f)

        if actual_ints_t1lads == 0:
            actual_ints_t1lads = 0.00001

        with open(f'{enr_results_path}/{enr_feature}_{ct}_T2LADs_actual_intersections.json', 'r') as f:
            actual_ints_t2lads = json.load(f)

        if actual_ints_t2lads == 0:
            actual_ints_t2lads = 0.00001

        # this loop was originally developed for TEs, but now 
        # works with all the features we assessed for enrichment
        for te_class in enr_dat_t1lad['TE_class'].unique():

            # calculate z-score for T1-LAD

            t1lad_sd = np.std(enr_dat_t1lad.query('TE_class == @te_class')['n_intersections'])
            t1lad_ints = actual_ints_t1lads[te_class]
            t1lad_obs_z = (t1lad_ints - \
                         np.mean(enr_dat_t1lad.query('TE_class == @te_class')['n_intersections'])) / t1lad_sd

            # calculate empirical p-value for T1-LAD

            if t1lad_obs_z > 0.00001:
                pt1lad = len(enr_dat_t1lad.query('(TE_class == @te_class) and (n_intersections > @t1lad_ints)')) /\
                len(enr_dat_t1lad.query('(TE_class == @te_class)'))
            else:
                pt1lad = len(enr_dat_t1lad.query('(TE_class == @te_class) and (n_intersections < @t1lad_ints)')) /\
                len(enr_dat_t1lad.query('(TE_class == @te_class)'))

            zscores.append(t1lad_obs_z)
            pvalues.append(pt1lad)
            lad_cts.append(ct)
            enr_elements.append(te_class)
            ladtypes.append('T1-LAD')

            # calculate z-score for T2-LAD

            t2lad_sd = np.std(enr_dat_t2lad.query('TE_class == @te_class')['n_intersections'])
            t2lad_ints = actual_ints_t2lads[te_class]
            t2lad_obs_z = (t2lad_ints - \
                         np.mean(enr_dat_t2lad.query('TE_class == @te_class')['n_intersections'])) / t2lad_sd

            # calculate empirical p-value for T2-LAD

            if t2lad_obs_z > 0.00001:
                pt2lad = len(enr_dat_t2lad.query('(TE_class == @te_class) and (n_intersections > @t2lad_ints)')) /\
                len(enr_dat_t2lad.query('(TE_class == @te_class)'))
            else:
                pt2lad = len(enr_dat_t2lad.query('(TE_class == @te_class) and (n_intersections < @t2lad_ints)')) /\
                len(enr_dat_t2lad.query('(TE_class == @te_class)'))

            zscores.append(t2lad_obs_z)
            pvalues.append(pt2lad)
            lad_cts.append(ct)
            enr_elements.append(te_class)
            ladtypes.append('T2-LAD')     

    zpdf = pd.DataFrame({
        'p_value':pvalues,
        'zscore':zscores,
        'LAD_type':ladtypes,
        'class':enr_elements,
        'celltype':lad_cts
    }).dropna().query('zscore != "-inf"').query('zscore != "inf"')

    zpdf['-log10p'] = -np.log10(zpdf['p_value'])
    zpdf['-log10p'] = zpdf['-log10p'].replace(np.inf, 1.0)
    return zpdf
        
def make_enrichment_df(feature_name, lad_or_kdd, cts):
    # df for merged plot

    zscores = []
    pvalues = []
    domain_cts = []
    classes = []
    domaintypes = []

    for ct in cts:

        t1dat = pd.read_csv(f'{lad_dat}/{feature_name}/{feature_name}_{ct}_T1{lad_or_kdd}s_enrichment_results.tsv',
                            sep='\t').replace(0.0, 0.00001)
        t2dat = pd.read_csv(f'{lad_dat}/{feature_name}/{feature_name}_{ct}_T2{lad_or_kdd}s_enrichment_results.tsv',
                            sep='\t').replace(0.0, 0.00001)
        
        class_dats = {
            'T1':t1dat,
            'T2':t2dat
        }

        with open(f'{lad_dat}/{feature_name}/{feature_name}_{ct}_T1{lad_or_kdd}s_actual_intersections.json', 'r') as f:
            actual_ints_t1 = json.load(f)

        if actual_ints_t1 == 0:
            actual_ints_t1 = 0.00001

        with open(f'{lad_dat}/{feature_name}/{feature_name}_{ct}_T2{lad_or_kdd}s_actual_intersections.json', 'r') as f:
            actual_ints_t2 = json.load(f)

        if actual_ints_t2 == 0:
            actual_ints_t2 = 0.00001

        for dat_class in t1dat['feature_class'].unique():
            
            for class_type in ['T1', 'T2']:
                
                dat = class_dats[class_type]

                # calculate z-score

                sd = np.std(dat.query('feature_class == @dat_class')['n_intersections'])
                if class_type == 'T1':
                    ints = actual_ints_t1[dat_class]
                else:
                    ints = actual_ints_t2[dat_class]
                obs_z = (ints - \
                         np.mean(dat.query('feature_class == @dat_class')['n_intersections'])) / sd

                # calculate empirical p-value
                
                if obs_z > 0.00001:
                    pval = len(dat.query('(feature_class == @dat_class) and (n_intersections > @ints)')) /\
                    len(t1dat.query('(feature_class == @dat_class)'))
                else:
                    pval = len(dat.query('(feature_class == @dat_class) and (n_intersections < @ints)')) /\
                    len(dat.query('(feature_class == @dat_class)'))

                zscores.append(obs_z)
                pvalues.append(pval)
                domain_cts.append(ct)
                classes.append(dat_class)
                domaintypes.append(f'{class_type}-{lad_or_kdd}')

    zpdf = pd.DataFrame({
        'p_value':pvalues,
        'zscore':zscores,
        f'{lad_or_kdd}_type':domaintypes,
        'class':classes,
        'celltype':domain_cts
    }).dropna().query('zscore != "-inf"').query('zscore != "inf"')

    zpdf['p_value'] = zpdf['p_value'].replace(0.0, 0.0001)
    zpdf['-log10p'] = -np.log10(zpdf['p_value'])
    zpdf['-log10p'] = zpdf['-log10p'].replace(np.inf, 1.0)  
    return(zpdf)
        
        
        

cts = ['Liver', 'CardiacMyocytes', 'EndoProgenitor', 'D5Midbrain',
       'H9ESC', 'Epicardium', 'MidHindgut', 'ParaxMesoderm',
       'EarlySomite', 'DefEctoderm', 'BorderEctoderm', 'D4Artery', 'APS']

# enrichment df for genes
zpdf_genes = get_enrichment_df(path_to_gene_results, 'genes')
zpdf_genes['TE_class'] = 'Genes'

# ATAC

def load_atac(f):
    return(BedTool.from_dataframe(pd.read_table(f'{path_to_atac}/GSE85330/{f}',
                                                            header=None, names=['chrom','start','end','name',
                                                                               'score','strand','score2',
                                                                               'score3','score4','score5'])[['chrom','start','end']]).sort().liftover('/pollard/data/functional_genomics/liftOver/hg19ToHg38.over.chain.gz').sort())

h9esc_atac_peaks_rep1 = load_atac('GSM2264826_H9_0_1.filterBL.bed.gz')
h9esc_atac_peaks_rep2 = load_atac('GSM2264827_H9_0_2.filterBL.bed.gz')

h9esc_atac_peaks = h9esc_atac_peaks_rep1.intersect(h9esc_atac_peaks_rep2).sort().merge()

mes_atac_peaks_rep1 = load_atac('GSM2264828_H9_2_1.filterBL.bed.gz')
mes_atac_peaks_rep2 = load_atac('GSM2264829_H9_2_2.filterBL.bed.gz')

mes_atac_peaks = mes_atac_peaks_rep1.intersect(mes_atac_peaks_rep2).sort().merge()

cm_atac_peaks_rep1 = load_atac('GSM2264832_H9_30_1.filterBL.bed.gz')
cm_atac_peaks_rep2 = load_atac('GSM2264833_H9_30_2.filterBL.bed.gz')

cm_atac_peaks = cm_atac_peaks_rep1.intersect(cm_atac_peaks_rep2).sort().merge()

atac_peaks = {}

atac_peaks['H9ESC'] = h9esc_atac_peaks
atac_peaks['ParaxMesoderm'] = mes_atac_peaks
atac_peaks['CardiacMyocytes'] = cm_atac_peaks

atac_kdd_zpdf = make_enrichment_df('ATAC','KDD',atac_peaks.keys())
atac_lad_zpdf = make_enrichment_df('ATAC','LAD',atac_peaks.keys())

bcomp_cts = ['H9ESC','CardiacMyocytes']

bcomp_kdd_zpdf = make_enrichment_df('bcomp', 'KDD', bcomp_cts)
bcomp_lad_zpdf = make_enrichment_df('bcomp', 'LAD', bcomp_cts)

cts = ['Liver', 'CardiacMyocytes', 'EndoProgenitor', 'D5Midbrain',
       'H9ESC', 'Epicardium', 'MidHindgut', 'ParaxMesoderm',
       'EarlySomite', 'DefEctoderm', 'BorderEctoderm', 'D4Artery', 'APS']

zpdf_genes['class'] = 'Genes'

zpdf = pd.concat([zpdf_genes, atac_lad_zpdf, bcomp_lad_zpdf], sort=False)

ct_replace = {
    'CardiacMyocytes':'Cardiac\nmyocytes',
    'EarlySomite':'Early somite',
    'H9ESC':'ESCs',
    'ParaxMesoderm':'Paraxial\nmesoderm',
    'DefEctoderm':'Definitive\nectoderm',
    'MidHindgut':'Mid-hindgut',
    'D5Midbrain':'Midbrain',
    'Liver':'Liver progenitors',
    'BorderEctoderm':'Border ectoderm',
    'EndoProgenitor':'Endothelial\nprogenitors',
    'D4Artery':'Artery progenitors',
    'APS':'Anterior primitive\nstreak'
}

zpdf_genes['celltype'].unique()

zpdf['celltype'].unique()

f, ax = plt.subplots(figsize=(4,3))
sns.scatterplot(y='class', x='zscore', size='-log10p', data=zpdf.replace(ct_replace), hue='celltype',
               style='LAD_type', palette='colorblind', ax=ax, sizes=(20, 200), alpha=0.6, linewidth=0)
plt.legend(ncol=2, bbox_to_anchor=(1, 1), frameon=False)
plt.axvline(0, linestyle='--', color='grey')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15446428571428572, 'fraction_numerical': 0.0260989010989011, 'mean_word_length': 2.9827680525164113, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '591201', 'license': 'MIT', 'n_tokens_mistral': 5828, 'n_tokens_neox': 5401, 'n_words': 702, 'path': 'Figure2_panels.ipynb', 'repo_name': 'ketringjoni/LAD_atlas', 'size': '14560'}