寫在前面：參考-# 給學徒ChIP-seq數據處理流程
目錄
軟件安裝
[公共數據的獲取---sratoolkit --prefetch ]
[得到fastq格式的測序數據----sratoolkit fastq-dump]
fastq數據的質量控制 ----fastqc腻贰；trim_galore 多參數
比對以及質控 ----bowtie建立索引;samtool格式轉換
討論去除PCR重復與否
[批處理合并多條lane的測序數據 ---samtools的merge命令]
使用macs2找peaks
[IGV可視化]
[IGV-bigWig和bedgraph文件詳情]
deeptools可視化
使用R包注釋
[使用Homer找motif]
流程-軟件
1 需要的軟件
下載雷绢，使用conda
source activate chipseq
# 可以用search先進行檢索
conda install -y sra-tools
conda install -y trim-galore samtools
conda install -y deeptools homer meme
conda install -y macs2 bowtie bowtie2
2 數據獲取
使用srr號+prefetch
默認下載地址：/home/yangjy/ncbi/public/sra
數據分析地址：/home/yangjy/chipseq_test/rawdata
移動一下：使用mv
(chipseq) [yangjy@GSCG01 ~]$ cd /home/yangjy/ncbi/public/sra
(chipseq) [yangjy@GSCG01 sra]$ ls
SRR391032.sra SRR391033.sra SRR391034.sra SRR391041.sra SRR391043.sra
(chipseq) [yangjy@GSCG01 sra]$ mv *.sra /home/yangjy/chipseq_test/rawdata
(chipseq) [yangjy@GSCG01 sra]$ ls #無文件媚朦，已移動
(chipseq) [yangjy@GSCG01 sra]$ cd #回到家目錄
(chipseq) [yangjy@GSCG01 ~]$ cd /home/yangjy/chipseq_test/rawdata
(chipseq) [yangjy@GSCG01 rawdata]$ ls
SRR391032.sra SRR391033.sra SRR391034.sra SRR391041.sra SRR391043.sra
(chipseq) [yangjy@GSCG01 rawdata]$ ls -lh #查看細節(jié)
total 2.3G
-rw-rw-r--. 1 yangjy yangjy 474M Jan 31 11:21 SRR391032.sra
-rw-rw-r--. 1 yangjy yangjy 473M Jan 31 12:15 SRR391033.sra
-rw-rw-r--. 1 yangjy yangjy 406M Jan 31 12:17 SRR391034.sra
-rw-rw-r--. 1 yangjy yangjy 322M Jan 31 12:29 SRR391041.sra
-rw-rw-r--. 1 yangjy yangjy 597M Jan 31 12:32 SRR391043.sra
批量建立文件夾 mkdir{sra,raw,qc,clean,align,peacks,motif}
3 得到fastq格式的測序數據
制作配置文件;即把srr號與基因組的名稱對應起來：
RNAPII_S5P_1 SRR391032
RNAPII_S5P_2 SRR391033
RNAPII_S2P_1 SRR391034
H2Aub1_1 SRR391041
H3K36me3_1 SRR391043
用fastq-dump把sra格式轉成fastq格式(fq格式)
ls *sra |while read id; do ~/miniconda3/envs/chipseq/bin/fastq-dump $id;done
實際情況
(chipseq) [yangjy@GSCG01 rawdata]$ which fastq-dump
~/miniconda3/envs/chipseq/bin/fastq-dump
(chipseq) [yangjy@GSCG01 rawdata]$ ls *sra |while read id; do ~/miniconda3/envs/chipseq/bin/fastq-dump $id;done
Read 13532944 spots for SRR391032.sra
Written 13532944 spots for SRR391032.sra
Read 13662812 spots for SRR391033.sra
Written 13662812 spots for SRR391033.sra
Read 26571408 spots for SRR391034.sra
Written 26571408 spots for SRR391034.sra
Read 9780305 spots for SRR391041.sra
Written 9780305 spots for SRR391041.sra
Read 12313493 spots for SRR391043.sra
Written 12313493 spots for SRR391043.sra
(chipseq) [yangjy@GSCG01 rawdata]$ ls
SRR391032.fastq SRR391033.fastq SRR391034.fastq SRR391041.fastq SRR391043.fastq
SRR391032.sra SRR391033.sra SRR391034.sra SRR391041.sra SRR391043.sra
(chipseq) [yangjy@GSCG01 rawdata]$ ls -lh
total 18G
-rw-rw-r--. 1 yangjy yangjy 2.6G Jan 31 14:56 SRR391032.fastq
-rw-rw-r--. 1 yangjy yangjy 474M Jan 31 11:21 SRR391032.sra
-rw-rw-r--. 1 yangjy yangjy 2.6G Jan 31 14:57 SRR391033.fastq
-rw-rw-r--. 1 yangjy yangjy 473M Jan 31 12:15 SRR391033.sra
-rw-rw-r--. 1 yangjy yangjy 5.0G Jan 31 14:58 SRR391034.fastq
-rw-rw-r--. 1 yangjy yangjy 406M Jan 31 12:17 SRR391034.sra
-rw-rw-r--. 1 yangjy yangjy 1.9G Jan 31 14:58 SRR391041.fastq
-rw-rw-r--. 1 yangjy yangjy 322M Jan 31 12:29 SRR391041.sra
-rw-rw-r--. 1 yangjy yangjy 2.9G Jan 31 14:59 SRR391043.fastq
-rw-rw-r--. 1 yangjy yangjy 597M Jan 31 12:32 SRR391043.sra
注意：
這里的SRRxxx.sra格式轉換后為*.fastq格式当叭，用--gzip就能輸出*.fastq.gz格式, 能夠節(jié)省空間的同時也不會給后續(xù)比對軟件造成壓力, 比對軟件都支持吹埠。
可以使用
ls *sra |while read id; do ~/miniconda3/envs/chipseq/bin/fastq-dump --gzip --split-3 $id -O ./fastq/;done
結果：
(chipseq) [yangjy@GSCG01 rawdata]$ cd fastq/
(chipseq) [yangjy@GSCG01 fastq]$ ls -lh
total 3.7G
-rw-rw-r--. 1 yangjy yangjy 736M Jan 31 16:21 SRR391032.fastq.gz
-rw-rw-r--. 1 yangjy yangjy 741M Jan 31 16:25 SRR391033.fastq.gz
-rw-rw-r--. 1 yangjy yangjy 855M Jan 31 16:31 SRR391034.fastq.gz
-rw-rw-r--. 1 yangjy yangjy 508M Jan 31 16:34 SRR391041.fastq.gz
-rw-rw-r--. 1 yangjy yangjy 879M Jan 31 16:39 SRR391043.fastq.gz
/home/yangjy/chipseq_test/fastq 移動后*.fastq.gz所在位置
4 fastq數據的質量控制
使用trim_galore軟件進行質控
先使用fastqc看一下凡怎；[對照]
ls *.gz | while read id ; do /home/yangjy/miniconda3/envs/chipseq/bin/fastqc $id -O ../qc/;done
可以先看一下質量:
(chipseq) [yangjy@GSCG01 fastq]$ cd ../qc/
(chipseq) [yangjy@GSCG01 qc]$ ls -lh
total 2.8M
-rw-rw-r--. 1 yangjy yangjy 244K Jan 31 17:34 SRR391032_fastqc.html# 打開
-rw-rw-r--. 1 yangjy yangjy 330K Jan 31 17:34 SRR391032_fastqc.zip
-rw-rw-r--. 1 yangjy yangjy 243K Jan 31 17:34 SRR391033_fastqc.html# 打開
-rw-rw-r--. 1 yangjy yangjy 331K Jan 31 17:34 SRR391033_fastqc.zip
-rw-rw-r--. 1 yangjy yangjy 220K Jan 31 17:35 SRR391034_fastqc.html# 打開
-rw-rw-r--. 1 yangjy yangjy 272K Jan 31 17:35 SRR391034_fastqc.zip
-rw-rw-r--. 1 yangjy yangjy 248K Jan 31 17:36 SRR391041_fastqc.html# 打開
-rw-rw-r--. 1 yangjy yangjy 327K Jan 31 17:36 SRR391041_fastqc.zip
-rw-rw-r--. 1 yangjy yangjy 265K Jan 31 17:37 SRR391043_fastqc.html# 打開
-rw-rw-r--. 1 yangjy yangjy 366K Jan 31 17:37 SRR391043_fastqc.zip
fastqc報告
[协屡！黃色為警告实蓬；× 紅色為報錯茸俭；]
這個時候選擇trim_galore軟件進行過濾，單端測序數據的代碼如下安皱；
ls *.gz | while read id ; do
nohup
/home/yangjy/miniconda3/envs/chipseq/bin/trim_galore $id -q 20 --phred33 --stringency 3 --length 20 -e 0.1 --paired fq --gzip -o ../clean/ &
done
較慢调鬓，于是放在后臺處理。
結果：
(base) [yangjy@GSCG01 clean]$ ls -lh
total 3.6G
-rw-rw-r--. 1 yangjy yangjy 20 Jan 31 21:16 fq_trimmed.fq.gz
-rw-rw-r--. 1 yangjy yangjy 1.5K Jan 31 21:16 fq_trimming_report.txt
-rw-rw-r--. 1 yangjy yangjy 2.7K Jan 31 21:03 SRR391032.fastq.gz_trimming_report.txt
-rw-rw-r--. 1 yangjy yangjy 703M Jan 31 21:03 SRR391032_trimmed.fq.gz
-rw-rw-r--. 1 yangjy yangjy 2.7K Jan 31 21:06 SRR391033.fastq.gz_trimming_report.txt
-rw-rw-r--. 1 yangjy yangjy 709M Jan 31 21:06 SRR391033_trimmed.fq.gz
-rw-rw-r--. 1 yangjy yangjy 2.4K Jan 31 21:10 SRR391034.fastq.gz_trimming_report.txt
-rw-rw-r--. 1 yangjy yangjy 843M Jan 31 21:10 SRR391034_trimmed.fq.gz
-rw-rw-r--. 1 yangjy yangjy 2.7K Jan 31 21:12 SRR391041.fastq.gz_trimming_report.txt
-rw-rw-r--. 1 yangjy yangjy 491M Jan 31 21:12 SRR391041_trimmed.fq.gz
-rw-rw-r--. 1 yangjy yangjy 3.1K Jan 31 21:16 SRR391043.fastq.gz_trimming_report.txt
-rw-rw-r--. 1 yangjy yangjy 844M Jan 31 21:16 SRR391043_trimmed.fq.gz
再次進行QC酌伊，看看處理情況袖迎；尤其時質量差的時候
5 比對以及質控 ----bowtie建立索引;samtool格式轉換
參考基因組bowtie2_indexes官網
參考基因組及注釋下載
下載小鼠 [fasta格式的壓縮文件] （迅雷下載較快,使用winSCP上傳）
解壓（tar）合并（cat）建立索引（bowtie2-build）
得到---> Output files: "mm10.*.bt2" # 輸出的文件格式：mm10.*.bt2
inflating: mm10.1.bt2
inflating: mm10.2.bt2
inflating: mm10.3.bt2
inflating: mm10.4.bt2
inflating: mm10.rev.1.bt2
inflating: mm10.rev.2.bt2
inflating: make_mm10.sh
使用命令：
cd 到align目錄下（相對于目錄）
ls ../clean/*.gz |while read id;
do
file=$(basename $id )
sample=${file%%.*}
echo $file $sample
bowtie2 -p 5 -x /home/yangjy/project/epi/align/reference/mm10 -U $id | samtools sort -O bam -@ 5 -o - > ${sample}.bam
done
參數解讀：
-p 5 指5個線程；
-x 索引文件，不是參考基因組本身燕锥，是使用bowtie built 后的文件辜贵，參考;-x 后面接著gene-index，需要指定路徑（/home/yangjy/project/epi/align/reference/）及其共用文件名（mm10）
-U 后面接著fastq格式的文件（單末端測序參數）
注：命令把bowtie2 生成的sam文件通過管道|傳遞到samtools归形，將sam轉換為bam文件托慨，省去中間sam文件的空間占用
雙末端測序使用：
bowtie2 -p 5 -3 5 --local -x mm10 -1 example_1.fastq -2 example_2.fastq -S example.sam
然后再將SAM 文件轉為 BAM 文件：
samtools sort example.sam > example.bam
處理中：
H2Aub1_1_trimmed.fq H2Aub1_1_trimmed # 處理輸出名稱
9208673 reads; of these:
9208673 (100.00%) were unpaired; of these:
239111 (2.60%) aligned 0 times
6227498 (67.63%) aligned exactly 1 time
2742064 (29.78%) aligned >1 times
97.40% overall alignment rate
[bam_sort_core] merging from 0 files and 5 in-memory blocks....
結果：
(chipseq) [yangjy@GSCG01 align]$ ls -l
total 3342040
drwxrwxr-x. 2 yangjy yangjy 122 Feb 2 20:51 index
-rw-rw-r--. 1 yangjy yangjy 526164473 Jan 21 14:16 H2Aub1_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 873489325 Jan 21 14:26 H3K36me3_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 765203056 Jan 21 14:51 RNAPII_S2P_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 620491240 Jan 21 14:59 RNAPII_S5P_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 636886282 Jan 21 15:02 RNAPII_S5P_2_trimmed.bam
6 對bam文件進行QC----- samtools 【samtools flagstat 】
cd ~/project/epi/align
ls *.bam |xargs -i samtools index {} # 需要構建索引
ls *.bam | while read id ;do (nohup samtools flagstat $id > $(basename $id ".bam").stat & );done
nohup 命令&
得到：后綴為.bai的文件
-rw-rw-r--. 1 yangjy yangjy 502M Jan 21 14:16 H2Aub1_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 3.0M Feb 2 21:44 H2Aub1_1_trimmed.bam.bai
-rw-rw-r--. 1 yangjy yangjy 834M Jan 21 14:26 H3K36me3_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 2.8M Feb 2 21:44 H3K36me3_1_trimmed.bam.bai
drwxrwxr-x. 2 yangjy yangjy 122 Feb 2 20:51 index
-rw-rw-r--. 1 yangjy yangjy 730M Jan 21 14:51 RNAPII_S2P_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 2.6M Feb 2 21:44 RNAPII_S2P_1_trimmed.bam.bai
-rw-rw-r--. 1 yangjy yangjy 592M Jan 21 14:59 RNAPII_S5P_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 3.7M Feb 2 21:44 RNAPII_S5P_1_trimmed.bam.bai
-rw-rw-r--. 1 yangjy yangjy 608M Jan 21 15:02 RNAPII_S5P_2_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 3.7M Feb 2 21:44 RNAPII_S5P_2_trimmed.bam.bai
以及：后綴為.stat的文件
total 3358004
-rw-rw-r--. 1 yangjy yangjy 526164473 Jan 21 14:16 H2Aub1_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 3109504 Feb 2 21:44 H2Aub1_1_trimmed.bam.bai
-rw-rw-r--. 1 yangjy yangjy 386 Feb 2 21:48 H2Aub1_1_trimmed.stat
-rw-rw-r--. 1 yangjy yangjy 873489325 Jan 21 14:26 H3K36me3_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 2912232 Feb 2 21:44 H3K36me3_1_trimmed.bam.bai
-rw-rw-r--. 1 yangjy yangjy 388 Feb 2 21:48 H3K36me3_1_trimmed.stat
drwxrwxr-x. 2 yangjy yangjy 122 Feb 2 20:51 index
-rw-rw-r--. 1 yangjy yangjy 765203056 Jan 21 14:51 RNAPII_S2P_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 2697544 Feb 2 21:44 RNAPII_S2P_1_trimmed.bam.bai
-rw-rw-r--. 1 yangjy yangjy 388 Feb 2 21:48 RNAPII_S2P_1_trimmed.stat
-rw-rw-r--. 1 yangjy yangjy 620491240 Jan 21 14:59 RNAPII_S5P_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 3805200 Feb 2 21:44 RNAPII_S5P_1_trimmed.bam.bai
-rw-rw-r--. 1 yangjy yangjy 388 Feb 2 21:48 RNAPII_S5P_1_trimmed.stat
-rw-rw-r--. 1 yangjy yangjy 636886282 Jan 21 15:02 RNAPII_S5P_2_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 3789176 Feb 2 21:44 RNAPII_S5P_2_trimmed.bam.bai
-rw-rw-r--. 1 yangjy yangjy 388 Feb 2 21:48 RNAPII_S5P_2_trimmed.stat
7 討論去除PCR重復與否
參考文章：
仔細探究picard的MarkDuplicates 是如何行使去除PCR重復reads功能的
仔細探究samtools的rmdup是如何行使去除PCR重復reads功能的
在bam文件夾下：
ls *.bam | while read id ;do (nohup samtools markdup -r $id $(basename $id ".bam").mergeBam & );done
得到：
(chipseq) [yangjy@GSCG01 align]$ ls -l *.rmdup.bam
-rw-rw-r--. 1 yangjy yangjy 474047646 Feb 2 21:59 H2Aub1_1_trimmed.rmdup.bam
-rw-rw-r--. 1 yangjy yangjy 830716373 Feb 2 22:00 H3K36me3_1_trimmed.rmdup.bam
-rw-rw-r--. 1 yangjy yangjy 632581765 Feb 2 22:00 RNAPII_S2P_1_trimmed.rmdup.bam
-rw-rw-r--. 1 yangjy yangjy 368497642 Feb 2 21:59 RNAPII_S5P_1_trimmed.rmdup.bam
-rw-rw-r--. 1 yangjy yangjy 371445769 Feb 2 21:59 RNAPII_S5P_2_trimmed.rmdup.bam
通過查看*.stat 文件：
(chipseq) [yangjy@GSCG01 align]$ cat H2Aub1_1_trimmed.stat #使用cat查看
9208673 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary 0條比對到第二個地方；
0 + 0 supplementary
0 + 0 duplicates
8969562 + 0 mapped (97.40% : N/A) #比對率97.40%
0 + 0 paired in sequencing
0 + 0 read1 因為是單端
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
來自于網頁（https://www.biostars.org/p/333465/）
45084184 + 0 in total (QC-passed reads + QC-failed reads)
4717987 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
45084184 + 0 mapped (100.00%:-nan%)
40366197 + 0 paired in sequencing
20273012 + 0 read1
20093185 + 0 read2
39146254 + 0 properly paired (96.98%:-nan%)
39644363 + 0 with itself and mate mapped
721834 + 0 singletons (1.79%:-nan%)
260722 + 0 with mate mapped to a different chr
235361 + 0 with mate mapped to a different chr (mapQ>=5)
解讀：
QC-passed reads的數量為9208673暇榴，未通過的0厚棵；意味著一共有9208673條reads
下面都是0？蔼紧？ 可能是單端測序的原因（雙末端測序下載應該有一側有數值）
注：可以使用 multiqc 對 *.stat 進行可視化
可以再對*.rmdup.bam 進行統(tǒng)計：【samtools flagstat 】
ls *.rmdup.bam |xargs -i samtools index {} 只有建立索引婆硬，才能進行--統(tǒng)計 ；得到后綴為.bai 的文件
ls *.rmdup.bam | while read id ;do (nohup samtools flagstat $id > $(basename $id ".bam").stat & );done
綜上兩個*.bam 文件：
(chipseq) [yangjy@GSCG01 align]$ ls -lh *.bam
-rw-rw-r--. 1 yangjy yangjy 502M Jan 21 14:16 H2Aub1_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 453M Feb 2 21:59 H2Aub1_1_trimmed.rmdup.bam #去除PCR重復的
-rw-rw-r--. 1 yangjy yangjy 834M Jan 21 14:26 H3K36me3_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 793M Feb 2 22:00 H3K36me3_1_trimmed.rmdup.bam
-rw-rw-r--. 1 yangjy yangjy 730M Jan 21 14:51 RNAPII_S2P_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 604M Feb 2 22:00 RNAPII_S2P_1_trimmed.rmdup.bam
-rw-rw-r--. 1 yangjy yangjy 592M Jan 21 14:59 RNAPII_S5P_1_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 352M Feb 2 21:59 RNAPII_S5P_1_trimmed.rmdup.bam
-rw-rw-r--. 1 yangjy yangjy 608M Jan 21 15:02 RNAPII_S5P_2_trimmed.bam
-rw-rw-r--. 1 yangjy yangjy 355M Feb 2 21:59 RNAPII_S5P_2_trimmed.rmdup.bam
由此奸例，可以對這些*.bam文件進行peak calling
此前彬犯，對【一個樣品分成了多個lane進行測序】情況可以把bam進行合并（merge）。
使用：samtools merge 參數查吊；
參看比對率： grep 'N/A' *.stat |grep %
(chipseq) [yangjy@GSCG01 align]$ grep 'N/A' *.stat |grep %
H2Aub1_1_trimmed.rmdup.stat:7744878 + 0 mapped (97.01% : N/A)
H2Aub1_1_trimmed.stat:8969562 + 0 mapped (97.40% : N/A)
H3K36me3_1_trimmed.rmdup.stat:11007740 + 0 mapped (98.82% : N/A)
H3K36me3_1_trimmed.stat:11737364 + 0 mapped (98.89% : N/A)
RNAPII_S2P_1_trimmed.rmdup.stat:18461592 + 0 mapped (96.32% : N/A)
RNAPII_S2P_1_trimmed.stat:25018815 + 0 mapped (97.26% : N/A)
RNAPII_S5P_1_trimmed.rmdup.stat:6042684 + 0 mapped (95.30% : N/A)
RNAPII_S5P_1_trimmed.stat:11801448 + 0 mapped (97.53% : N/A)
RNAPII_S5P_2_trimmed.rmdup.stat:6123727 + 0 mapped (96.43% : N/A)
RNAPII_S5P_2_trimmed.stat:12182337 + 0 mapped (98.18% : N/A)
注意格式：常用的mapping軟件bwa & sam格式簡介
samtools view
(base) [yangjy@GSCG01 dum]$ samtools view H2Aub1.merge.rmdup.bam | cut -f 1 | sort | uniq -c | wc -l
18090450
8 MACS 找call peak
有了*.bam文件后就可以進行call了
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。
處理：
cd bam 文件目錄
使用命令：
macs2 callpeak -c Control_1_trimmed.bam -t H3K36me3_1_trimmed.bam -f BAM -B -g mm -n H3K36me3
# macs2 callpeak -c <Control.bam> -t <test.bam> -f BAM -B -g mm -n testname
注：chipseq環(huán)境下的macs一直報錯；
退出conda環(huán)境渣锦，使用命令---> 成功
(base) [yangjy@GSCG01 align1]$ conda list macs2
packages in environment at /home/yangjy/miniconda3:
#########
Name Version Build Channel
macs2 2.2.7.1 py38h0213d0e_1 https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
處理時的信息參考
得到文件：
-rw-rw-r--. 1 yangjy yangjy 565332244 Feb 3 20:40 H3K36me3_control_lambda.bdg
-rw-rw-r--. 1 yangjy yangjy 99001 Feb 3 20:39 H3K36me3_model.r
-rw-rw-r--. 1 yangjy yangjy 264716 Feb 3 20:40 H3K36me3_peaks.narrowPeak
-rw-rw-r--. 1 yangjy yangjy 302627 Feb 3 20:40 H3K36me3_peaks.xls
-rw-rw-r--. 1 yangjy yangjy 179908 Feb 3 20:40 H3K36me3_summits.bed
-rw-rw-r--. 1 yangjy yangjy 689424712 Feb 3 20:40 H3K36me3_treat_pileup.bdg
文件解讀：
1 Excel文件
前幾行為運行命令以及信息
后面從左到右分別是染色體硝岗，peak的起始終止位置，長度袋毙，峰值位置型檀，高度，pvalue值記憶富集程度
(base) [yangjy@GSCG01 peaks]$ cat H2Aub1_peaks.xls
This file is generated by MACS version 2.1.0.20150731
Command line: callpeak -c Control.merge.bam -t H2Aub1.merge.bam -f BAM -B -g mm -n H2Aub1 --outdir ../peaks
ARGUMENTS LIST:
name = H2Aub1
format = BAM
ChIP-seq file = ['H2Aub1.merge.bam']
control file = ['Control.merge.bam']
effective genome size = 1.87e+09
band width = 300
model fold = [5, 50]
qvalue cutoff = 5.00e-02
Larger dataset will be scaled towards smaller dataset.
Range for calculating regional lambda is: 1000 bps and 10000 bps
Broad region calling is off
tag size is determined as 49 bps
total tags in treatment: 22199430
tags after filtering in treatment: 17515747
maximum duplicate tags at the same position in treatment = 1
Redundant rate in treatment: 0.21
total tags in control: 14660284
tags after filtering in control: 12331048
maximum duplicate tags at the same position in control = 1
Redundant rate in control: 0.16
d = 83
alternative fragment length(s) may be 83 bps
chr start end length abs_summit pileup -log10(pvalue) fold_enrichment -log10(qvalue) name
chr1 4492658 4492755 98 4492668 11.97 9.08008 6.48398 4.02077 H2Aub1_peak_1
chr1 4492851 4493214 364 4493158 14.78 16.26556 10.20087 8.94022 H2Aub1_peak_2
chr1 4493287 4493599 313 4493468 9.15 9.39303 6.56103 4.02077 H2Aub1_peak_3
chr1 4496519 4496628 110 4496554 8.45 7.41898 5.67787 2.72721 H2Aub1_peak_4
chr1 4497051 4497263 213 4497103 9.86 9.39303 7.01601 4.02077 H2Aub1_peak_5
2 BED文件
染色體听盖，peak的起始終止位置胀溺，
(base) [yangjy@GSCG01 mergeBam]$ less RNAPII_S5P_summits.bed
chr1 3150842 3150843 RNAPII_S5P_peak_1 8.17172
chr1 3953635 3953636 RNAPII_S5P_peak_2 4.27895
chr1 4085551 4085552 RNAPII_S5P_peak_3 9.22983
............................
chrY 90834489 90834490 RNAPII_S5P_peak_56933 7.80872
chrY 90835530 90835531 RNAPII_S5P_peak_56934 6.31523
chrY 90837388 90837389 RNAPII_S5P_peak_56935 6.23419
chrY 90838538 90838539 RNAPII_S5P_peak_56936 12.47469
chrY 90839706 90839707 RNAPII_S5P_peak_56937 5.36220
chrY 90840114 90840115 RNAPII_S5P_peak_56938 6.03569
可以使用wc -l *.bed 查看一下每個基因有多少個peaks
0 Control_summits.bed
1089 H2Aub1_summits.bed
40846 H3K36me3_summits.bed
26027 Ring1B_summits.bed
41816 RNAPII_8WG16_summits.bed
19973 RNAPII_S2P_summits.bed
56942 RNAPII_S5P_summits.bed
72262 RNAPII_S7P_summits.bed
MockIP是control，所以它自己跟自己比較皆看，是沒有peaks的;
NAMEpeaks.xls: 以表格形式存放peak信息仓坞，雖然后綴是xls，但其實能用文本編輯器打開腰吟，和bed格式類似无埃，但是以1為基，而bed文件是以0為基.也就是說xls的坐標都要減一才是bed文件的坐標
NAMEpeaks.narrowPeak NAMEpeaks.broadPeak 類似蝎困。后面4列表示為录语， integer score for display， fold-change禾乘，-log10pvalue澎埠，-log10qvalue，relative summit position to peak start始藕。內容和NAMEpeaks.xls基本一致蒲稳，適合用于導入R進行分析
NAMEsummits.bed：記錄每個peak的peak summits，也就是記錄極值點的位置伍派。MACS建議用該文件尋找結合位點的motif江耀。
NAME_model.r，能通過NAME_model.r作圖诉植，得到是基于你提供數據的peak模型
參考：http://www.reibang.com/p/2b8e2ea26665
9 將BAM文件轉化為bw文件祥国，使用deeptools
deeptools提供
bamCoverage和bamCompare進行格式轉換，為了能夠比較不同的樣本晾腔，需要對先將基因組分成等寬分箱(bin)舌稀，統(tǒng)計每個分箱的read數，最后得到描述性統(tǒng)計值灼擂。對于兩個樣本壁查，描述性統(tǒng)計值可以是兩個樣本的比率，或是比率的log2值剔应，或者是差值睡腿。如果是單個樣本语御，可以用SES方法進行標準化
bamCoverage的基本用法
`bamCoverage -b test.bam -o test.bw
activate chipseq
bamCoverage -e 170 -bs 10 -b ap2_chip_rep1_2_sorted.bam -o ap2_chip_rep1_2.bw
# ap2_chip_rep1_2_sorted.bam是前期比對得到的BAM文件
得到的bw文件就可以送去IGV/Jbrowse進行可視化。
參數使用 -e/--extendReads和-bs/--binSize 即拓展原來的read長度席怪，且設置分箱的大小
如果需要以100為分箱应闯，并且標準化到1x，且僅統(tǒng)計某一條染色體區(qū)域的正鏈何恶，輸出格式為bedgraph
bamCoverage -e 170 -bs 100 -of bedgraph -r Chr4:12985884:12997458 --normalizeTo1x 100000000 -b 02-read-alignment/ap2_chip_rep1_1_sorted.bam -o chip.bedgraph
bamCompare和bamCoverage類似
需要提供兩個樣本孽锥，并且采用SES方法進行標準化，多了--ratio參數细层。
首先把bam文件轉為bw文件惜辑，詳情：wig、bigWig和bedgraph文件詳解
cd bam文件目錄下
activate chipseq
ls *.bam |xargs -i samtools index {} # 如果沒有索引疫赎，需要先建立索引
ls *.bam |while read id;do
nohup bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.bw &
done
參數 --normalizeUsing CPM
以上的bed盛撑，bdg 以及bw 文件可以載入IGV進行可視化；
head *.bed 查看
通過：awk '{print $1":"$2"-"$3}' RNAPII_S5P_summits.bed |head可以在屏幕顯示出 chr 和坐標信息捧搞；輸入IGV來檢查這個位置的peaks
bw文件反應bam的測序深度
10 查看TSS附件信號強度：使用deeptools
需要下載參考的注釋文件tss的bed格式
下載方法
首先對單一樣本繪圖：
使用computeMatrix函數對tss上下10kb(10000)計算信號強度
mkdir -p ~/project/epi/tss
cd ~/project/epi/tss
computeMatrix reference-point --referencePoint TSS -p 15 \
-b 10000 -a 10000 \
-R /home/yangjy/project/epi/peaks/mm10.tss.bed \
-S /home/yangjy/project/epi/align/mergeBam/H2Aub1.bw \
--skipZeros -o matrix1_test_TSS.gz \
--outFileSortedRegions regions1_test_genes.bed 輸出文件
進行批量處理：
##
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
bed=/home/yangjy/project/epi/peaks/mm10.tss.bed
for id in /home/yangjy/project/epi/align/mergeBam/*.bw ;
do
echo $id
file=$(basename $id )
sample=${file%%.*}
echo $sample
computeMatrix reference-point --referencePoint TSS -p 15 \
-b 2000 -a 2000 \
-R $bed \
-S $id \
--skipZeros -o matrix1_${sample}_TSS_2K.gz \
--outFileSortedRegions regions1_${sample}_TSS_2K.bed
plotHeatmap -m matrix1_${sample}_TSS_2K.gz -out ${sample}_Heatmap_2K.png
plotHeatmap -m matrix1_${sample}_TSS_2K.gz -out ${sample}_Heatmap_2K.pdf --plotFileFormat pdf --dpi 720
plotProfile -m matrix1_${sample}_TSS_2K.gz -out ${sample}_Profile_2K.png
plotProfile -m matrix1_${sample}_TSS_2K.gz -out ${sample}_Profile_2K.pdf --plotFileFormat pdf --perGroup --dpi 720
done
輸出：
-rw-rw-r--. 1 yangjy yangjy 10108967 Jan 27 09:46 Ring1B_Heatmap_10K.pdf
-rw-rw-r--. 1 yangjy yangjy 1093343 Jan 27 09:44 Ring1B_Heatmap_10K.png
-rw-rw-r--. 1 yangjy yangjy 4434818 Jan 27 00:33 Ring1B_Heatmap_2K.pdf
-rw-rw-r--. 1 yangjy yangjy 639686 Jan 27 00:32 Ring1B_Heatmap_2K.png
-rw-rw-r--. 1 yangjy yangjy 392676 Jan 27 09:47 Ring1B_Profile_10K.pdf
-rw-rw-r--. 1 yangjy yangjy 38805 Jan 27 09:46 Ring1B_Profile_10K.png
-rw-rw-r--. 1 yangjy yangjy 385755 Jan 27 00:33 Ring1B_Profile_2K.pdf
-rw-rw-r--. 1 yangjy yangjy 35869 Jan 27 00:33 Ring1B_Profile_2K.png
圖片：
Control_Profile_2K.png
為了統(tǒng)計全基因組范圍的peak在基因特征的分布情況抵卫，也就是需要用到computeMatrix計算，用plotHeatmap以熱圖的方式對覆蓋進行可視化胎撇，用plotProfile以折線圖的方式展示覆蓋情況介粘。
computeMatrix具有兩個模式:scale-region和reference-point。前者用來信號在一個區(qū)域內分布晚树，后者查看信號相對于某一個點的分布情況姻采。無論是那個模式，都有有兩個參數是必須的爵憎，-S是提供bigwig文件慨亲，-R是提供基因的注釋信息。還有更多個性化的可視化選項宝鼓。
11 使用R包對找到的peaks文件進行注釋(還不熟悉）
12 peaks相關基因集的注釋
homer軟件來 尋找motif
homer軟件處理的原理：參考
先準備好格式化文件刑棵，再使用命令；
1 提供包含基因組坐標的文件 比如 peak文件或BED文件
2 分析peak文件中富集的motif：
findMotifsGenome.pl <peak/BED file> <genome> <output directory> -size
簡單例子：
findMotifsGenome.pl ERpeaks.txt hg18 ER_MotifOutput/ -size 200 -mask
-mask 使用repeated-mask序列
-size 設置motif長度
-len : motif的長度設置愚铡，默認8蛉签，10，12沥寥，越大越消耗計算資源碍舍。處理后得到的文件
homerMotifs.motifs <#> : these are the output files from the de novo motif finding, separated by motif length, and represent separate runs of the algorithm.
homerMotifs.all.motifs : Simply the concatenated file composed of all the homerMotifs.motifs<#> files.
motifFindingParameters.txt : 記錄執(zhí)行findMotifsGenome.pl的命令
knownResults.txt : 記錄motif的統(tǒng)計數據，text file(open in EXCEL).
seq.autonorm.tsv : autonormalization statistics for lower-order oligo normalization.
homerResults.html : de novo motif finding的格式化輸出
- 文件解讀：
1 txt文件： 記錄執(zhí)行findMotifsGenome.pl的命令以及motif的統(tǒng)計數據营曼，
(chipseq) [yangjy@GSCG01 motif]$ cat H2Aub1.annLog.txt
Peak file = homer_peaks.tmp
Genome = mm10
Organism = mouse
Peak/BED file conversion summary:
BED/Header formatted lines: 0
peakfile formatted lines: 1089
Duplicated Peak IDs: 0
Peak File Statistics:
Total Peaks: 1089
Redundant Peak IDs: 0
Peaks lacking information: 0 (need at least 5 columns per peak)
Peaks with misformatted coordinates: 0 (should be integer)
Peaks with misformatted strand: 0 (should be either +/- or 0/1)
Peak file looks good!
Reading Positions...
-----------------------
Finding Closest TSS...
Annotating:....................
Annotation Number of peaks Total size (bp) Log2 Ratio (obs/exp) LogP enrichment (+values depleted)
3UTR 38.0 20746162 2.146 -29.987
miRNA 0.0 31126 -0.018 0.013
ncRNA 22.0 3439613 3.950 -42.243
TTS 20.0 27176332 0.830 -4.492
pseudo 0.0 540203 -0.291 0.224
Exon 205.0 34540955 3.842 -376.610
Intron 400.0 947608369 0.028 -1.128
Intergenic 163.0 1564039797 -1.990 464.277
Promoter 183.0 30063639 3.878 -339.028
5UTR 58.0 2354993 5.895 -184.382
snoRNA 0.0 19 -0.000 0.000
rRNA 0.0 5631 -0.003 0.002
NOTE: If this part takes more than 2 minutes, there is a good chance
your machine ran out of memory: consider hitting ctrl+C and rerunning
the command with "-noann"
To capture annotation stats in a file, use "-annStats <filename>" next time
Annotating:....................
Annotation Number of peaks Total size (bp) Log2 Ratio (obs/exp) LogP enrichment (+values depleted)
3UTR 38.0 20746162 2.146 -30.001
Other 0.0 7154386 -1.989 2.964
RC? 0.0 10979 -0.007 0.005
RNA 0.0 113962 -0.066 0.047
miRNA 0.0 31126 -0.018 0.013
ncRNA 22.0 3439613 3.950 -42.252
TTS 20.0 27176332 0.831 -4.497
LINE 0.0 521060815 -8.076 240.272
LINE? 0.0 8168 -0.005 0.003
srpRNA 0.0 43307 -0.026 0.018
SINE 0.0 193760064 -6.452 83.282
RC 0.0 65629 -0.039 0.027
tRNA 0.0 266147 -0.151 0.110
DNA? 0.0 142070 -0.082 0.059
pseudo 0.0 540203 -0.291 0.224
DNA 1.0 28358286 -3.553 9.244
Exon 205.0 34540955 3.842 -376.698
Intron 243.0 601725627 -0.035 1.055
Intergenic 103.0 784206321 -1.656 135.581
Promoter 183.0 30063639 3.879 -339.106
5UTR 58.0 2354993 5.895 -184.408
snoRNA 0.0 19 -0.000 0.000
LTR? 0.0 192563 -0.111 0.080
scRNA 0.0 577291 -0.309 0.239
CpG-Island 206.0 3107258 7.324 -865.124
Low_complexity 1.0 18489411 -2.936 5.514
LTR 3.0 292929648 -5.336 115.531
Simple_repeat 6.0 55283218 -1.931 10.524
snRNA 0.0 236859 -0.135 0.098
Unknown 0.0 1234764 -0.596 0.511
SINE? 0.0 29758 -0.018 0.012
Satellite 0.0 3685416 -1.338 1.526
rRNA 0.0 165655 -0.096 0.069
Counting Tags in Peaks from each directory...
Organism: mouse
Loading Gene Informaiton...
readline() on closed filehandle IN at /home/yangjy/miniconda3/envs/chipseq/bin/annotatePeaks.pl line 3741, <IN> line 2178.
readline() on closed filehandle IN at /home/yangjy/miniconda3/envs/chipseq/bin/annotatePeaks.pl line 3751.
Outputing Annotation File...
Done annotating peaks file
2 xls文件：記錄染色體，起始位點愚隧，鏈+/-蒂阱， peaks落在什么位置（UTR锻全， promoter-TSS，exon 录煤， intron）鳄厌，距TSS最近PromoterID的注釋距離...等。
(chipseq) [yangjy@GSCG01 motif]$ head H2Aub1.peakAnn.xls
PeakID (cmd=annotatePeaks.pl homer_peaks.tmp mm10) Chr Start End Strand Peak Score Focus Ratio/Region Size Annotation Detailed Annotation Distance to TSS Nearest PromoterID Entrez ID Nearest Unigene Nearest Refseq Nearest Ensembl Gene Name Gene Alias Gene Description Gene Type
H2Aub1.log_peak_2 chr1 4493157 4493158 + 0 NA exon (NM_001289464, exon 3 of 4) exon (NM_001289464, exon 3 of 4) 4197NM_001289464
H2Aub1.log_peak_62 chr10 19849898 19849899 + 0 NA intron (NM_029529, intron 1 of 1) intron (NM_029529, intron 1 of 1) 1561 NM_029529
H2Aub1.log_peak_141 chr11 96354139 96354140 + 0 NA promoter-TSS (NR_155303) promoter-TSS (NR_155303) -745NR_155303
H2Aub1.log_peak_1020 chr9 32542410 32542411 + 0 NA promoter-TSS (NM_008026) promoter-TSS (NM_008026) -958NM_008026
H2Aub1.log_peak_269 chr13 112651870 112651871 + 0 NA intron (NM_010029, intron 1 of 20) intron (NM_010029, intron 1 of 20) 440 NM_001145885
H2Aub1.log_peak_129 chr11 93099621 93099622 + 0 NA 5' UTR (NM_028296, exon 1 of 9) 5' UTR (NM_028296, exon 1 of 9) 331 NM_028296
H2Aub1.log_peak_442 chr18 64346742 64346743 + 0 NA intron (NM_194268, intron 1 of 1) intron (NM_194268, intron 1 of 1) 6378 NM_194268
3 homer*.motifs文件：
每個motif信息是一塊妈踊，均以>開頭了嚎。>所在行的信息以tab分隔。 motif首行信息解釋：
>ACCTGCTTGAAA 19-ACCTGCTTGAAA 10.757736 -456.639967 0 T:120.0(0.46%),B:1.8(0.01%),P:1e-198
0.997 0.001 0.001 0.001
0.001 0.995 0.001 0.003
0.001 0.997 0.001 0.001
0.001 0.001 0.001 0.997
0.003 0.001 0.995 0.001
0.001 0.997 0.001 0.001
0.001 0.001 0.001 0.997
0.003 0.001 0.003 0.993
0.001 0.001 0.997 0.001
一致性序列：如圖上的>ACCTGCTTGAAA
Motif名稱：如圖上的19-ACCTGCTTGAAA
比值檢測概率的log值：10.757736
P-value得log值： -456.639967
占位符：如上圖得0廊营，不具有任何信息
逗號分隔得富集信息歪泳，如：T:120.0(0.46%),B:1.8(0.01%),P:1e-198
其中富集信息的解讀：
T:120.0(0.46%),B:1.8(0.01%),P:1e-198
- T表示帶有該motif的目標序列在總的目標序列（target）中得百分比
- B表示帶有該motif的背景序列在總的背景序列（background）中得百分比
- P表示最終富集的p-value
逗號分隔的motif統(tǒng)計信息，如:
Tpos:68.2,Tstd:51.5,Bpos:123.7,Bstd:43.3,StrandBias:3.0,Multiplicity:1.00
- Tpos：motif在目標序列中的平均位置
- Tstd：motif在目標序列中位置的標準差
- Bpos：motif在背景序列中的平均位置
- Bstd：motif在背景序列中位置的標準差
- StrandBias: 鏈偏好性露筒，在正義鏈上的motif數與反義鏈motif數的比值的log
- Multiplicity：具有一個或多個結合位點的序列中每個序列的平均出現次數;(大于等于1）
其他可視化軟件：
1 ngs.plot
bam文件建立索引
根據功能原件的基因組坐標索引bam文件
計算功能原件區(qū)域富集信號豐度
參考畫圖：富集輪廓圖和熱圖
2 將peakAnn.xls文件導入到本地呐伞，通過excel透視功能進行可視化
ATAC參考：https://www.it610.com/article/1224820226680524800.htm
callpeaks文件 解讀: https://blog.csdn.net/sunyu_03/article/details/82633799
熱圖軟件 deeptools：https://blog.csdn.net/xiaomotong123/article/details/108734727
