1. 環(huán)境配置

#!/bin/bash
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh 
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes

conda  create -n chipseq  python=2 bwa
conda info --envs
source activate chipseq
# 可以用search先進(jìn)行檢索
conda search trim_galore
## 保證所有的軟件都是安裝在 chipseq 這個(gè)環(huán)境下面
conda install -y sra-tools  
conda install -y trim-galore  samtools
conda install -y deeptools homer  meme
conda install -y macs2 bowtie bowtie2 

2. 數(shù)據(jù)下載

(1)prefetch下載

首先找到要下載的數(shù)據(jù)（SRR...號(hào)）
創(chuàng)建project文件夾

cd project/
#將下列數(shù)據(jù)寫(xiě)入vim文件中
vim SraAccList.txt
SRR1266976
SRR1266977
SRR1266978
SRR1266979
SRR1266980
SRR1266981

利用chipseq環(huán)境中prefetch 軟件進(jìn)行下載
如果不指定下載目錄 則默認(rèn)下載目錄：~/ncbi/public/sra/

source activate chipseq
mkdir {sra,bedgraph,fastq,rmdup,tss,clean,align,peaks,motif,qc/{raw,trimed},annotation} #在project目錄下新建目錄
#將數(shù)據(jù)下載于sra目錄下
cat SraAccList.txt | while read id;
do
(nohup prefetch -O ./sra $id &)#保存在sra目錄下
done
#control +c可以后臺(tái)下載

prefetch --option-file SraAccList.txt #或者批量下載可以點(diǎn)擊下載SraAccList.txt文件

prefetch常用選項(xiàng)

(2)Aspera Connect下載

Aspera—ascp下載SRA數(shù)據(jù)

SRA數(shù)據(jù)庫(kù)下載：數(shù)據(jù)的存放地址是ftp://ftp.sra.ebi.ac.uk/vol1
舉例：下載SRR1577019文件锤悄，首先我需要找到地址碟绑，去ftp://ftp.sra.ebi.ac.uk/vol1芦昔，一層層尋找伊者，直至找到素挽，ftp://ftp.sra.ebi.ac.uk/vol1/srr/SRR157/009/SRR1577019

去geo官網(wǎng)搜索srr然后下載獲取SRR_Acc_List.txt 然后觀察是err+7位數(shù)的褐捻，所以直接使用7位數(shù)的代碼偏螺。（7位數(shù)y取最后一位痪宰，8位數(shù)y取最后2位叼架，6位數(shù)把y變量刪掉，應(yīng)該是這樣衣撬。）一般來(lái)說(shuō)乖订，NCBI的sra文件前面的地址都是一樣的 ~/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk，那么寫(xiě)腳本批量下載也就不難了具练！最后那個(gè) ./ 表示下載后的路徑乍构。

1. 1批量下載srr文件（推薦）

cd ~/
mkdir -p project/{sra,bedgraph,fastq,rmdup,tss,clean,align,peaks,motif,qc/{raw,trimed},annotation}#創(chuàng)建文件夾

cat SraAccList.txt | while read id
do
x=$(echo $id | cut -b1-6)
y=$(echo $id | cut -b10-10)
echo $id
ascp -QT -k 1 -l 300m -P33001  -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/srr/$x/00$y/$id/ ~/project/sra/
done

1.2 單個(gè)下載srr文件

ascp -QT -k 1  -l 500m -P 33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/srr/SRR126/000/SRR1266980 ./sra

2. 下載fastq文件

cat SRR_Acc_List.txt | while read id
do
x=$(echo $id | cut -b1-6)
y=$(echo $id | cut -b10-10)
echo $id
ascp -QT -k 1 -l 300m -P33001  -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/$x/00$y/$id/ ./sra
done

如果下載的fastq數(shù)據(jù)是2個(gè)文件 *_1.fastq 和 *_2.fastq表明是雙端測(cè)序。

Aspera用法如下：
ascp [參數(shù)] 目標(biāo)文件 保存路徑
-v verbose mode 實(shí)時(shí)知道程序在干啥
-T 取消加密扛点，否則有時(shí)候數(shù)據(jù)下載不了
-i 提供私鑰文件的地址
-l 設(shè)置最大傳輸速度哥遮，一般200m到500m，如果不設(shè)置陵究，反而速度會(huì)比較低眠饮，可能有個(gè)較低的默認(rèn)值
-k 斷點(diǎn)續(xù)傳，一般設(shè)置為值1
-Q 一般加上它
-P 提供SSH port铜邮，端口一般是33001

3. sra文件轉(zhuǎn)fastq文件：

• sra轉(zhuǎn)化為fastq文件可以使用sratoolkit中的fastq-dump命令仪召。

fastq-dump --split-3 filename

其中 --split-3 參數(shù)代表著如果是單端測(cè)序就生成一個(gè) 、.fastq文件松蒜，如果是雙端測(cè)序就生成*_1.fastq 和*_2.fastq 文件扔茅；即單端測(cè)序輸出1個(gè)fastq文件，雙端測(cè)序輸出2個(gè)fastq文件秸苗。
進(jìn)入到sra文件中,將SRR文件夾去掉召娜，數(shù)據(jù)文件直接保存在sra目錄中，我們可以用下述代碼進(jìn)行批量的格式轉(zhuǎn)化：
首先制作如下對(duì)應(yīng)的文本

SRR8444250   ES_Input
SRR8444251   ES_Smad2_ChIPseq
SRR8444256   EB_AC_Input
SRR8444257   EB-SB_Smad2_ChIPseq
SRR8444258   EB-AC_Smad2_ChIPseq

然后

## 下面需要用循環(huán)
cd ~/project/sra/
#把list.txt文件復(fù)制到sra目錄下
source activate chipseq
## 下面用到的 list.txt 文件惊楼，就是上面自行制作的對(duì)應(yīng)的文件玖瘸。
cat list.txt | while read id
do 
echo $id
arr=($id)
srr=${arr[0]} #通過(guò)中括號(hào)獲取數(shù)組中的某一個(gè)元素：${arr[0]} 得到第一個(gè)元素吐句，${arr[0]} 第二個(gè)
sample=${arr[1]}
fastq-dump -A $sample -O ~/project/fastq/ --gzip --split-3  $srr
done 
#-A -$sample是輸出文件名字
 #--gzip打包節(jié)省空間，且對(duì)后續(xù)分析不影響店读，如果有影響就把這個(gè)參數(shù)去掉嗦枢，-O輸出到目錄

