<jupyter_start><jupyter_text>List of cells shared by a colocalized SNP and fithichip loop<jupyter_code>summary[summary.fithichip == 1]
shared_sgloops = summary[summary.fithichip == 1].groupby(['rs_id', 'near_gene'])
shared_sgloops = summary[summary.fithichip == 1].groupby(['rs_id', 'near_gene'])

shared_sgloops = shared_sgloops.cell.value_counts()
shared_sgloops
shared_sgloops = summary[summary.fithichip == 1].groupby(['rs_id', 'near_gene'])
shared_sgloops = summary[summary.fithichip == 1].groupby(['rs_id', 'near_gene'])

shared_sgloops = shared_sgloops.cell.unique()
shared_sgloops = shared_sgloops.to_frame().reset_index()
shared_sgloops = shared_sgloops.merge(snp_info, on='rs_id')
new_order = ['rs_id',
 'geneName',
 'near_gene',
 'cell',
 'chr',
 'pos',
 'pp_H0_Coloc_Summary',
 'pp_H1_Coloc_Summary',
 'pp_H2_Coloc_Summary',
 'pp_H3_Coloc_Summary',
 'pp_H4_Coloc_Summary',
 'variant_id',
 'dist',
 'pvalue',
 'FDR',
 'slope_snp',
 'ref',
 'alt',
 'AC',
 'AF',
 'AN',
 'slope_se_snp',
 'slope_gwas',
 'slope_se_gwas',
 'pval_nominal',
 'SampleSize']
shared_sgloops = shared_sgloops.loc[:, new_order]
shared_sgloops.rename(columns={'geneName': 'egene'}, inplace=True)
shared_sgloops.set_index(['rs_id', 'egene'])<jupyter_output><empty_output><jupyter_text>Make the master table<jupyter_code>import os 
import sys
import pybedtools as pbt
import pandas as pd
import numpy as np
import subprocess as sp
import json
os.chdir('/mnt/BioHome/jreyna/jreyna/projects/dchallenge/')
pbt.set_bedtools_path('/mnt/BioApps/bedtools/bin/')
bgzip = '/mnt/BioApps/tabix/tabix-0.2.6/bgzip'
tabix = '/mnt/BioApps/tabix/tabix-0.2.6/tabix'

bedpe_6cols = ['chrA', 'startA', 'endA', 'chrB', 'startB', 'endB']
bedpe_10cols = ['chrA', 'startA', 'endA', 'chrB', 'startB', 'endB', 'name', 'score', 'strand1', 'strand2']

# default values for the command line
sys.argv = [0] * 7
sys.argv[1] =  'results/main/2021_Nikhil_eQTL/Results/Colocalization/T1D_34012112_Gaulton/'
sys.argv[1] += 'BLUEPRINT_eQTL_Monocyte/FINAL_Summary_Coloc_Gene_SNP_Pairs.bed'
sys.argv[2] = 'results/refs/ensembl/gencode.v19.annotation.bed'
sys.argv[3] = 'results/main/2021_Nikhil_eQTL/Data/FitHiChIP_Loops/CM/FitHiChIP_L/FitHiChIP.interactions_FitHiC_Q0.01.bed'
sys.argv[4] = 'results/refs/spp/SPP_D-Challenge_networks.xlsx'
sys.argv[5] = 'results/refs/hg19/hg19.chrom.sizes'
sys.argv[6] = 'results/main/loop_analysis/washU/'

# parsing the commandline arguments
coloc_fn = sys.argv[1]
genes_fn = sys.argv[2]
loop_fn = sys.argv[3]
spp_fn = sys.argv[4]
gs_fn = sys.argv[5]
outdir = sys.argv[6]

# setting the output file names
os.makedirs(outdir, exist_ok=True)
summary_fn = os.path.join(outdir, 'sgloop_summary.xlsx')
sg_pairs_fn = os.path.join(outdir, 'gs_pairs.longrange.bed')
sg_loops_fn = os.path.join(outdir, 'gs_loops.longrange.bed')<jupyter_output><empty_output><jupyter_text>Load the colocalization data<jupyter_code># load the colocalization data
coloc = pd.read_table(coloc_fn)

# extract the most significant according the H4 
coloc_sig_df = coloc[coloc['pp_H4_Coloc_Summary'] > 0.75]
coloc_sig_df.rename(columns={'pos': 'end'}, inplace=True)
coloc_sig_df.loc[:, 'start'] = coloc_sig_df.loc[:, 'end'] - 1
coloc_sig_full = coloc_sig_df.copy(deep=True)
coloc_sig_df = coloc_sig_df[['chr', 'start', 'end', 'rs_id', 'variant_id']]
coloc_sig_df = coloc_sig_df.loc[~coloc_sig_df.duplicated(subset='rs_id'),]
coloc_sig_pbt = pbt.BedTool.from_dataframe(coloc_sig_df.iloc[:, 0:4]).sort()
#csnp_slop_pbt = coloc_sig_pbt.slop(b=500000, g=gs_fn)<jupyter_output><empty_output><jupyter_text>Load the gene data<jupyter_code># load the gencode coords
cols = ['chrom', 'start', 'end', 'strand', 'type', 'gene_id', 'gene_name']
gencode = pd.read_table(genes_fn, header=None, names=cols)

