这里是佳奥！
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至此结束。
学习无止境，我们下一个篇章再见！