4. 使用trim_galore軟件對(duì)fastq進(jìn)行質(zhì)量控制-過(guò)濾

（1）基礎(chǔ)知識(shí)

trim_galore使用手冊(cè)

cutadapt 軟件可以對(duì)NGS數(shù)據(jù)進(jìn)行質(zhì)量過(guò)濾
FastQC 軟件可以查看NGS數(shù)據(jù)的質(zhì)量分布
trim_galore 將這兩個(gè)軟件封裝到一起，使用起來(lái)更加方便屯断。

• Trim galore文虏，是可以自動(dòng)檢測(cè)adapter
• 去除reads 3’端的低質(zhì)量堿基 # 自動(dòng)調(diào)用cutadapt
• 去除adapter序列 # 自動(dòng)去除3’端的adapter,可以通過(guò)設(shè)定--Illumina，-- small_rna殖演，--nextera參數(shù)來(lái)指定對(duì)應(yīng)的adapter類型
• 去除長(zhǎng)度太短的序列 #通過(guò)設(shè)定--length參數(shù)氧秘，小于設(shè)定值被去除
• 其它過(guò)濾

參考：https://cloud.tencent.com/developer/article/1625193
和http://www.reibang.com/p/7a3de6b8e503?from=groupmessage
這兩個(gè)鏈接要看。

說(shuō)明

1. Trim Galore是對(duì)FastQC和Cutadapt的包裝趴久。
2. trim_galore適用于所有高通量測(cè)序丸相，包括RRBS(Reduced Representation Bisulfite-Seq ), Illumina、Nextera和smallRNA測(cè)序平臺(tái)的雙端和單端數(shù)據(jù)彼棍。
3. 主要功能包括兩步：
   第一步 首先去除低質(zhì)量堿基灭忠，然后去除3' 末端的adapter(如果沒(méi)有指定具體的adapter，程序會(huì)自動(dòng)檢測(cè)前1million的序列）
   第二步 對(duì)比前12-13bp的序列是否符合以下類型的adapter :

• Illumina: AGATCGGAAGAGC # 如果輸入?yún)?shù) --Illumina,就會(huì)默認(rèn)trimmed前13bp的adapter
• Small RNA: TGGAATTCTCGG # 同上 如果輸入?yún)?shù) --small_rna,就會(huì)默認(rèn)trimmed前12bp的adapter
• Nextera: CTGTCTCTTATA # 同上 如果輸入?yún)?shù) --nextera,就會(huì)默認(rèn)trimmed前12bp的adapter

對(duì)于單端測(cè)序數(shù)據(jù)座硕，基本用法如下
trim_galore --quality 20 -a AGATCGGAAGAGC --length 20 -o out_dir input.fq
其中-a后面可以跟著序列（-aAGATCGGAAGAGC）

對(duì)于雙端測(cè)序數(shù)據(jù)弛作，基本用法如下
trim_galore --paired --quality 20 -a AGATCGGAAGAGC -a2 AGATCGGAAGAGC --length 20 -o out_dir R1.fq.gz R2.fq.gz

參數(shù)說(shuō)明

--quality/-q<INT>：設(shè)定Phred quality score閾值，默認(rèn)為Phred 20 切除質(zhì)量得分低于設(shè)定值的序列华匾。
--phred33/64：使用ASCII+33/64質(zhì)量得分作為Phred得分映琳，選擇-phred33或者-phred64，表示測(cè)序平臺(tái)使用的Phred quality score蜘拉。具體怎么選擇萨西，看你用什么測(cè)序平臺(tái)。
（需要確認(rèn):Sanger/Illumina 1.9+ encoding為33; Illumina 1.5 encoding為64）
--adapter/-a ：輸入adapter序列旭旭。也可以不輸入谎脯，Trim Galore!會(huì)自動(dòng)尋找可能性最高的平臺(tái)對(duì)應(yīng)的adapter。自動(dòng)搜選的平臺(tái)三個(gè)您机，也直接顯式輸入這三種平臺(tái)穿肄，即--illumina、--nextera和--small_rna际看。
-s/--stringency<INT>：設(shè)定可以忍受的前后adapter重疊的堿基數(shù)咸产，默認(rèn)為1（非常苛刻）仲闽∧砸纾可以適度放寬，因?yàn)楹笠粋€(gè)adapter幾乎不可能被測(cè)序儀讀到。
--length<INT>：設(shè)定輸出reads長(zhǎng)度閾值小于設(shè)定值會(huì)被拋棄屑彻，默認(rèn)值20bp验庙；即小于20bp的被去除。注意社牲，在pe150下粪薛，不要設(shè)置太高，可以50或36（默認(rèn)20）
--max_length : 設(shè)置長(zhǎng)度大于此值被丟棄
-e <ERROR RATE> ：默認(rèn)0.1
--paired：對(duì)于雙端測(cè)序結(jié)果搏恤，一對(duì)reads中违寿，如果有一個(gè)被剔除，那么另一個(gè)會(huì)被同樣拋棄熟空，而不管是否達(dá)到標(biāo)準(zhǔn)藤巢。
--retain_unpaired：對(duì)于雙端測(cè)序結(jié)果，一對(duì)reads中息罗，如果一個(gè)read達(dá)到標(biāo)準(zhǔn)掂咒，但是對(duì)應(yīng)的另一個(gè)要被拋棄，達(dá)到標(biāo)準(zhǔn)的read會(huì)被單獨(dú)保存為一個(gè)文件迈喉。
--gzip和--dont_gzip：清洗后的數(shù)據(jù)zip打包或者不打包绍刮。
-o/--output_dir：輸入目錄 [需要提前建立目錄，否則運(yùn)行會(huì)報(bào)錯(cuò)]弊添。
-- trim-n : 移除read一端的reads
-j :使用線程數(shù), 注意假設(shè)已使用Python 3并安裝了Pigz录淡，那么內(nèi)核設(shè)置為4捌木，實(shí)際使用內(nèi)核是15油坝，因此，最高設(shè)為4.
--fastqc #當(dāng)分析結(jié)束后刨裆，使用默認(rèn)選項(xiàng)對(duì)結(jié)果文件進(jìn)行fastqc分析

phred 33 和phred 64對(duì)照表

