1. 环境配置

#!/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先进行检索
conda search trim_galore
## 保证所有的软件都是安装在 chipseq 这个环境下面
conda install -y sra-tools  
conda install -y trim-galore  samtools
conda install -y deeptools homer  meme
conda install -y macs2 bowtie bowtie2 

2. 数据下载

(1)prefetch下载

首先找到要下载的数据（SRR...号）
创建project文件夹

cd project/
#将下列数据写入vim文件中
vim SraAccList.txt
SRR1266976
SRR1266977
SRR1266978
SRR1266979
SRR1266980
SRR1266981

利用chipseq环境中prefetch 软件进行下载
如果不指定下载目录 则默认下载目录：~/ncbi/public/sra/

source activate chipseq
mkdir {sra,bedgraph,fastq,rmdup,tss,clean,align,peaks,motif,qc/{raw,trimed},annotation} #在project目录下新建目录
#将数据下载于sra目录下
cat SraAccList.txt | while read id;
do
(nohup prefetch -O ./sra $id &)#保存在sra目录下
done
#control +c可以后台下载

prefetch --option-file SraAccList.txt #或者批量下载可以点击下载SraAccList.txt文件

prefetch常用选项

(2)Aspera Connect下载

Aspera—ascp下载SRA数据

SRA数据库下载：数据的存放地址是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官网搜索srr然后下载获取SRR_Acc_List.txt 然后观察是err+7位数的，所以直接使用7位数的代码。（7位数y取最后一位，8位数y取最后2位，6位数把y变量删掉，应该是这样。）一般来说，NCBI的sra文件前面的地址都是一样的 ~/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk，那么写脚本批量下载也就不难了！最后那个 ./ 表示下载后的路径。

1. 1批量下载srr文件（推荐）

cd ~/
mkdir -p project/{sra,bedgraph,fastq,rmdup,tss,clean,align,peaks,motif,qc/{raw,trimed},annotation}#创建文件夹

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 单个下载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数据是2个文件 *_1.fastq 和 *_2.fastq表明是双端测序。

Aspera用法如下：
ascp [参数] 目标文件 保存路径
-v verbose mode 实时知道程序在干啥
-T 取消加密，否则有时候数据下载不了
-i 提供私钥文件的地址
-l 设置最大传输速度，一般200m到500m，如果不设置，反而速度会比较低，可能有个较低的默认值
-k 断点续传，一般设置为值1
-Q 一般加上它
-P 提供SSH port，端口一般是33001

3. sra文件转fastq文件：

• sra转化为fastq文件可以使用sratoolkit中的fastq-dump命令。

fastq-dump --split-3 filename

其中 --split-3 参数代表着如果是单端测序就生成一个 、.fastq文件，如果是双端测序就生成*_1.fastq 和*_2.fastq 文件；即单端测序输出1个fastq文件，双端测序输出2个fastq文件。
进入到sra文件中,将SRR文件夹去掉，数据文件直接保存在sra目录中，我们可以用下述代码进行批量的格式转化：
首先制作如下对应的文本

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

然后

## 下面需要用循环
cd ~/project/sra/
#把list.txt文件复制到sra目录下
source activate chipseq
## 下面用到的 list.txt 文件，就是上面自行制作的对应的文件。
cat list.txt | while read id
do 
echo $id
arr=($id)
srr=${arr[0]} #通过中括号获取数组中的某一个元素：${arr[0]} 得到第一个元素，${arr[0]} 第二个
sample=${arr[1]}
fastq-dump -A $sample -O ~/project/fastq/ --gzip --split-3  $srr
done 
#-A -$sample是输出文件名字
 #--gzip打包节省空间，且对后续分析不影响，如果有影响就把这个参数去掉，-O输出到目录

4. 使用trim_galore软件对fastq进行质量控制-过滤

（1）基础知识

trim_galore使用手册

cutadapt 软件可以对NGS数据进行质量过滤
FastQC 软件可以查看NGS数据的质量分布
trim_galore 将这两个软件封装到一起，使用起来更加方便。

• Trim galore，是可以自动检测adapter
• 去除reads 3’端的低质量碱基 # 自动调用cutadapt
• 去除adapter序列 # 自动去除3’端的adapter,可以通过设定--Illumina，-- small_rna，--nextera参数来指定对应的adapter类型
• 去除长度太短的序列 #通过设定--length参数，小于设定值被去除
• 其它过滤

