这里是佳奥！

ChIP-seq分析的最后一篇。让我们开始吧。

序列motif往往是DNA上的反复出现的模式，并被假设拥有生物学功能。而且，经常是一些具有序列特异性的蛋白的结合位点（如，转录因子）或者是涉及到重要生物过程的（如，RNA 起始，RNA 终止， RNA 剪切等等）。

1 homer寻找motif

由于需要下载很大的关联数据库，这里仅作展示。

##软件安装
conda install -c bioconda homer

perl ~/miniconda3/envs/chipseq/share/homer-4.9.1-5/configureHomer.pl  -install mm10 
ls -lh  ~/miniconda3/envs/chipseq/share/homer-4.9.1-5/data/
##我们上游分析是基于mm10找到的peaks文件
##下载成功后会多出 ~/miniconda3/envs/chipseq/share/homer-4.9.1-5/data/genomes/mm10/ 文件夹, 共 4.9G
##这个文件夹取决于你把homer这个软件安装到了什么地方。

## 或者用下面代码安装:
cd ~/biosoft
mkdir homer &&  cd homer
wget http://homer.salk.edu/homer/configureHomer.pl 
perl configureHomer.pl -install
perl configureHomer.pl -install hg19 

homer软件找motif整合了两个方法：依赖于数据库的查询，和de novo的推断。

都是读取ChIP-seq数据上游分析得到的bed格式的peaks文件。

##使用homer：http://homer.ucsd.edu/homer/ngs/peakMotifs.html

cd  ~/project/epi/motif  
for id in /home/jmzeng/project/epi/peaks/*.bed;
do
echo $id
file=$(basename $id )
sample=${file%%.*} 
echo $sample  
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' $id >homer_peaks.tmp  
findMotifsGenome.pl homer_peaks.tmp mm10 ${sample}_motifDir -len 8,10,12
annotatePeaks.pl    homer_peaks.tmp mm10  1>${sample}.peakAnn.xls 2>${sample}.annLog.txt 
done 

把上面的代码保存为脚本runMotif.sh，然后运行：bash runMotif.sh 1>motif.log &

这样不仅找了motif，还把peaks注释了一下。得到的后缀名为peakAnn.xls 的文件。

这和使用R包得到的结果相似。

2 meme寻找motif

需要通过bed格式的peaks的坐标来获取fasta序列。

2.1 首先用.bed文件生成fasta文件

##usage: Rscript  peakView.R  peaks.bed  IP.sorted.bam  input.sorted.bam  10
#options(echo=TRUE) # if you want see commands in output file
args <- commandArgs(trailingOnly = TRUE)
if(length(args) != 1 ){
    print(" usage: Rscript  peakAnno.R  peaks.bed ")
}

bedPeaksFile        = args[1] ;

##自这开始，.bed文件要和R Project文件在同一目录下
bedFiles=list.files(pattern = '*.bed')
[1] "Control_summits.bed"          "H2Aub1_summits.bed"          
[3] "H3K36me3_summits.bed"         "Ring1B_summits.bed"          
[5] "RNAPII_8WG16_summits.bed"     "RNAPII_S2P_summits.bed"      
[7] "RNAPII_S5P_summits.bed"       "RNAPII_S5PRepeat_summits.bed"
[9] "RNAPII_S7P_summits.bed"   

BiocManager::install("BSgenome.Mmusculus.UCSC.mm10")
library(BSgenome.Mmusculus.UCSC.mm10)
library(ChIPpeakAnno)

##生成.fa文件
bedPeaksFile=bedFiles[2]##第二个文件即H2Aub1_summits.bed，要下一个就[3]
sampleName=strsplit(bedPeaksFile,'\\.')[[1]][1]
peak <- toGRanges(bedPeaksFile, format="BED") 
keepChr= !grepl('_',seqlevels(peak))
#seqlevels(peak, force=TRUE) <- seqlevels(peak)[keepChr]

seq <- getAllPeakSequence(peak, upstream=20, downstream=20, genome=Mmusculus)
write2FASTA(seq,  paste0(sampleName,'.fa'))

2.2 MEME（web）

https://meme-suite.org/meme/

QQ截图20220814161012.png QQ截图20220814161049.png

默认，在选择文件处上传fa文件。

QQ截图20220814161159.png

随后Start Search提交

QQ截图20220814161222.png

等待，出现该页面说明运行结束，查看HTML结果，或下载output压缩文件。

QQ截图20220814161316.png QQ截图20220814161400.png

可以查看到motif结果。

2.3 meme(Linux)

##摘抄自网页结果
meme ./seqs-centered -oc meme_out -mod zoops -nmotifs 3 -minw 6 -maxw 15 -bfile ./background -dna -searchsize 100000 -time 2919 -revcomp -nostatus

ChIP-seq至此结束。

学习无止境，我们下一个篇章再见！