1. 如何判斷phred 33 和phred 64?
   有時(shí)候得到的原始fastq文件澈圈，無(wú)法知道質(zhì)量值體系，你就無(wú)法進(jìn)行質(zhì)量值的過(guò)濾帆啃，我們可以在正常情況下瞬女，按照上面的表對(duì)回去，統(tǒng)計(jì)一下幾條reads的最大和最小質(zhì)量值的區(qū)間努潘，就可以知道到底是phred 33 還是phred 64體系诽偷。
2. 根據(jù)上圖中Phred+33與Phred+64所使用的質(zhì)量字符范圍的不同，可以對(duì)fastq文件中質(zhì)量得分的編碼方式做出判斷(第四行是測(cè)序質(zhì)量疯坤，查看第四行符號(hào)是在Phred+64還是Phred+33范圍)报慕。圖中顯示，ASCII值小于等于58（相應(yīng)的質(zhì)量得分小于等于25）對(duì)應(yīng)的字符只有在Phred+33的編碼中被使用压怠，所有Phred+64所使用的字符的ASCII值都大于等于59眠冈。在通常情況下，ASCII值大于等于74的字符只出現(xiàn)在Phred+64中菌瘫。利用這些信息即可在程序中進(jìn)行判斷蜗顽。（近幾年應(yīng)該都是Phred+33）

（2）數(shù)據(jù)處理

• 處理單個(gè)雙端測(cè)序fastq文件處理

trim_galore -q 25 --phred33 --stringency 3 --length 36  --paired CK-4_1.fq.gz CK-4_2.fq.gz --gzip -o ./cleandata/trim_galoredata/

• 批量處理雙端測(cè)序數(shù)據(jù)

