Setup the data¶
Focus on the 200bp window around the peaks.
In [5]:
from basepair.config import get_data_dir
import pyBigWig
import matplotlib.pyplot as plt
from concise.utils.fasta import read_fasta
from pybedtools import BedTool
from tqdm import tqdm
import pandas as pd
import numpy as np
In [6]:
from basepair.data import get_sox2_data
In [7]:
ddir = get_data_dir()
Data¶
Peak sites¶
In [4]:
# data:
fasta_file_200 = f"{ddir}/raw/chipnexus/scratch_jzeitlinger_collab/meme_chip_analysis/Sox2_123b_1_ppr.IDR0.05.filt.summit_centered_200bp.fasta"
bed_file_200 = f"{ddir}/raw/chipnexus/scratch_jzeitlinger_collab/meme_chip_analysis/Sox2_123b_1_ppr.IDR0.05.filt.summit_centered_200bp.narrowPeak"
# Sox2_123b_1_ppr.IDR0.05.filt.summit_centered_200bp.fasta
fas = read_fasta(fasta_file_200)
bed = BedTool(bed_file_200)
assert len(bed) == len(fas)
ChipNexus¶
In [17]:
dfc = get_sox2_data()
Averaged data signal¶
In [21]:
plt.plot(np.arange(-100.5, 100.5), np.stack(dfc.cuts_pos).mean(axis=0), label='pos')
plt.plot(np.arange(-100.5, 100.5), np.stack(dfc.cuts_neg).mean(axis=0), label='neg')
plt.xlabel("Peak position")
plt.legend()
plt.title("Average signal");
Raw data at random sites¶
In [8]:
n=2*6
fig = plt.figure(figsize=(20, 5))
for i,idx in enumerate(dfc.sample(n).index):
plt.subplot(n//6, 6, i+1)
plt.bar(np.arange(-100.5, 100.5), dfc.iloc[idx].cuts_pos, label='pos')
plt.bar(np.arange(-100.5, 100.5), dfc.iloc[idx].cuts_neg, label='neg')
if i==0:
plt.title("Average signal")
plt.xlabel("Peak position")
plt.legend()
plt.tight_layout()
In [9]:
dfc.sample(10)
Out[9]:
In [24]:
len(dfc.iloc[0].cuts_pos)
Out[24]:
In [25]:
len(dfc.iloc[0].seq)
Out[25]:
In [27]:
len(dfc)
Out[27]:
Motif sites¶
In [8]:
motif_sites = {"motif1": f"{ddir}/processed/chipnexus/motifs/sox2/homer_200bp/de-novo/motif1.peaks200.txt",
"motif2": f"{ddir}/processed/chipnexus/motifs/sox2/homer_200bp/de-novo/motif2.peaks200.txt"}
In [12]:
motifs = {"motif1": f"{ddir}/processed/chipnexus/motifs/sox2/homer_200bp/de-novo/motif1.motif",
"motif2": f"{ddir}/processed/chipnexus/motifs/sox2/homer_200bp/de-novo/motif2.motif"}
In [13]:
from pybedtools import BedTool
from pysam import Fastafile
import pandas as pd
from concise.utils.pwm import PWM,load_motif_db
import concise
import matplotlib.pyplot as plt
from plotnine import *
In [14]:
from basepair.motif.homer import load_motif, read_motif_hits
pwm_list = [load_motif(fname) for k,fname in motifs.items()]
In [16]:
for i,pwm in enumerate(pwm_list):
pwm.plotPWM((5,1.5))
plt.title(f"motif{i+1}")
In [17]:
for i,pwm in enumerate(pwm_list):
pwm.plotPWMInfo((5,1.5))
plt.title(f"motif{i+1}")
In [18]:
df = pd.concat([read_homer_motif_hits(fname).assign(motif_id=k) for k,fname in motif_sites.items()])
In [19]:
df.head()
Out[19]:
In [20]:
for i,k in enumerate(motif_sites):
plt.hist(df.center[df.motif_id==k], bins=20, histtype="step", label=k)
plt.xlabel("Relative position");
plt.legend()
plt.title("Relative position wrt sequence");
Signal distribution at motif sites¶
In [21]:
dfcm = df.merge(dfc, left_on='interval_id', right_on='seq_id', suffixes=('_motif', '_interval'))
In [22]:
dfcm.shape
Out[22]:
In [23]:
dfcm.head()
Out[23]:
In [24]:
cuts_arr_pos = np.stack(dfcm.cuts_pos)
cuts_arr_neg = np.stack(dfcm.cuts_neg)
In [25]:
from tqdm import trange
In [40]:
def extract_signal_pos(row, flip_strand=True, pad=10):
out = row.cuts_pos[(row.start_motif-pad):(row.end_motif+pad)]
if flip_strand and row.strand == "-":
return out[::-1]
else:
return out
In [27]:
def skip_row(row, pad=10):
return row.start_motif < pad or row.end_motif + pad > 200
In [41]:
def extract_signal_neg(row, flip_strand=True, pad=10):
out = row.cuts_neg[(row.start_motif-pad):(row.end_motif+pad)]
if flip_strand and row.strand == "-":
return out[::-1]
else:
return out
In [49]:
def extract_signal(row, flip_strand=True, pad=10):
out_pos = row.cuts_pos[(row.start_motif-pad):(row.end_motif+pad)]
out_neg = row.cuts_neg[(row.start_motif-pad):(row.end_motif+pad)]
if flip_strand and row.strand == "-":
return np.array([out_neg[::-1], out_pos[::-1]])
else:
return np.array([out_pos, out_neg])
In [29]:
pad=20
In [52]:
cuts_motif = np.stack([extract_signal(dfcm.iloc[i], pad=pad) for i in range(len(dfcm)) if not skip_row(dfcm.iloc[i], pad=pad)])
In [42]:
cuts_motif_pos = np.stack([extract_signal_pos(dfcm.iloc[i], pad=pad) for i in trange(len(dfcm)) if not skip_row(dfcm.iloc[i], pad=pad)])
cuts_motif_neg = np.stack([extract_signal_neg(dfcm.iloc[i], pad=pad) for i in trange(len(dfcm)) if not skip_row(dfcm.iloc[i], pad=pad)])
In [31]:
len(cuts_motif_pos)
Out[31]:
In [32]:
# We haven't lost more than 80% of the sequences
assert len(cuts_motif_pos) > 0.8*len(dfcm)
In [54]:
motif_id = np.array([dfcm.iloc[i].motif_id for i in range(len(dfcm)) if not skip_row(dfcm.iloc[i], pad=pad)])
In [55]:
cuts_motif_pos.shape
Out[55]:
In [58]:
# Strand flipped and reversed
plt.figure(figsize=(10,3))
plt.subplot(121)
plt.plot(cuts_motif[motif_id=='motif1',0].mean(axis=0), label='pos')
plt.plot(cuts_motif[motif_id=='motif1',1].mean(axis=0), label='neg')
plt.title("motif1")
plt.subplot(122)
plt.plot(cuts_motif[motif_id=='motif2',0].mean(axis=0), label='pos')
plt.plot(cuts_motif[motif_id=='motif2',1].mean(axis=0), label='neg')
plt.legend()
plt.title("motif2");
In [44]:
# Strand flipped
plt.figure(figsize=(10,3))
plt.subplot(121)
plt.plot(cuts_motif_pos[motif_id=='motif1'].mean(axis=0), label='pos')
plt.plot(cuts_motif_neg[motif_id=='motif1'].mean(axis=0), label='neg')
plt.title("motif1")
plt.subplot(122)
plt.plot(cuts_motif_pos[motif_id=='motif2'].mean(axis=0), label='pos')
plt.plot(cuts_motif_neg[motif_id=='motif2'].mean(axis=0), label='neg')
plt.legend()
plt.title("motif2");
In [35]:
# Strand not flipped
plt.figure(figsize=(10,3))
plt.subplot(121)
plt.plot(cuts_motif_pos[motif_id=='motif1'].mean(axis=0), label='pos')
plt.plot(cuts_motif_neg[motif_id=='motif1'].mean(axis=0), label='neg')
plt.title("motif1")
plt.subplot(122)
plt.plot(cuts_motif_pos[motif_id=='motif2'].mean(axis=0), label='pos')
plt.plot(cuts_motif_neg[motif_id=='motif2'].mean(axis=0), label='neg')
plt.legend()
plt.title("motif2");
In [36]:
for i,pwm in enumerate(pwm_list):
pwm.plotPWMInfo((5,1.5))
plt.title(f"motif{i+1}")