参考：https://cloud.tencent.com/developer/article/1625193
和https://www.jianshu.com/p/7a3de6b8e503?from=groupmessage
这两个链接要看。

说明

1. Trim Galore是对FastQC和Cutadapt的包装。
2. trim_galore适用于所有高通量测序，包括RRBS(Reduced Representation Bisulfite-Seq ), Illumina、Nextera和smallRNA测序平台的双端和单端数据。
3. 主要功能包括两步：
   第一步 首先去除低质量碱基，然后去除3' 末端的adapter(如果没有指定具体的adapter，程序会自动检测前1million的序列）
   第二步 对比前12-13bp的序列是否符合以下类型的adapter :

• Illumina: AGATCGGAAGAGC # 如果输入参数 --Illumina,就会默认trimmed前13bp的adapter
• Small RNA: TGGAATTCTCGG # 同上 如果输入参数 --small_rna,就会默认trimmed前12bp的adapter
• Nextera: CTGTCTCTTATA # 同上 如果输入参数 --nextera,就会默认trimmed前12bp的adapter

对于单端测序数据，基本用法如下
trim_galore --quality 20 -a AGATCGGAAGAGC --length 20 -o out_dir input.fq
其中-a后面可以跟着序列（-aAGATCGGAAGAGC）

对于双端测序数据，基本用法如下
trim_galore --paired --quality 20 -a AGATCGGAAGAGC -a2 AGATCGGAAGAGC --length 20 -o out_dir R1.fq.gz R2.fq.gz

参数说明

--quality/-q<INT>：设定Phred quality score阈值，默认为Phred 20 切除质量得分低于设定值的序列。
--phred33/64：使用ASCII+33/64质量得分作为Phred得分，选择-phred33或者-phred64，表示测序平台使用的Phred quality score。具体怎么选择，看你用什么测序平台。
（需要确认:Sanger/Illumina 1.9+ encoding为33; Illumina 1.5 encoding为64）
--adapter/-a ：输入adapter序列。也可以不输入，Trim Galore!会自动寻找可能性最高的平台对应的adapter。自动搜选的平台三个，也直接显式输入这三种平台，即--illumina、--nextera和--small_rna。
-s/--stringency<INT>：设定可以忍受的前后adapter重叠的碱基数，默认为1（非常苛刻）。可以适度放宽，因为后一个adapter几乎不可能被测序仪读到。
--length<INT>：设定输出reads长度阈值小于设定值会被抛弃，默认值20bp；即小于20bp的被去除。注意，在pe150下，不要设置太高，可以50或36（默认20）
--max_length : 设置长度大于此值被丢弃
-e <ERROR RATE> ：默认0.1
--paired：对于双端测序结果，一对reads中，如果有一个被剔除，那么另一个会被同样抛弃，而不管是否达到标准。
--retain_unpaired：对于双端测序结果，一对reads中，如果一个read达到标准，但是对应的另一个要被抛弃，达到标准的read会被单独保存为一个文件。
--gzip和--dont_gzip：清洗后的数据zip打包或者不打包。
-o/--output_dir：输入目录 [需要提前建立目录，否则运行会报错]。
-- trim-n : 移除read一端的reads
-j :使用线程数, 注意假设已使用Python 3并安装了Pigz，那么内核设置为4，实际使用内核是15，因此，最高设为4.
--fastqc #当分析结束后，使用默认选项对结果文件进行fastqc分析

phred 33 和phred 64对照表

1. 如何判断phred 33 和phred 64?
   有时候得到的原始fastq文件，无法知道质量值体系，你就无法进行质量值的过滤，我们可以在正常情况下，按照上面的表对回去，统计一下几条reads的最大和最小质量值的区间，就可以知道到底是phred 33 还是phred 64体系。
2. 根据上图中Phred+33与Phred+64所使用的质量字符范围的不同，可以对fastq文件中质量得分的编码方式做出判断(第四行是测序质量，查看第四行符号是在Phred+64还是Phred+33范围)。图中显示，ASCII值小于等于58（相应的质量得分小于等于25）对应的字符只有在Phred+33的编码中被使用，所有Phred+64所使用的字符的ASCII值都大于等于59。在通常情况下，ASCII值大于等于74的字符只出现在Phred+64中。利用这些信息即可在程序中进行判断。（近几年应该都是Phred+33）