# extract just the genes
genes_df = gencode.loc[gencode.type.isin(['gene'])]
genes_df = genes_df.loc[~genes_df.duplicated(subset='gene_id'), :]
genes_df.loc[:, 'chrom'] = genes_df['chrom'].astype(str)
genes_df = genes_df.iloc[:, [0,1,2,6,5]]
genes_pbt = pbt.BedTool.from_dataframe(genes_df).sort()
genes_df.head()<jupyter_output><empty_output><jupyter_text>Find the closest gene<jupyter_code>closest_gene = coloc_sig_pbt.closest(genes_pbt, d=True)
closest_gene = closest_gene.to_dataframe()
closest_gene = closest_gene.iloc[:, [3,7,8,9]]
closest_gene.columns = ['rs_id', 'cls_gname', 'cls_id', 'cls_dist']
closest_gene.head()
uniq_cls_gname = closest_gene.groupby(['rs_id']).cls_gname.apply(lambda x: ','.join(x))
uniq_cls_ids = closest_gene.groupby(['rs_id']).cls_id.apply(lambda x: ','.join(x))
uniq_cls_dist = closest_gene.groupby(['rs_id']).cls_dist.apply(lambda x: ','.join([str(i) for i in x]))
uniq_cls = pd.merge(uniq_cls_gname, uniq_cls_ids, left_index=True, right_index=True)
uniq_cls = pd.merge(uniq_cls, uniq_cls_dist, left_index=True, right_index=True)
uniq_cls.reset_index(inplace=True)
uniq_cls<jupyter_output><empty_output><jupyter_text>Find all genes +/- 500kb<jupyter_code># get a list of gene names within +- 500kb of the SNPs
fivekb_gnames = coloc_sig_pbt.slop(b=500000, g=gs_fn)
fivekb_gnames = fivekb_gnames.map(genes_pbt, c=4, o='collapse')
fivekb_gnames = fivekb_gnames.to_dataframe()
fivekb_gnames = fivekb_gnames.iloc[:, [3,4]]
fivekb_gnames.columns = ['rs_id', 'gene_name']

# get a list of gene ids within +- 500kb of the SNPs
fivekb_gids = coloc_sig_pbt.slop(b=500000, g=gs_fn)
fivekb_gids = fivekb_gids.map(genes_pbt, c=5, o='collapse')
fivekb_gids = fivekb_gids.to_dataframe()
fivekb_gids = fivekb_gids.iloc[:, [3,4]]
fivekb_gids.columns = ['rs_id', 'gene_id']

# merge the two above results
fivekb_genes = fivekb_gnames.merge(fivekb_gids, on='rs_id')
fivekb_genes.columns = ['rs_id', '5kb_gname', '5kb_gid']
# get eQTL's
eqtl_fn = 'results/main/2021_Nikhil_eQTL/Data/eqtl_sqtl_summ_stats/BLUEPRINT_eQTL/Monocyte.txt.gz'
eqtls = pd.read_table(eqtl_fn)
eqtls.columns = ['eqtl_gname', 'nvar', 'shape1', 'shape2', 'dummy',
                 'sid', 'dist', 'npval', 'slope', 'ppval', 'bpval', 'qval']<jupyter_output><empty_output><jupyter_text>Generate the master table<jupyter_code>master = coloc_sig_full.copy()

# add sid which is the chr:position of the SNP
master['sid'] = master['chr'].str.replace('chr', '') + ':' +  master['end'].astype(str)

# add the closest gene
master = master.merge(uniq_cls, on='rs_id', how='left')

# add the +/- fivekb genes 
master = master.merge(fivekb_genes, on='rs_id', how='left')

# add the eQTL data
eqtl_genes = master.merge(eqtls[['sid', 'eqtl_gname']], on='sid')
eqtl_genes = eqtl_genes.groupby('rs_id').eqtl_gname.unique()
eqtl_genes = eqtl_genes.apply(lambda x: ','.join(x))
master = master.merge(eqtl_genes, on='rs_id', how='left')

new_colnames = [
 'rs_id',
 'variant_id',
 'chr',
 'start',
 'end',           
 'geneName',  
 'eqtl_gname',
 'cls_gname',

 '5kb_gname',

 'pp_H0_Coloc_Summary',
 'pp_H1_Coloc_Summary',
 'pp_H2_Coloc_Summary',
 'pp_H3_Coloc_Summary',
 'pp_H4_Coloc_Summary',           
 'dist',
 'pvalue',
 'FDR',
 'slope_snp',
 'ref',
 'alt',
 'AC',
 'AF',
 'AN',
 'slope_se_snp',
 'slope_gwas',
 'slope_se_gwas',
 'pval_nominal',
 'SampleSize',
 'sid',
 'cls_id',
 'cls_dist',
 '5kb_gid'
]
master = master.loc[:, new_colnames]
master.rename(columns={'geneName': 'coloc_gname',
                       'end': 'pos',
                       'eqtl_gname': 'eqtl_gnames', 
                       'cls_gname': 'cls_gnames', 
                       'cls_id': 'cls_ids'}, inplace=True)
master.drop(['start'], axis=1, inplace=True)
master<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16773504273504272, 'fraction_numerical': 0.017094017094017096, 'mean_word_length': 6.37832512315271, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 30, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '392368', 'license': 'MIT', 'n_tokens_mistral': 3289, 'n_tokens_neox': 3172, 'n_words': 486, 'path': 'scripts/loop_analysis/archive/Annotation_SGLoops.ipynb', 'repo_name': 'joreynajr/dchallenge', 'size': '7488'}