cd ~/project/clean
files=~/project/fastq/*.gz
ls $files | while read id
do
sample=$(basename $id)
fq=${sample%%_?.fastq*}
echo $fq
trim_galore -q 25 --phred33 --length 25 --stringency 4 --paired ~/project/fastq/${fq}_1.fastq.gz ~/project/fastq/${fq}_2.fastq.gz --gzip -o ~/project/clean/
done

單端測(cè)序結(jié)果處理

cd ~/project/clean
files=~/project/fastq/*.gz
ls $files | while read id
do
sample=$(basename $id)
echo $sample
trim_galore -q 25 --phred33 --length 25 --stringency 4 ~/project/fastq/${sample} -o ~/project/clean/ 
done 
#具體稍作調(diào)整

雙端測(cè)序結(jié)果主文件為 *_1_val_1.fq.gz 和 *_2_val_2.fq.gz布卡，單端測(cè)序結(jié)果為 *_trimmed.fq。

(3)trim_galore處理結(jié)束雇盖，利用 --fastqc進(jìn)行測(cè)序質(zhì)量分析

cd ~/project/qc #切換到qc目錄
#比較經(jīng)trim_galore處理前后數(shù)據(jù)質(zhì)量
ls ~/project/fastq/*.gz | xargs fastqc -t 10 -o  ./raw
ls ~/project/clean/*.gz | xargs fastqc -t 10 -o ./trimed
#-t：調(diào)用核心數(shù)目
#-q：安靜運(yùn)行忿等，運(yùn)行過(guò)程中不會(huì)生成報(bào)告，在結(jié)束時(shí)將報(bào)告生成一個(gè)文件
#-o ./：文件輸出當(dāng)前目錄
#*.gz：輸入文件

查看QC結(jié)果
單個(gè)查看：鼠標(biāo)雙擊打開(kāi)html文件查看
批量查看：使用MultiQC軟件: multiqc *fastqc.zip

• 使用 MultiQC 將多個(gè) FastQC 的結(jié)果匯集到一個(gè)文件：

multiqc ../fastqc -n SRAxxx -o ../multiqcResults/
# multiqc 后面直接跟上 fastqc 結(jié)果所在的文件夾
# -n 輸出結(jié)果的文件名前綴
# -o 設(shè)置輸出結(jié)果的文件夾
#(有時(shí)間再看)

• Fastqc結(jié)果報(bào)告關(guān)注重點(diǎn)：

1）.basic statistics
2）.per base sequence quality
3）.per base sequcence content
4）.adaptor content
5）.sequence duplication levels

最主要的幾個(gè)指標(biāo)是

• GC含量
• Q20和Q30的比例以及是否存在接頭（adaptor）
• index以及其他物種序列的污染等崔挖。

從頁(yè)面左側(cè)的的【summary】中可以看出有哪些選項(xiàng)沒(méi)有通過(guò)这弧，可以看出此數(shù)據(jù)的測(cè)序質(zhì)量。從【Basic Statics】看出數(shù)據(jù)的序列數(shù)量虚汛，測(cè)序平臺(tái)以及GC含量等相關(guān)信息

圖片.png

Per base sequence quality : 每個(gè)read各位置堿基的測(cè)序質(zhì)量匾浪。
橫軸堿基的位置，縱軸是質(zhì)量分?jǐn)?shù)卷哩，
Quality score= -10log10p
（p代表錯(cuò)誤率）蛋辈，所以當(dāng)質(zhì)量分?jǐn)?shù)為40的時(shí)候，p就是0.0001将谊，質(zhì)量算高了冷溶。
[紅色線]代表中位數(shù)，
[藍(lán)色線]代表平均數(shù)
[黃色]是25%-75%區(qū)間尊浓，
[權(quán)]是10%-90%區(qū)間逞频。
若任一位置的下 [四分位數(shù)] 低于10或者 [中位數(shù)] 低于25---->出現(xiàn) “警告”
若任一位置的下 [四分位數(shù)] 低于5或者 [中位數(shù)] 低于20，出現(xiàn)“失敗栋齿，F(xiàn)ail”

圖片.png

此圖中

• 橫軸是測(cè)序序列第1個(gè)堿基到第51個(gè)堿基
• 縱軸是質(zhì)量得分苗胀，Q = -10*log10（error P）即20表示1%的錯(cuò)誤率，30表示0.1%
• 每1個(gè)boxplot瓦堵，都是該位置的所有序列的測(cè)序質(zhì)量的一個(gè)統(tǒng)計(jì):
  *上面的bar是90%分位數(shù)
  *下面的bar是10%分位數(shù)基协，
  *箱子的中間的橫線是50%分位數(shù)(上邊是75%分位數(shù)，下邊是25%分位數(shù))
• 圖中藍(lán)色的細(xì)線是各個(gè)位置的平均值的連線(一般要求此圖中菇用，所有位置的10%分位數(shù)大于20,即Q20過(guò)濾)
  澜驮。。上面的這個(gè)測(cè)序結(jié)果惋鸥，需要把后面的27bp以后的序列切除杂穷，從而保證后續(xù)分析的正確性
  【Failure 報(bào)錯(cuò) 如果任何堿基質(zhì)量低于5,或者是任何中位數(shù)低于20】

圖片.png

Per tile sequence quality 每tile的堿基質(zhì)量
【藍(lán)色】代表測(cè)序質(zhì)量很高，【暖色】代表測(cè)序質(zhì)量不高

圖片.png

剩下的自己百度

5. 使用bowtie2比對(duì)到基因組并排序

直接用bowtie2進(jìn)行比對(duì)和統(tǒng)計(jì)比對(duì)率, 需要提前下載參考基因組然后使用命令構(gòu)建索引卦绣，或者直接就下載索引文件：
Bowtie2 是將測(cè)序reads與長(zhǎng)參考序列比對(duì)工具耐量。適用于將長(zhǎng)度大約為50到100或1000字符的reads與相對(duì)較長(zhǎng)的基因組（如哺乳動(dòng)物）進(jìn)行比對(duì)。
Bowtie2 通常是比較基因組學(xué)管道的第一步迎卤，包括識(shí)別變體（variation calling）拴鸵，ChIP-seq，RNA-seq，BS-seq劲藐。Bowtie2 和Bowtie 也高度整合在一些工具中八堡，包括TopHat（快速拼接RNA-seq reads 的 mapper），Crossbow（重測(cè)序數(shù)據(jù)分析云的軟件工具）聘芜，Myrna（對(duì)齊RNA-seq reads和分析差異基因表達(dá)的云計(jì)算軟件工具）

5.1 下載官方索引

#下載小鼠參考基因組的索引和注釋文件, 這里用常用的mm10
# 索引大小為3.2GB兄渺， 不建議自己下載基因組構(gòu)建，可以直接下載索引文件汰现，代碼如下：
mkdir referece && cd reference
wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip #（官方索引）
unzip mm10.zip
#我的已下載保存在 ~/biosoft/referece 里
#人的h19/h38已構(gòu)建好存放在~/biosoft/referece

• 自建索引

#這里以構(gòu)建人,hg38為例 
#構(gòu)建參考基因組索引
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz 
tar -zxvf hg38.fa.gz  
cd ~/biosoft/reference/hg38/#在目標(biāo)文件夾構(gòu)建不用指定目錄
bowtie2-build --threads 10 -f hg38.fa hg38

索引構(gòu)建參考

5.2 bowtie2 使用方法

雙端
bowtie2 -p 10 -x genome_index -1 input_1.fq -2 input_2.fq | samtools sort -O bam -@ 10 -o - > output.bam

• 批量處理挂谍，運(yùn)行bowtie2 獲取 SAM 文件

cd ~/project/align
bowtie2_index=~/biosoft/reference/mm10/mm10#參考基因組索引路徑.
bin_bowtie2='/home/chen/biosoft/bowtie/bowtie2-2.2.9/bowtie2'
files_1=~/project/clean/*_1.fq.gz  #或files=cd ~/project/clean/;ls *.fq.gz
files_2=~/project/clean/*_2.fq.gz  

ls $files_1 | while read id
do 
fq1=$id
for i in `ls $files_2 | grep "${fq1%%_1_val*}"`
do 
fq2=$i
file=$(basename $fq1)#理解basename這個(gè)函數(shù)
fq=${file%%_1_val*}
done
echo $fq1 
echo $fq2
bowtie2 -p 12 -3 5 --local -x $bowtie2_index -1 $fq1 -2 $fq2 | samtools sort -O bam -@ 6 -o ~/project/align/${fq}.bam
done
#依賴已經(jīng)添加了，$bin_bowtie2可以直接換成bowtie2
#比對(duì)很慢
#-p/--threads NTHREADS
#-3/--trim3 <int>: Trim <int> bases from 3' (right) end of each read before alignment (default: 0).
#--local: 我的理解是不是精確比對(duì)瞎饲，允許個(gè)別堿基錯(cuò)配
#-x: The basename of the index for the reference genome. 
#bowtie2輸出格式默認(rèn)為SAM口叙，可以指定為多種其它格式，比如說(shuō)BAM（SAM的二進(jìn)制文件）等等
#利用samtools sort 進(jìn)行排序并將sam格式轉(zhuǎn)換為bam格式

• 也可以單獨(dú)將SAM 文件轉(zhuǎn)為 BAM 文件
  這個(gè) BAM 文件就是后面 MACS2 的輸入文件嗅战。
  samtools sort example.sam > example.bam

igvtools count -z 5 -w 25 *.bam *.bam.tdf hg38

將.bam文件轉(zhuǎn)換為.tdf文件妄田，導(dǎo)入IGV可查看peaks

單端
"bowtie2 -p 10 -x genome_index -U input.fq | samtools sort -O bam -@ 10 -o - > output.bam

cd ~/project/align
bowtie2_index=~/biosoft/reference/mm10/mm10 #指定參考基因組索引路徑.
ls ~/project/clean/*.fq.gz | while read id
do 
file=$(basename $id )
sample=${file%%.*}
echo $file $sample
## 比對(duì)過(guò)程3分鐘一個(gè)樣本
bowtie2 -p 5 --local -x  $bowtie2_index -U  $id | samtools sort  -O bam  -@ 5 -o ~/project/align/${sample}.bam
done 

需要注意的是：genome_index 指的是用于bowtie2的索引文件（如下圖），而不是參考基因組本身驮捍，構(gòu)建過(guò)程參考前文5.1部分疟呐。genome_index 需要指定路徑及其共用文件名，比如我的索引文件放在~/biosoft/referece/mm10目錄下东且，但是需要輸入的參數(shù)為~/biosoft/referece/mm10/mm10启具。最后一個(gè)mm10指的是共用文件名。

#例如我的索引文件
chen@chen:~/biosoft/referece/mm10$ ls -l
總用量 3678284
-rw-r--r-- 1 chen chen 888464705 5月   2  2012 mm10.1.bt2
-rw-r--r-- 1 chen chen 663195880 5月   2  2012 mm10.2.bt2
-rw-r--r-- 1 chen chen      6119 5月   2  2012 mm10.3.bt2
-rw-r--r-- 1 chen chen 663195875 5月   2  2012 mm10.4.bt2
-rw-r--r-- 1 chen chen 888464705 5月   3  2012 mm10.rev.1.bt2
-rw-r--r-- 1 chen chen 663195880 5月   3  2012 mm10.rev.2.bt2

bowtie2必須參數(shù)

bowtie2可選參數(shù)

5.3對(duì)bam文件進(jìn)行QC

cd ~/project/epi/align
ls  *.bam  | xargs -i samtools index {} #生成.bai文件不知如何看
ls  *.bam  | while read id ;do (nohup samtools flagstat $id > $(basename $id ".bam").stat & );done
#生成 .stat文件

5.4合并bam文件

參考這里
因?yàn)橐粋€(gè)樣品分成了多個(gè)lane進(jìn)行測(cè)序珊泳，所以在進(jìn)行peaks calling的時(shí)候鲁冯，需要把bam進(jìn)行合并,即一個(gè)樣分成了多個(gè)SRR...數(shù)據(jù),需要合并。

## 如果不用循環(huán)：
## samtools merge  control.merge.bam Control_1_trimmed.bam Control_2_trimmed.bam
## 通常我們用批處理旨椒。
cd ~/project/
mkdir mergeBam
cd ~/project/align
ls *.bam | sed 's/_[0-9]_trimmed.bam//g' | sort -u | while read id;do samtools merge ../mergeBam/$id.merge.bam $id*.bam ;done

5.5去除PCR重復(fù)

5.5.1 Sambamba去除PCR重復(fù)（推薦此方法）

cd ~/project/rmdup
source activate chipseq
files_bam=~/project/align/*.bam
files_name=$(basename -s .bam $files_bam)#-s表示去除后綴
ls $files_bam | while read id 
do
echo $id
sambamba markdup -r -p -t 12 $id $(basename -s .bam $id).rmdup.bam
done
#-t是線程數(shù)晓褪，-p顯示過(guò)程，-r remove duplicates instead of just marking them
#處理后同時(shí)會(huì)給出索引文件

5.5.2 samtools markdup去除PCR重復(fù)（不推薦综慎，還需要多一步才能進(jìn)行）

cd ~/project/rmdup
source activate chipseq
files_bam=~/project/align/*.bam
ls $files_bam  | while read id ;do  samtools markdup -r $id $(basename $id ".bam").rmdup.bam ;done#去除PCR重復(fù), -r是直接去除重復(fù), 不加是直接標(biāo)記
ls $files_bam  | xargs -i samtools index {} 

5.6計(jì)算比對(duì)率

#以比對(duì)后的bam文件進(jìn)行計(jì)算
cd ~/project/align/
files_bam=~/project/align/*.bam
ls $files_bam  | while read id 
do 
samtools flagstat $id > $(basename -s .bam $id).stat
done
#獲取比對(duì)率

6 使用macs2進(jìn)行peaks calling

macs2包含一系列的子命令，其中最主要的就是callpeak勤庐， 官方提供了使用實(shí)例

macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01

• -t: 實(shí)驗(yàn)組的輸出結(jié)果- -c: 對(duì)照組的輸出結(jié)果
• -f: -t和-c提供文件的格式示惊，可以是”ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT”, “ELANDMULTIPET” (for pair-end tags), “SAM”, “BAM”, “BOWTIE”, “BAMPE” “BEDPE” 任意一個(gè)。如果不提供這項(xiàng)愉镰，就是自動(dòng)檢測(cè)選擇米罚。
• -g: 基因組大小， 默認(rèn)提供了hs, mm, ce, dm選項(xiàng)丈探， 不在其中的話录择，比如說(shuō)擬南芥，就需要自己提供了。
• -n: 輸出文件的前綴名
• -B: 會(huì)保存更多的信息在bedGraph文件中隘竭，如fragment pileup, control lambda, -log10pvalue and -log10qvalue scores
• -q: q值塘秦，也就是最小的PDR閾值， 默認(rèn)是0.05动看。q值是根據(jù)p值利用BH計(jì)算尊剔，也就是多重試驗(yàn)矯正后的結(jié)果。
• -p： 這個(gè)是p值菱皆，指定p值后MACS2就不會(huì)用q值了须误。
• -m: 和MFOLD有關(guān)，而MFOLD和MACS預(yù)構(gòu)建模型有關(guān)仇轻，默認(rèn)是5：50京痢，MACS會(huì)先尋找100多個(gè)peak區(qū)構(gòu)建模型，一般不用改篷店。

cd ~/project/peaks/  
files_bams=~/project/rmdup/*.bam  #如果沒(méi)有去除PCR重復(fù)就用files_bams=~/project/align/*.bam

#control.bam換成自己的對(duì)照組
ls $files_bams | while read id
do
sample=$(basename -s .bam $id)
echo $sample
macs2 callpeak -c ~/project/rmdup/control.bam -t $id -f BAM --bdg -g mm -n $sample.bam --outdir ~/project/peaks/
done  
#雙端測(cè)序結(jié)果用BAMPE??
#control.bam是對(duì)照組
#輸出文件.bdg可以直接導(dǎo)入IGV可視化
#為何指定了輸出文件夾历造，但還是輸出在.bam所在的文件夾里？
#沒(méi)有control文件也可以不用加船庇，系統(tǒng)會(huì)用算法模擬一個(gè)吭产，-t 后可接兩個(gè)bam文件，因?yàn)檎綄?shí)驗(yàn)一般測(cè)兩個(gè)重復(fù)



#或者用下面這個(gè)循環(huán)
for i in files_bams
do 
macs2 callpeak -c control.bam -t $i -n $i -f BAM -g hs --outdir ~/project/peaks --bdg -q 0.05
done

#沒(méi)有control數(shù)據(jù)如下處理
for i in files_bams
do 
macs2 callpeak -f BAM -t $i -n $i -g hs --outdir ~/project/peaks --bdg -q 0.05
done
# -f 指定輸入文件的格式鸭轮，如SAM臣淤、BAM、BED等
# -c 對(duì)照組窃爷，可以接多個(gè)數(shù)據(jù)邑蒋，空格分隔。
# -t 實(shí)驗(yàn)組按厘，ChIP-seq數(shù)據(jù)医吊。可以接多個(gè)數(shù)據(jù)逮京，空格分隔卿堂。
# -n 輸出文件名的前綴
# -g 有效基因組大小。軟件有幾個(gè)默認(rèn)值懒棉，如hs指human草描，mm指mouse
# --outdir 輸出結(jié)果的存儲(chǔ)路徑
# --bdg 輸出 bedgraph 格式的文件
# -q FDR閾值, 默認(rèn)為0.05
# NAME 是使用 -n 設(shè)置的輸出文件前綴
# NAME_control_lambda.bdg 和 NAME_treat_pileup.bdg 
## bdg即為bedGraph格式，可以導(dǎo)入U(xiǎn)CSC或者轉(zhuǎn)換為bigwig格式策严。
# 兩個(gè)bdg文件：
  #- treat_pileup
  #- control_lambda
#2. NAME_model.r           # R代碼
#3. NAME_peaks.narrowPeak  # BED 6+4 格式
#4. NAME_peaks.xls         # 包含 peak 信息的 xls 文件穗慕，注意 chrStart 是1-based
#5. NAME_peaks_summits.bed 
# BED 格式，包含peak的summits位置妻导，第5列是-log10pvalue逛绵。如果想找motif怀各，推薦使用此文件。

NAME 是使用 -n 設(shè)置的輸出文件前綴

1. NAME_control_lambda.bdg 和 NAME_treat_pileup.bdg
   bdg即為bedGraph格式术浪，可以導(dǎo)入U(xiǎn)CSC或者轉(zhuǎn)換為bigwig格式瓢对。
2. NAME_model.r # R代碼
3. NAME_peaks.narrowPeak # BED 6+4 格式
4. NAME_peaks.xls # 包含 peak 信息的 xls 文件，注意 chrStart 是1-based
5. NAME_peaks_summits.bed
   BED 格式添吗，包含peak的summits位置沥曹，第5列是-log10pvalue。
   如果想找motif碟联，推薦使用此文件妓美。

其中需要注意，NAME_model.r 以 R 代碼的形式保存了“雙峰模型”鲤孵。

在命令行中輸入：Rscript NAME_model.r會(huì)得到名為 NAME_model.pdf 的文件↓：

圖片.png

首先我們要記住一點(diǎn)——所有數(shù)據(jù)分析的流程都基于實(shí)驗(yàn)原理贵试。

做 ChIP-seq，就是為了確認(rèn)蛋白-DNA結(jié)合位點(diǎn)凯正。測(cè)序獲得的 read 只是跟隨著蛋白一起沉淀下來(lái)的 DNA 片段的末端毙玻，并不是真實(shí)的 DNA 和蛋白結(jié)合的位置。怎么找到蛋白和DNA真正的結(jié)合位置廊散？
單端測(cè)序桑滩，對(duì)于測(cè)序的同一片段，reads 在正負(fù)鏈上的分布是平均的允睹，并且正負(fù)鏈的 reads 之間會(huì)有一個(gè)間距运准。近似的形成“雙峰”。因此缭受，把來(lái)自于同一片段位于不同鏈上的 reads 相向移動(dòng)胁澳，就能有效地降低信號(hào)中的噪音，找到結(jié)合的位置米者。而對(duì) 雙端測(cè)序 而言韭畸，它本身測(cè)的就是文庫(kù)的兩端，因此不用建立模型和偏倚塘雳。我們只需要對(duì) MACS 設(shè)置參數(shù) --nodel 就能略過(guò)雙峰模型建立的過(guò)程.

7 Peak 可視化

7.1 bam文件轉(zhuǎn)化為bw(bigWig)或bedgraph文件

BAM轉(zhuǎn)換為bigWig或bedGraph
BAM文件是SAM的二進(jìn)制轉(zhuǎn)換版陆盘，應(yīng)該都知道。那么bigWig格式是什么败明？bigWig是wig或bedGraph的二進(jìn)制版，存放區(qū)間的坐標(biāo)軸信息和相關(guān)計(jì)分(score)太防，主要用于在基因組瀏覽器上查看數(shù)據(jù)的連續(xù)密度圖妻顶，可用wigToBigWig從wiggle進(jìn)行轉(zhuǎn)換酸员。
bedGraph和wig格式是什么? USCS的幫助文檔稱這兩個(gè)格式數(shù)是已經(jīng)過(guò)時(shí)的基因組瀏覽器圖形軌展示格式，前者展示稀松型數(shù)據(jù)讳嘱，后者展示連續(xù)性數(shù)據(jù)幔嗦。目前推薦使用bigWig/bigBed這兩種格式取代前兩者。
為什么要用bigWig? 主要是因?yàn)锽AM文件比較大沥潭，直接用于展示時(shí)對(duì)服務(wù)器要求較大邀泉。因此在GEO上僅會(huì)提供bw,即bigWig下載，便于下載和查看钝鸽。如果真的感興趣汇恤，則可以下載原始數(shù)據(jù)進(jìn)行后續(xù)分析。

deeptools

使用deeptool進(jìn)行可視化
deeptools提供bamCoverage和bamCompare進(jìn)行格式轉(zhuǎn)換拔恰，為了能夠比較不同的樣本因谎，需要對(duì)先將基因組分成等寬分箱(bin)，統(tǒng)計(jì)每個(gè)分箱的read數(shù)颜懊，最后得到描述性統(tǒng)計(jì)值。對(duì)于兩個(gè)樣本，描述性統(tǒng)計(jì)值可以是兩個(gè)樣本的比率锋勺，或是比率的log2值沙庐，或者是差值。如果是單個(gè)樣本咸这，可以用SES方法進(jìn)行標(biāo)準(zhǔn)化夷恍。
bamCoverage的基本用法

source activate chipseq
bamCoverage -e 170 -bs 10 -b sample.bam -o sample.bw
#sample.bam是前期比對(duì)得到的BAM文件

得到的bw文件就可以送去IGV/Jbrowse進(jìn)行可視化。 這里的參數(shù)僅使用了-e/--extendReads和-bs/--binSize即拓展了原來(lái)的read長(zhǎng)度炊苫，且設(shè)置分箱的大小裁厅。其他參數(shù)還有--filterRNAstrand {forward, reverse}: 僅統(tǒng)計(jì)指定正鏈或負(fù)鏈
--region/-r CHR:START:END: 選取某個(gè)區(qū)域統(tǒng)計(jì)--smoothLength: 通過(guò)使用分箱附近的read對(duì)分箱進(jìn)行平滑化, 如果為了其他結(jié)果進(jìn)行比較，還需要進(jìn)行標(biāo)準(zhǔn)化侨艾，deeptools提供了如下參數(shù)：

--scaleFactor: 縮放系數(shù)
--normalizeUsingRPKMReads: Per Kilobase per Million mapped reads (RPKM)標(biāo)準(zhǔn)化
--normalizeTo1x: 按照1x測(cè)序深度(reads per genome coverage, RPGC)進(jìn)行標(biāo)準(zhǔn)化
--ignoreForNormalization： 指定那些染色體不需要經(jīng)過(guò)標(biāo)準(zhǔn)化
如果需要以100為分箱执虹，并且標(biāo)準(zhǔn)化到1x，且僅統(tǒng)計(jì)某一條染色體區(qū)域的正鏈唠梨，輸出格式為bedgraph,那么命令行可以這樣寫(xiě)

• bam轉(zhuǎn)換為bedgraph文件

cd ~/project/bedgraph/
source activate chipseq
sample=~/project/rmdup/*.bam
ls $sample | while read id
do
echo $id
sample_name=$(basename -s .bam $id)
bamCoverage -bs 100 --normalizeUsing CPM -of bedgraph -b $id -o $sample_name.bedgraph
done 

bamCompare和bamCoverage類似袋励，只不過(guò)需要提供兩個(gè)樣本，并且采用SES方法進(jìn)行標(biāo)準(zhǔn)化当叭，于是多了--ratio參數(shù)茬故。

• 批量把bam文件轉(zhuǎn)為bw文件，詳情：http://www.bio-info-trainee.com/1815.html

cd ~/project/bw/
source activate chipseq
sample=~/project/rmdup/*.bam
ls $sample | xargs -i samtools index {} #轉(zhuǎn)換前需要先構(gòu)建index文件以.bai為后綴
#其實(shí)在Sambamba去除PCR重復(fù)時(shí)已經(jīng)自動(dòng)構(gòu)建好了蚁鳖，存放在rmdup目錄下
ls $sample | while read id
do
echo $id
sample_name=$(basename -s .bam $id)
bamCoverage --normalizeUsing CPM -b $id -o $sample_name.bw
done 

轉(zhuǎn)換的很慢
參數(shù)：--outFileName, -o Output file name.
--outFileFormat, -of Possible choices: bigwig, bedgraph. Output file type: Either “bigwig” or “bedgraph”.
具體可參考這里

7.2 IGV可視化

• 使用 IGV 查看 BedGraph 文件
  bw文件就可以送去IGV/Jbrowse進(jìn)行可視化
  進(jìn)入IGV目錄磺芭，運(yùn)行IGV

cd ~/biosoft/igv/IGV_Linux_2.9.4/
bash igv.sh

bash igv.sh即可打開(kāi)桌面版IGV軟件。
這里直接上傳 bg/bw 文件：File->Load from File->選擇sample.bg并打開(kāi)醉箕。

7.3 bam文件轉(zhuǎn)化為bed文件

#將bam文件轉(zhuǎn)換為bed文件
bedtools bamtobed -i ~/project/align/SRR065175.bam  > ~/project/bed/SRR065175.bed

7.4 查看轉(zhuǎn)錄起始點(diǎn)TSS附近信號(hào)強(qiáng)度

deeptools軟件里的computeMatrix命令可以查看TSS附近peaks強(qiáng)度钾腺，輸入文件為 .bw徙垫，并給一個(gè)參考的注釋文件（即 .bed文件下載方法見(jiàn)這里及這里我的 .bed文件已下載存儲(chǔ)在~/biosoft/refseq_bed目錄里。三個(gè)文件分別是hg19.bed, hg38.bed, mm10.bed）放棒，從而讓軟件查看樣本在TSS（轉(zhuǎn)錄起始位點(diǎn)）附近是否有富集姻报。

• 首先對(duì)單一樣本繪圖：

cd ~/project/tss 
##both -R and -S can accept multiple files 
sample=~/project/bw/*.bw
ls $sample | while read id
do
echo $(basename $id)
computeMatrix reference-point  --referencePoint TSS -p 15 -b 2000 
-a 2000 -R ~/biosoft/refseq_bed/mm10.bed 
-S $id --skipZeros  
-o matrix1_test_TSS.gz
--outFileSortedRegions regions1_test_genes.bed
done
#輸出文件在當(dāng)前tss目錄里 輸出2個(gè)文件，-b和-a 是TSS上下游距離间螟，-S是輸入上文中bam轉(zhuǎn)換成的.bw文件

##both plotHeatmap and plotProfile will use the output from   computeMatrix
plotHeatmap -m matrix1_test_TSS.gz  -out test_Heatmap.png
plotHeatmap -m matrix1_test_TSS.gz  -out test_Heatmap.pdf --plotFileFormat pdf  --dpi 720  
plotProfile -m matrix1_test_TSS.gz  -out test_Profile.png
plotProfile -m matrix1_test_TSS.gz  -out test_Profile.pdf --plotFileFormat pdf --perGroup --dpi 720 

• (批量處理)畫(huà)10K附近吴旋，即TSS前后10K范圍，可通過(guò)修改參數(shù)-a和-b來(lái)修改厢破。

cd ~/project/tss 
bed=~/biosoft/refseq_bed/mm10.bed
for id in ~/project/bw/*.bw
do 
echo $id
file=$(basename $id)
sample=${file%%.*} 
echo $sample
computeMatrix reference-point -S $id -R $bed  -p 15 -a 2000 -b 2000
--referencePoint TSS --skipZeros -o ~/project/tss/matrix1_${sample}_TSS.gz
--outFileSortedRegions regions1_${sample}_TSS.bed
done
#輸出的gz文件用于plotHeatmap, plotProfile命令作圖荣瑟。
#輸出文件在當(dāng)前tss目錄里 輸出2個(gè)文件，-b和-a 是TSS上下游距離溉奕，-S是輸入上文中bam轉(zhuǎn)換成的.bw文件

##both plotHeatmap and plotProfile will use the output from   computeMatrix
for i in ~/project/tss/*.gz
do
file=$(basename $i)
sample=${file%%.*}
echo $sample
plotHeatmap -m ${sample}.gz  -out ${sample}_Heatmap_10K.png
plotHeatmap -m ${sample}.gz  -out ${sample}_Heatmap_10K.pdf --plotFileFormat pdf  --dpi 720  
plotProfile -m ${sample}.gz  -out ${sample}_Profile_10K.png
plotProfile -m ${sample}.gz  -out ${sample}_Profile_10K.pdf --plotFileFormat pdf --perGroup --dpi 720 
done 

用plotHeatmap以熱圖的方式對(duì)覆蓋進(jìn)行可視化褂傀，用plotProfile以折線圖的方式展示覆蓋情況。
computeMatrix具有兩個(gè)模式:scale-region和reference-point加勤。前者用來(lái)信號(hào)在一個(gè)區(qū)域內(nèi)分布仙辟，后者查看信號(hào)相對(duì)于某一個(gè)點(diǎn)的分布情況。無(wú)論是那個(gè)模式鳄梅，都有有兩個(gè)參數(shù)是必須的叠国，-S是提供bigwig文件，-R是提供基因的注釋信息戴尸。還有更多個(gè)性化的可視化選項(xiàng)粟焊。

7.5 使用R包對(duì)找到的peaks文件進(jìn)行注釋

bedPeaksFile         = '8WG16_summits.bed'; 
bedPeaksFile
## loading packages
require(ChIPseeker)
require(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
require(clusterProfiler) 
peak <- readPeakFile( bedPeaksFile )  
keepChr= !grepl('_',seqlevels(peak))
seqlevels(peak, pruning.mode="coarse") <- seqlevels(peak)[keepChr]
peakAnno <- annotatePeak(peak, tssRegion=c(-3000, 3000), 
                         TxDb=txdb, annoDb="org.Mm.eg.db") 
peakAnno_df <- as.data.frame(peakAnno)

可以載入IGV看看效果，檢測(cè)軟件找到的peaks是否真的合理孙蒙，還可以配合rmarkdown來(lái)出自動(dòng)化報(bào)告项棠。
也可以使用其它代碼進(jìn)行下游分析； https://github.com/jmzeng1314/NGS-pipeline/tree/master/CHIPseq

matrix1_BMSC_d00_h3k27me3_R2_trimmed_TSS_2K_Heatmap_10K.png

7.6 peaks相關(guān)基因集的注釋

都是得到感興趣基因集挎峦，然后注釋香追，分析方法等同于GEO數(shù)據(jù)挖掘課程或者轉(zhuǎn)錄組下游分析： https://github.com/jmzeng1314/GEO （有配套視頻）

8 用homer軟件尋找motif

conda install -c bioconda homer ,homer安裝在 ~/anaconda3/envs/chipseq/share/homer下，然后使用其附帶的配置腳本下載人和小鼠的基因組數(shù)據(jù)庫(kù)（5G左右比較大,我的已經(jīng)下載好存放在~/anaconda3/envs/chipseq/share/homer/data/genomes/）坦胶。

perl ~/anaconda3/envs/chipseq/share/homer/configureHomer.pl -install hg19#優(yōu)先使用這行命令下載
#上面安裝不了試試這行命令：perl ~/anaconda3/envs/chipseq/share/homer/configureHomer.pl  -install mm10 
## 我們上游分析是基于mm10找到的peaks文件

## 下載成功后會(huì)多出 ~/anaconda3/envs/chipseq/share/homer/data/genomes/mm10/ 文件夾, 共 4.9G
## 這個(gè)文件夾取決于你把homer這個(gè)軟件安裝到了什么地方透典。

## 或者用下面代碼安裝:這種安裝會(huì)裝在~/biosoft/homer文件夾下，運(yùn)行會(huì)報(bào)錯(cuò)顿苇，把/genomes/mm10/ 和config文件復(fù)制粘貼到anaconda3文件夾下的homer即可峭咒。
cd ~/biosoft/homer
wget http://homer.salk.edu/homer/configureHomer.pl 
perl configureHomer.pl -install
perl configureHomer.pl -install hg19

homer軟件找motif整合了兩個(gè)方法，包括依賴于數(shù)據(jù)庫(kù)的查詢纪岁，和de novo的推斷,都是讀取ChIP-seq數(shù)據(jù)上游分析得到的bed格式的peaks文件凑队。

運(yùn)行homer軟件需要輸入 .bed文件
但是使用起來(lái)很簡(jiǎn)單：http://homer.ucsd.edu/homer/ngs/peakMotifs.html
findMotifsGenome.pl <peak/BED file> <genome> <output directory> -size # [options]
例如

findMotifsGenome.pl ERpeaks.txt hg18 ER_MotifOutput/ -size 200 -mask

• homer motif及peaks注釋，具體方法如下：

cd  ~/project/motif  
files=~/project/peaks/*.bed
ls $files | while read id
do
echo $id
file=$(basename -s .bed $id)
echo $file
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' $id > homer_peaks.tmp
findMotifsGenome.pl homer_peaks.tmp mm10 ~/project/motif/${file}  #自動(dòng)創(chuàng)建目錄
sample=$(basename -s trimmed.rmdup.bam_summits.bed $id)
echo ${sample}_Annotation
annotatePeaks.pl $id mm10 >~/project/annotation/${sample}.peakAnn.xls 2>~/project/annotation/${sample}.annLog.txt   #注釋也一起做了
#homer也可以peaks注釋
done

#或者用下面這行
for id in ~/chen/project/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幔翰，然后運(yùn)行：nohup bash runMotif.sh 1>motif.log &
不僅僅找了motif顽决，還順便把peaks注釋了一下短条。得到的后綴為peakAnn.xls 的文件就可以看到和使用R包注釋的結(jié)果是差不多的导匣。
結(jié)果分析參考這里
還可以使用meme來(lái)找motif才菠，需要通過(guò)bed格式的peaks的坐標(biāo)來(lái)獲取fasta序列。MEME贡定，鏈接：http://meme-suite.org/

其它高級(jí)分析

比如可以 比較不同的peaks文件赋访，代碼見(jiàn)：https://github.com/jmzeng1314/NGS-pipeline/blob/master/CHIPseq/step6-ChIPpeakAnno-Venn.R

• 本文參考自
  https://zhuanlan.zhihu.com/p/90180058
  http://www.reibang.com/p/1384173c353b
  腳本文件在：https://github.com/jmzeng1314/NGS-pipeline/tree/master/CHIPseq