（2）数据处理

• 处理单个双端测序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/

• 批量处理双端测序数据

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

单端测序结果处理

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 
#具体稍作调整

双端测序结果主文件为 *_1_val_1.fq.gz 和 *_2_val_2.fq.gz，单端测序结果为 *_trimmed.fq。

(3)trim_galore处理结束，利用 --fastqc进行测序质量分析

cd ~/project/qc #切换到qc目录
#比较经trim_galore处理前后数据质量
ls ~/project/fastq/*.gz | xargs fastqc -t 10 -o  ./raw
ls ~/project/clean/*.gz | xargs fastqc -t 10 -o ./trimed
#-t：调用核心数目
#-q：安静运行，运行过程中不会生成报告，在结束时将报告生成一个文件
#-o ./：文件输出当前目录
#*.gz：输入文件

查看QC结果
单个查看：鼠标双击打开html文件查看
批量查看：使用MultiQC软件: multiqc *fastqc.zip

• 使用 MultiQC 将多个 FastQC 的结果汇集到一个文件：

multiqc ../fastqc -n SRAxxx -o ../multiqcResults/
# multiqc 后面直接跟上 fastqc 结果所在的文件夹
# -n 输出结果的文件名前缀
# -o 设置输出结果的文件夹
#(有时间再看)

• Fastqc结果报告关注重点：

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

最主要的几个指标是

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

从页面左侧的的【summary】中可以看出有哪些选项没有通过，可以看出此数据的测序质量。从【Basic Statics】看出数据的序列数量，测序平台以及GC含量等相关信息

图片.png

Per base sequence quality : 每个read各位置碱基的测序质量。
横轴碱基的位置，纵轴是质量分数，
Quality score= -10log10p
（p代表错误率），所以当质量分数为40的时候，p就是0.0001，质量算高了。
[红色线]代表中位数，
[蓝色线]代表平均数
[黄色]是25%-75%区间，
[权]是10%-90%区间。
若任一位置的下 [四分位数] 低于10或者 [中位数] 低于25---->出现 “警告”
若任一位置的下 [四分位数] 低于5或者 [中位数] 低于20，出现“失败，Fail”

图片.png

此图中

• 横轴是测序序列第1个碱基到第51个碱基
• 纵轴是质量得分，Q = -10*log10（error P）即20表示1%的错误率，30表示0.1%
• 每1个boxplot，都是该位置的所有序列的测序质量的一个统计:
  *上面的bar是90%分位数
  *下面的bar是10%分位数，
  *箱子的中间的横线是50%分位数(上边是75%分位数，下边是25%分位数)
• 图中蓝色的细线是各个位置的平均值的连线(一般要求此图中，所有位置的10%分位数大于20,即Q20过滤)
  。。上面的这个测序结果，需要把后面的27bp以后的序列切除，从而保证后续分析的正确性
  【Failure 报错 如果任何碱基质量低于5,或者是任何中位数低于20】

图片.png

Per tile sequence quality 每tile的碱基质量
【蓝色】代表测序质量很高，【暖色】代表测序质量不高

图片.png

剩下的自己百度

5. 使用bowtie2比对到基因组并排序

直接用bowtie2进行比对和统计比对率, 需要提前下载参考基因组然后使用命令构建索引，或者直接就下载索引文件：
Bowtie2 是将测序reads与长参考序列比对工具。适用于将长度大约为50到100或1000字符的reads与相对较长的基因组（如哺乳动物）进行比对。
Bowtie2 通常是比较基因组学管道的第一步，包括识别变体（variation calling），ChIP-seq，RNA-seq，BS-seq。Bowtie2 和Bowtie 也高度整合在一些工具中，包括TopHat（快速拼接RNA-seq reads 的 mapper），Crossbow（重测序数据分析云的软件工具），Myrna（对齐RNA-seq reads和分析差异基因表达的云计算软件工具）

5.1 下载官方索引

#下载小鼠参考基因组的索引和注释文件, 这里用常用的mm10
# 索引大小为3.2GB， 不建议自己下载基因组构建，可以直接下载索引文件，代码如下：
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已构建好存放在~/biosoft/referece

• 自建索引

#这里以构建人,hg38为例 
#构建参考基因组索引
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz 
tar -zxvf hg38.fa.gz  
cd ~/biosoft/reference/hg38/#在目标文件夹构建不用指定目录
bowtie2-build --threads 10 -f hg38.fa hg38

索引构建参考

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

• 批量处理，运行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这个函数
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
#依赖已经添加了，$bin_bowtie2可以直接换成bowtie2
#比对很慢
#-p/--threads NTHREADS
#-3/--trim3 <int>: Trim <int> bases from 3' (right) end of each read before alignment (default: 0).
#--local: 我的理解是不是精确比对，允许个别碱基错配
#-x: The basename of the index for the reference genome. 
#bowtie2输出格式默认为SAM，可以指定为多种其它格式，比如说BAM（SAM的二进制文件）等等
#利用samtools sort 进行排序并将sam格式转换为bam格式

• 也可以单独将SAM 文件转为 BAM 文件
  这个 BAM 文件就是后面 MACS2 的输入文件。
  samtools sort example.sam > example.bam

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

将.bam文件转换为.tdf文件，导入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
## 比对过程3分钟一个样本
bowtie2 -p 5 --local -x  $bowtie2_index -U  $id | samtools sort  -O bam  -@ 5 -o ~/project/align/${sample}.bam
done 

需要注意的是：genome_index 指的是用于bowtie2的索引文件（如下图），而不是参考基因组本身，构建过程参考前文5.1部分。genome_index 需要指定路径及其共用文件名，比如我的索引文件放在~/biosoft/referece/mm10目录下，但是需要输入的参数为~/biosoft/referece/mm10/mm10。最后一个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必须参数

bowtie2可选参数

5.3对bam文件进行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文件

参考这里
因为一个样品分成了多个lane进行测序，所以在进行peaks calling的时候，需要把bam进行合并,即一个样分成了多个SRR...数据,需要合并。

## 如果不用循环：
## 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重复

5.5.1 Sambamba去除PCR重复（推荐此方法）

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是线程数，-p显示过程，-r remove duplicates instead of just marking them
#处理后同时会给出索引文件

5.5.2 samtools markdup去除PCR重复（不推荐，还需要多一步才能进行）

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重复, -r是直接去除重复, 不加是直接标记
ls $files_bam  | xargs -i samtools index {} 

5.6计算比对率

#以比对后的bam文件进行计算
cd ~/project/align/
files_bam=~/project/align/*.bam
ls $files_bam  | while read id 
do 
samtools flagstat $id > $(basename -s .bam $id).stat
done
#获取比对率

6 使用macs2进行peaks calling

macs2包含一系列的子命令，其中最主要的就是callpeak， 官方提供了使用实例

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

• -t: 实验组的输出结果- -c: 对照组的输出结果
• -f: -t和-c提供文件的格式，可以是”ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT”, “ELANDMULTIPET” (for pair-end tags), “SAM”, “BAM”, “BOWTIE”, “BAMPE” “BEDPE” 任意一个。如果不提供这项，就是自动检测选择。
• -g: 基因组大小， 默认提供了hs, mm, ce, dm选项， 不在其中的话，比如说拟南芥，就需要自己提供了。
• -n: 输出文件的前缀名
• -B: 会保存更多的信息在bedGraph文件中，如fragment pileup, control lambda, -log10pvalue and -log10qvalue scores
• -q: q值，也就是最小的PDR阈值， 默认是0.05。q值是根据p值利用BH计算，也就是多重试验矫正后的结果。
• -p： 这个是p值，指定p值后MACS2就不会用q值了。
• -m: 和MFOLD有关，而MFOLD和MACS预构建模型有关，默认是5：50，MACS会先寻找100多个peak区构建模型，一般不用改。

cd ~/project/peaks/  
files_bams=~/project/rmdup/*.bam  #如果没有去除PCR重复就用files_bams=~/project/align/*.bam

#control.bam换成自己的对照组
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  
#双端测序结果用BAMPE??
#control.bam是对照组
#输出文件.bdg可以直接导入IGV可视化
#为何指定了输出文件夹，但还是输出在.bam所在的文件夹里？
#没有control文件也可以不用加，系统会用算法模拟一个，-t 后可接两个bam文件，因为正式实验一般测两个重复



#或者用下面这个循环
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

#没有control数据如下处理
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 对照组，可以接多个数据，空格分隔。
# -t 实验组，ChIP-seq数据。可以接多个数据，空格分隔。
# -n 输出文件名的前缀
# -g 有效基因组大小。软件有几个默认值，如hs指human，mm指mouse
# --outdir 输出结果的存储路径
# --bdg 输出 bedgraph 格式的文件
# -q FDR阈值, 默认为0.05
# NAME 是使用 -n 设置的输出文件前缀
# NAME_control_lambda.bdg 和 NAME_treat_pileup.bdg 
## bdg即为bedGraph格式，可以导入UCSC或者转换为bigwig格式。
# 两个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 设置的输出文件前缀

1. NAME_control_lambda.bdg 和 NAME_treat_pileup.bdg
   bdg即为bedGraph格式，可以导入UCSC或者转换为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会得到名为 NAME_model.pdf 的文件↓：

图片.png

首先我们要记住一点——所有数据分析的流程都基于实验原理。

做 ChIP-seq，就是为了确认蛋白-DNA结合位点。测序获得的 read 只是跟随着蛋白一起沉淀下来的 DNA 片段的末端，并不是真实的 DNA 和蛋白结合的位置。怎么找到蛋白和DNA真正的结合位置？
单端测序，对于测序的同一片段，reads 在正负链上的分布是平均的，并且正负链的 reads 之间会有一个间距。近似的形成“双峰”。因此，把来自于同一片段位于不同链上的 reads 相向移动，就能有效地降低信号中的噪音，找到结合的位置。而对 双端测序 而言，它本身测的就是文库的两端，因此不用建立模型和偏倚。我们只需要对 MACS 设置参数 --nodel 就能略过双峰模型建立的过程.

7 Peak 可视化

7.1 bam文件转化为bw(bigWig)或bedgraph文件

BAM转换为bigWig或bedGraph
BAM文件是SAM的二进制转换版，应该都知道。那么bigWig格式是什么？bigWig是wig或bedGraph的二进制版，存放区间的坐标轴信息和相关计分(score)，主要用于在基因组浏览器上查看数据的连续密度图，可用wigToBigWig从wiggle进行转换。
bedGraph和wig格式是什么? USCS的帮助文档称这两个格式数是已经过时的基因组浏览器图形轨展示格式，前者展示稀松型数据，后者展示连续性数据。目前推荐使用bigWig/bigBed这两种格式取代前两者。
为什么要用bigWig? 主要是因为BAM文件比较大，直接用于展示时对服务器要求较大。因此在GEO上仅会提供bw,即bigWig下载，便于下载和查看。如果真的感兴趣，则可以下载原始数据进行后续分析。

deeptools

使用deeptool进行可视化
deeptools提供bamCoverage和bamCompare进行格式转换，为了能够比较不同的样本，需要对先将基因组分成等宽分箱(bin)，统计每个分箱的read数，最后得到描述性统计值。对于两个样本，描述性统计值可以是两个样本的比率，或是比率的log2值，或者是差值。如果是单个样本，可以用SES方法进行标准化。
bamCoverage的基本用法

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

得到的bw文件就可以送去IGV/Jbrowse进行可视化。 这里的参数仅使用了-e/--extendReads和-bs/--binSize即拓展了原来的read长度，且设置分箱的大小。其他参数还有--filterRNAstrand {forward, reverse}: 仅统计指定正链或负链
--region/-r CHR:START:END: 选取某个区域统计--smoothLength: 通过使用分箱附近的read对分箱进行平滑化, 如果为了其他结果进行比较，还需要进行标准化，deeptools提供了如下参数：

--scaleFactor: 缩放系数
--normalizeUsingRPKMReads: Per Kilobase per Million mapped reads (RPKM)标准化
--normalizeTo1x: 按照1x测序深度(reads per genome coverage, RPGC)进行标准化
--ignoreForNormalization： 指定那些染色体不需要经过标准化
如果需要以100为分箱，并且标准化到1x，且仅统计某一条染色体区域的正链，输出格式为bedgraph,那么命令行可以这样写

• bam转换为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类似，只不过需要提供两个样本，并且采用SES方法进行标准化，于是多了--ratio参数。

• 批量把bam文件转为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 {} #转换前需要先构建index文件以.bai为后缀
#其实在Sambamba去除PCR重复时已经自动构建好了，存放在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 

转换的很慢
参数：--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进行可视化
  进入IGV目录，运行IGV

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

bash igv.sh即可打开桌面版IGV软件。
这里直接上传 bg/bw 文件：File->Load from File->选择sample.bg并打开。

7.3 bam文件转化为bed文件

#将bam文件转换为bed文件
bedtools bamtobed -i ~/project/align/SRR065175.bam  > ~/project/bed/SRR065175.bed

7.4 查看转录起始点TSS附近信号强度

deeptools软件里的computeMatrix命令可以查看TSS附近peaks强度，输入文件为 .bw，并给一个参考的注释文件（即 .bed文件下载方法见这里及这里我的 .bed文件已下载存储在~/biosoft/refseq_bed目录里。三个文件分别是hg19.bed, hg38.bed, mm10.bed），从而让软件查看样本在TSS（转录起始位点）附近是否有富集。

• 首先对单一样本绘图：

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
#输出文件在当前tss目录里 输出2个文件，-b和-a 是TSS上下游距离，-S是输入上文中bam转换成的.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 

• (批量处理)画10K附近，即TSS前后10K范围，可通过修改参数-a和-b来修改。

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命令作图。
#输出文件在当前tss目录里 输出2个文件，-b和-a 是TSS上下游距离，-S是输入上文中bam转换成的.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以热图的方式对覆盖进行可视化，用plotProfile以折线图的方式展示覆盖情况。
computeMatrix具有两个模式:scale-region和reference-point。前者用来信号在一个区域内分布，后者查看信号相对于某一个点的分布情况。无论是那个模式，都有有两个参数是必须的，-S是提供bigwig文件，-R是提供基因的注释信息。还有更多个性化的可视化选项。

7.5 使用R包对找到的peaks文件进行注释

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看看效果，检测软件找到的peaks是否真的合理，还可以配合rmarkdown来出自动化报告。
也可以使用其它代码进行下游分析； https://github.com/jmzeng1314/NGS-pipeline/tree/master/CHIPseq

matrix1_BMSC_d00_h3k27me3_R2_trimmed_TSS_2K_Heatmap_10K.png

7.6 peaks相关基因集的注释

都是得到感兴趣基因集，然后注释，分析方法等同于GEO数据挖掘课程或者转录组下游分析： https://github.com/jmzeng1314/GEO （有配套视频）

8 用homer软件寻找motif

conda install -c bioconda homer ,homer安装在 ~/anaconda3/envs/chipseq/share/homer下，然后使用其附带的配置脚本下载人和小鼠的基因组数据库（5G左右比较大,我的已经下载好存放在~/anaconda3/envs/chipseq/share/homer/data/genomes/）。

perl ~/anaconda3/envs/chipseq/share/homer/configureHomer.pl -install hg19#优先使用这行命令下载
#上面安装不了试试这行命令：perl ~/anaconda3/envs/chipseq/share/homer/configureHomer.pl  -install mm10 
## 我们上游分析是基于mm10找到的peaks文件

## 下载成功后会多出 ~/anaconda3/envs/chipseq/share/homer/data/genomes/mm10/ 文件夹, 共 4.9G
## 这个文件夹取决于你把homer这个软件安装到了什么地方。

## 或者用下面代码安装:这种安装会装在~/biosoft/homer文件夹下，运行会报错，把/genomes/mm10/ 和config文件复制粘贴到anaconda3文件夹下的homer即可。
cd ~/biosoft/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软件需要输入 .bed文件
但是使用起来很简单：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}  #自动创建目录
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，然后运行：nohup bash runMotif.sh 1>motif.log &
不仅仅找了motif，还顺便把peaks注释了一下。得到的后缀为peakAnn.xls 的文件就可以看到和使用R包注释的结果是差不多的。
结果分析参考这里
还可以使用meme来找motif，需要通过bed格式的peaks的坐标来获取fasta序列。MEME，链接：http://meme-suite.org/

其它高级分析

比如可以 比较不同的peaks文件，代码见：https://github.com/jmzeng1314/NGS-pipeline/blob/master/CHIPseq/step6-ChIPpeakAnno-Venn.R

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