上一篇介绍了使用ChIPseeker对peaks的分布和邻近基因的注释。除此之外，下游分析通常还包括鉴定我们感兴趣的蛋白质结合的motif。
motif是比较有特征的短序列，会多次出现的，一般认为它的生物学意义重大，做完CHIP-seq分析之后，一般都会寻找motif 。
motif的定义：recurring pattern. eg, sequence motif, structure motif or network motif
DNA sequence motif: short, recurring patterns in DNA that are presumed to have a biological function.从定义可以看出，其实motif这个单词就是形容一种反复出现的模式，而序列motif往往是DNA上的反复出现的模式，并被假设拥有生物学功能。而且，经常是一些具有序列特异性的蛋白的结合位点（如，转录因子）或者是涉及到重要生物过程的（如，RNA 起始，RNA 终止， RNA 剪切等等）
motif最先是通过实验的方法发现的，换句话说，不是说有了ChIP-seq才有了motif分析，起始很早人们就开始研究motif了！例如，'TATAAT’ box在1975年就被pribnow发现了，它与上游的‘TTGACA’motif是RNA聚合酶结合位点的特异性序列。而且，当时的人们就知道，不是所有的结合位点都一定完美地与motif匹配，大部分都只匹配了12个碱基中的7-9个。结合位点与motif的匹配程度往往也与蛋白质与DNA的结合强弱有关。
1.首先对bed文件进行处理
用MACS2软件进行peak calling，也就是找比对结果(bam文件)的peaks，MACS2找到的peaks会存放在生成的*.bed文件里。homer软件找motif需要读取我们这里得到的bed格式的peaks文件。homer软件不仅可以找到motif，还可以注释peaks。
homer软件在读取bed文件时，需要提取对应的列作为输入文件
我们要对MACS找到的peaks记录文件，还需提取对应的列给HOMER作为输入文件，提取操作为：
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' sample_peaks.bed >sample_homer.bed
不理解没关系，我们拿我跑的chipseq数据做示范，看awk对bed文件做了什么
使用命令less -S MUT_rep1_summits.bed 看下awk处理之前bed文件长什么样子：
less -S idr_peak.narrowPeak
chr1 121483892 121485575 . 1000 . -1 114.83000 -1 1257 4.128463 4.128463 121483919
chr19 27735849 27736995 . 1000 . -1 66.73870 -1 572 3.411804 3.636552 27736099
使用awk对bed文件进行处理，处理完的文件命名为
# summit 列 为相对于peak起始位点偏移多少碱基是 summit(可选步骤，更严格筛选步骤)
awk '{m1=$13+$16; m2=$17+$20; if((m1-m2) <= 50 && (m1-m2) >= -50) {m=int((m1+m2)/2); print $1, m, m+1}}' idr_peak.narrowPeak | wc -l
# 构造 融合的 peaks 结果
awk 'BEGIN{FS="\t";OFS="\t"}{m1=$13+$16; m2=$17+$20; if((m1-m2) <= 50 && (m1-m2) >= -50) {m=int((m1+m2)/2); print $1,$2,$3,$5,$6,$7,$8,$9,m-$2}}' idr_peak.narrowPeak | sort -k1,1 -k2,2n | awk 'BEGIN{FS="\t";OFS="\t"}{print $1,$2,$3,"MUT_peak_"NR,$4,$5,$6,$7,$8,$9}' >MUT_IDR_results.reshapedAndSorted.peaks.bed
# 所以将原始peak summit在 两次重复中相差50nt以内的summit 合并取其中点构造 summit 信息
awk 'BEGIN{FS="\t";OFS="\t"}{m1=$13+$16; m2=$17+$20; if((m1-m2) <= 50 && (m1-m2) >= -50) {m=int((m1+m2)/2); print $1,m,m+1,$5,$6,$7,$8,$9}}' idr_peak.narrowPeak |sort -k1,1 -k2,2n | awk 'BEGIN{FS="\t";OFS="\t"}{print $1,$2,$3,"CEBPB_bSite_summit_"NR,$4,$5,$6,$7,$8}' >MUT_IDR_results.reshapedAndSorted.summit.bed
#取MUT_IDR_results.reshapedAndSorted.summit.bed的summit上下15bp，前3000个peak
sort -k7,7nr MUT_IDR_results.reshapedAndSorted.peaks.bed | head -n 3000 | cut -f 1,2,3 | awk 'BEGIN{FS="\t";OFS="\t"}{print $1,$2-15,$3+15,"MUT_peak_summit_eb15_"NR}' >MUT_IDR_results.reshapedAndSorted.top3000summit.extendB15.bed
#现在文件变成了这种格式
less -S MUT_IDR_results.reshapedAndSorted.top3000summit.extendB15.bed
chr10 42383091 42383888 MUT_peak_summit_eb15_1
chr10 42388920 42389779 MUT_peak_summit_eb15_2
chr10 42399913 42400846 MUT_peak_summit_eb15_3
把bed文件变成这个样子，用于Homer软件的读取,我们称这样的文件叫做 Homer Peak/Positions file
我们可以看到Homer Peak/Positions file有4列：
第一列: 染色体
第二列: 起始位置
第三列: 终止位置
第四列: 链的方向(+/- or 0/1, where 0=”+”, 1=”-“)
2.取fasta格式
bedtools getfasta -fi /home/zhangyihui/my_data/ref_data/hg19/hg19.fa -bed MUT_IDR_results.reshapedAndSorted.top3000summit.extendB15.bed -fo top3000summit.extendB15.fa -name
3.将fasta文件上传到meme获取motif序列
1.打开MEME主页
2.依次点击
image.png
3.传入fasta序列
image.png
4.
image.png
5.
image.png
网友评论