<jupyter_start><jupyter_text>Intersect loops with SNPs and Genes<jupyter_code># load the colocalization data
coloc_fn = 'results/main/2021_Nikhil_eQTL/Results/Colocalization/T1D_34012112_Gaulton/'
coloc_fn += 'BLUEPRINT_eQTL_Monocyte/FINAL_Summary_Coloc_Gene_SNP_Pairs.bed'
coloc = pd.read_table(coloc_fn)
# extract the most significant according the H4 
coloc_sig = coloc[coloc['pp_H4_Coloc_Summary'] > 0.75]
coloc_sig = coloc_sig[['chr', 'pos', 'rs_id', 'variant_id']]
coloc_sig.rename(columns={'pos': 'end'}, inplace=True)
coloc_sig.loc[:, 'start'] = coloc_sig.loc[:, 'end'] - 1
coloc_sig = coloc_sig[['chr', 'start', 'end', 'rs_id', 'variant_id']]
coloc_sig = pbt.BedTool.from_dataframe(coloc_sig)
print('In total there are {} colocalized SNPs.'.format(len(coloc_sig)))
# load the loop data
loop_fn = 'results/main/2021_Nikhil_eQTL/Data/FitHiChIP_Loops/CD4N/FitHiChIP_S/FitHiChIP.interactions_FitHiC_Q0.01.bed'
loops = pd.read_table(loop_fn)
tmp_loops = loops[['chr1', 's1', 'e1', 'chr2', 's2', 'e2']]
tmp_loops.rename(columns={'p': 'score'}, inplace=True)
tmp_loops.loc[:, 'name'] = '.'
tmp_loops.loc[:, 'score'] = loops['p']
tmp_loops.loc[:, 'strand1'] = '.'
tmp_loops.loc[:, 'strand2'] = '.'
loops = pbt.BedTool.from_dataframe(tmp_loops)
print('FitHiChIP found {} significant loops.'.format(loops.count()))
# extract loops where an anchor overlaps with a SNP
snp_loops = loops.pairtobed(coloc_sig)
snp_loops_tmp = snp_loops.to_dataframe(disable_auto_names=True, header=None)
nsnp_loops = sum(snp_loops_tmp.duplicated(subset=[0,1,2,3,4,5]))
print('There are {} loops with colocalized SNPs.'.format(nsnp_loops))
# load the gene coords
genes_fn = 'results/refs/ensembl/Homo_sapiens.GRCh38.104.bed'
genes = pbt.BedTool(genes_fn)
# extract loops where an anchor overlaps with a SNP
gs_loops = snp_loops.pairtobed(genes)
gs_loops = gs_loops.to_dataframe(disable_auto_names=True, header=None)
gs_loops
gs_loops.columns = bedpe_10cols + ['chrS', 'startS', 'endS', 'rs_id', 'variant_id'] + \
                                  ['chrG', 'startG', 'endG', 'gene_name', 'gene_num', 'gene_id']
gs_loops<jupyter_output><empty_output><jupyter_text>Summary of loop intersections   Number of loops per SNP<jupyter_code>gs_loops.rs_id.value_counts()<jupyter_output><empty_output><jupyter_text>Number of loops per gene<jupyter_code>gs_loops.gene_name.value_counts()<jupyter_output><empty_output><jupyter_text>Convert all results to output and WashU visualization results<jupyter_code>def bedpe_to_WashU_longrange(fn, df):
    """
        Convert from a loop bedpe file into WashU longrange, 
        includes bgzip and tabix of the fn. 
        
        Params
        -------
        fn: str
            path to the longrange output file (without gz)
            
        df: dataframe
            columns 1-6 are as expected and column 8 is the p or q-value. 
            
        Output
        ------
        gzfn: str
            path to the longrange with bgzip compression
        tabix_fn: str
            path to the index of the longrange file
            
    """

    # parsing the data into WashU longrage format
    data = []
    for sr in df.values.tolist():

        # calculate the -log(FDR)
        qval = -np.log(sr[7])

        # get the first pair data
        second_pair_str = '{}:{}-{},{:.5f}'.format(*sr[3:6], qval)
        first_row = sr[0:3] + [second_pair_str]

        # get the second pair data
        first_pair_str = '{}:{}-{},{:.5f}'.format(*sr[0:3], qval)
        second_row = sr[3:6] + [first_pair_str]

        # add each data row
        data.append(first_row)
        data.append(second_row)

    data = sorted(data, key=lambda x: (x[0], x[1], x[2]))

    # writing out the data
    with open(fn, 'w') as f:
        for line in data:
            info = [str(x) for x in line]
            info = '\t'.join(info)
            f.write(info + '\n')
            
    # run bgzip
    cmd = '{} {}'.format(bgzip, lrange_fn)
    print(cmd)
    job = sp.Popen(cmd, stderr=sp.PIPE,stdout=sp.PIPE, shell=True)

    out, err = job.communicate()
    print('out:', out.decode())
    print('err:', err.decode())
    
    # run tabix
    lrange_gzfn = os.path.join(outdir, 'gs_loops.longrange.bed.gz')
    cmd = '{} {}'.format(tabix, lrange_gzfn)
    print(cmd)
    job = sp.Popen(cmd, stderr=sp.PIPE,stdout=sp.PIPE, shell=True)

    out, err = job.communicate()
    print('out:', out.decode())
    print('err:', err.decode())

    print('Created the gzfn: {}'.format(fn + '.gz'))
    print('Created the tabix: {}'.format(fn + '.gz.tbi'))
# make the longrange file
lrange_fn = os.path.join(outdir, 'gs_loops.longrange.bed')
bedpe_to_WashU_longrange(lrange_fn, gs_loops)
gs_loops.gene_name.unique()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1467483108108108, 'fraction_numerical': 0.013724662162162162, 'mean_word_length': 3.5504322766570606, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '569292', 'license': 'MIT', 'n_tokens_mistral': 1843, 'n_tokens_neox': 1775, 'n_words': 361, 'path': 'scripts/loop_analysis/archive/HiChip_Loop_with_Colocs_and_Novel_Genes.ipynb', 'repo_name': 'joreynajr/dchallenge', 'size': '4736'}