<jupyter_start><jupyter_text>How to: assign TF Motifs to ChIP-seq peaksThis tutorial demonstrates one way to assign CTCF motifs to CTCF ChIP-seq peaks using bioframe.<jupyter_code>import bioframe
import os
import subprocess
import numpy as np
import pandas as pd 
import matplotlib.pyplot as plt
from scipy.stats import pearsonr, spearmanr
base_dir = '/tmp/bioframe_tutorial_data_v2/'
assembly = 'GRCh38'<jupyter_output><empty_output><jupyter_text>Load CTCF ChIP-seq peaks for HFF from ENCODEThis approach makes use of the `narrowPeak` schema for bioframe.read_table .<jupyter_code>ctcf_peaks = bioframe.read_table("https://www.encodeproject.org/files/ENCFF401MQL/@@download/ENCFF401MQL.bed.gz", schema='narrowPeak')
ctcf_peaks[0:5]<jupyter_output><empty_output><jupyter_text>Get CTCF motifs from JASPAR<jupyter_code>### CTCF motif: http://jaspar.genereg.net/matrix/MA0139.1/
jaspar_url = 'http://expdata.cmmt.ubc.ca/JASPAR/downloads/UCSC_tracks/2022/hg38/'
jaspar_motif_file = 'MA0139.1.tsv.gz'
ctcf_motifs = bioframe.read_table(jaspar_url+jaspar_motif_file,schema='jaspar',skiprows=1)    
ctcf_motifs[0:4]<jupyter_output><empty_output><jupyter_text>Overlap peaks & motifs<jupyter_code>df_peaks_motifs = bioframe.overlap(ctcf_peaks,ctcf_motifs, suffixes=('_1','_2'))<jupyter_output><empty_output><jupyter_text>There are often multiple motifs overlapping one ChIP-seq peak:<jupyter_code>motifs_per_peak = df_peaks_motifs.groupby(["chrom_1", "start_1","end_1"]).size().values
plt.hist(motifs_per_peak,np.arange(0,np.max(motifs_per_peak)))
plt.xlabel('number of overlapping motifs per peak')
plt.ylabel('number of peaks')
plt.semilogy();<jupyter_output><empty_output><jupyter_text>assign the strongest motif to each peak<jupyter_code># since idxmax does not currently take NA, fill with -1
df_peaks_motifs['pval_2'] = df_peaks_motifs['pval_2'].fillna(-1) 
idxmax_peaks_motifs = df_peaks_motifs.groupby(["chrom_1", "start_1","end_1"])["pval_2"].idxmax().values
df_peaks_maxmotif = df_peaks_motifs.loc[idxmax_peaks_motifs]
df_peaks_maxmotif['pval_2'].replace(-1,np.nan,inplace=True)<jupyter_output><empty_output><jupyter_text>stronger peaks tend to have stronger motifs:<jupyter_code>plt.rcParams['font.size']=12
df_peaks_maxmotif['fc_1']  = df_peaks_maxmotif['fc_1'].values.astype('float')
plt.scatter(df_peaks_maxmotif['fc_1'].values, 
            df_peaks_maxmotif['pval_2'].values, 5, alpha=0.5,lw=0)
plt.xlabel('ENCODE CTCF peak strength, fc')
plt.ylabel('JASPAR CTCF motif strength \n (-log10 pval *100)')
plt.title('corr: '+str(np.round(df_peaks_maxmotif['fc_1'].corr(df_peaks_maxmotif['pval_2']),2)));<jupyter_output><empty_output><jupyter_text>We can also ask the reverse question: how many motifs overlap a ChIP-seq peak?<jupyter_code>df_motifs_peaks = bioframe.overlap(ctcf_motifs,ctcf_peaks,how='left', suffixes=('_1','_2'))
m = df_motifs_peaks.sort_values('pval_1')
plt.plot( m['pval_1'].values[::-1] ,
    np.cumsum(pd.isnull(m['chrom_2'].values[::-1])==0)/np.arange(1,len(m)+1))
plt.xlabel('pval')
plt.ylabel('probability motif overlaps a peak');<jupyter_output><empty_output><jupyter_text>filter peaks overlapping blacklisted regionsdo any of our peaks overlap blacklisted genomic regions?<jupyter_code>blacklist = bioframe.read_table('https://www.encodeproject.org/files/ENCFF356LFX/@@download/ENCFF356LFX.bed.gz',
                                schema='bed3')
blacklist[0:3]<jupyter_output><empty_output><jupyter_text>there appears to be a small spike in the number of peaks close to blacklist regions<jupyter_code>import matplotlib.pyplot as plt
closest_to_blacklist = bioframe.closest(ctcf_peaks,blacklist)
plt.hist(closest_to_blacklist['distance'].astype('Float64').astype('float'),np.arange(0,1e4,100));<jupyter_output><empty_output><jupyter_text>to be safe, let's remove anything +/- 1kb from a blacklisted region<jupyter_code># first let's select the columns we want for our final dataframe of peaks with motifs
df_peaks_maxmotif = df_peaks_maxmotif[
                            ['chrom_1','start_1','end_1','fc_1',
                             'chrom_2','start_2','end_2','pval_2','strand_2']]
#  then rename columns for convenience when subtracting
for i in df_peaks_maxmotif.keys():
    if '_1' in i: df_peaks_maxmotif.rename(columns={i:i.split('_')[0]},inplace=True)

# now subtract, expanding the blacklist by 1kb        
df_peaks_maxmotif_clean = bioframe.subtract(df_peaks_maxmotif,bioframe.expand(blacklist,1000))<jupyter_output><empty_output><jupyter_text>there it is! we now have a dataframe containing positions of CTCF ChIP peaks, including the strongest motif underlying that peak, and after conservativefiltering for proximity to blacklisted regions<jupyter_code>df_peaks_maxmotif_clean.iloc[0:3]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12728418399495905, 'fraction_numerical': 0.023104389834068473, 'mean_word_length': 7.325174825174825, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 49, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '142445', 'license': 'MIT', 'n_tokens_mistral': 1855, 'n_tokens_neox': 1765, 'n_words': 302, 'path': 'docs/tutorials/tutorial_assign_motifs_to_peaks.ipynb', 'repo_name': 'gfudenberg/bioframe', 'size': '4761'}