[实战1] Polycomb associates genome
2021-02-01 本文已影响0人
shine_9457
写在前面：参考-# 给学徒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 #查看细节
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格式, 能够节省空间的同时也不会给后续比对软件造成压力, 比对软件都支持。
可以使用
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 进行统计：【samtools flagstat 】
ls *.rmdup.bam |xargs -i samtools index {} 只有建立索引，才能进行--统计 ；得到后缀为.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环境下的macs一直报错；
退出conda环境，使用命令---> 成功
(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模型
参考：https://www.jianshu.com/p/2b8e2ea26665
9 将BAM文件转化为bw文件，使用deeptools
deeptools提供
bamCoverage和bamCompare进行格式转换，为了能够比较不同的样本，需要对先将基因组分成等宽分箱(bin)，统计每个分箱的read数，最后得到描述性统计值。对于两个样本，描述性统计值可以是两个样本的比率，或是比率的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，且仅统计某一条染色体区域的正链，输出格式为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
为了统计全基因组范围的peak在基因特征的分布情况，也就是需要用到computeMatrix计算，用plotHeatmap以热图的方式对覆盖进行可视化，用plotProfile以折线图的方式展示覆盖情况。
computeMatrix具有两个模式:scale-region和reference-point。前者用来信号在一个区域内分布，后者查看信号相对于某一个点的分布情况。无论是那个模式，都有有两个参数是必须的，-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 : 记录执行findMotifsGenome.pl的命令
knownResults.txt : 记录motif的统计数据，text file(open in EXCEL).
seq.autonorm.tsv : autonormalization statistics for lower-order oligo normalization.
homerResults.html : de novo motif finding的格式化输出
- 文件解读：
1 txt文件： 记录执行findMotifsGenome.pl的命令以及motif的统计数据，
(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统计信息，如:
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文件
计算功能原件区域富集信号丰度
参考画图：富集轮廓图和热图
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
