Chipseq分析步骤（1)

chipseq步骤简记
第一步 环境搭建

#####01 环境搭建######
#linux前期准备，R语言去看R脚本
wget -c https://mirrors.bfsu.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh
## 软链接即可
cd  ~
ln  -s  /home/t_linux/Miniconda3-latest-Linux-x86_64.sh  ./
## 安装，安装过程只需要输入 yes 或者按 Enter
bash  Miniconda3-latest-Linux-x86_64.sh 
## 重新激活环境
source  ~/.bashrc 
## 查看 conda 的帮助文档conda --help
## 配置镜像
# 下面四行配置北京大学的conda的channel地址（首选）
conda config --add channels https://mirrors.pku.edu.cn/anaconda/pkgs/main/
conda config --add channels https://mirrors.pku.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.pku.edu.cn/anaconda/cloud/bioconda/
conda config --set show_channel_urls yes
# 删除defaults频道
sed -i '/defaults/d' ~/.condarc
## 配置镜像成功
# 查看配置结果
cat ~/.condarc

conda env list
###创建一个chipseq的环境
conda create -n chipseq 
#激活chipseq环境
conda activate chipseq  
####安装软件####如果安装不上就有两种可能？1.镜像源问题，删除原有镜像源，添加新的。2.conda版本问题，使用命令更新conda:conda update -n base conda;conda update --all
conda install -y sra-tools
conda install fastqc
conda install multiqc
conda install trim-galore
conda install -c bioconda bowtie2
conda install -c bioconda samtools
conda install -c bioconda picard
python --version
conda install -y python=3.8
conda install -c bioconda macs2
conda install -c bioconda deeptools

conda activate chipseq

第一步 准备参考基因组文件

#构建参考基因组索引Dec. 2013 (GRCh38/hg38) 2022
#方法一
mkdir chipseq_human_hg38
cd chipseq_human_hg38
nohup wget https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/latest/hg38.fa.gz > output.log &
gunzip hg38.fa.gz
nohup bowtie2-build hg38.fa hg38 > refchipseq.log &
jobs
#方法二：直接下载这个https://bowtie-bio.sourceforge.net/bowtie2/manual.shtml 别人构建好的索引，但只有hg19
mkdir human_hg19
cd human_hg19
nohup wget https://genome-idx.s3.amazonaws.com/bt/hg19.zip > output.log &
unzip hg19.zip 
rm -f hg19.zip 
#成功构建后，将得到6个索引文件

第二步 根据需要准备chipseq文件

mkdir 1rawdata 2clean 3mapping 4rmdup 5peakcalling
##写循环运行的代码
cd 1rawdata
cat >srr.list
#chipseq，需要去网上下载的文件，或者是输入自己的数据名称
SRR12515038
SRR12515039
SRR12515040
SRR12515041
#按ctrl C退出
conda activate chipseq
##1 rawdata and QC
cat srr.list |while read id;do ( nohup prefetch $id & );done
cat nohup.out
#查看是否successfully下载，并且移到同一目录下，删除空白目录
cat srr.list |while read id;do ( nohup mv /home/1rawdata/$id/$id.sra /home/1rawdata & );done
cat srr.list |while read id;do ( rm -d $id );done

ls *sra
fastq-dump --gzip --split-3 *.sra
rm *sra
fastqc *.fastq.gz
conda activate multiqc
multiqc *zip

##2 clean 过滤再QC
cd ~/2clean
#单端测序文件
for i in `ls ~/1rawdata/*.fastq.gz`; do i=${i/.fastq.gz/}; echo "trim_galore --phred33 -q 20 --length 36 --stringency 3 --fastqc -o /home/2clean ${i}.fastq.gz"; done>trim_galore.command
cat trim_galore.command
sh trim_galore.command

#双端测序文件
for i in `ls *R1.fastq.gz`; do i=${i/R1.fastq.gz/}; echo "trim_galore --phred33 -q 20 --length 36 --stringency 3 --fastqc --paired -o ~/datawpl/2clean ${i}R1.fastq.gz ${i}R2.fastq.gz"; done>trim_galore.command
cat trim_galore.command
echo "multiqc ~/rawdata/*zip -o ~/rawdata" >> trim_galore.command
cat trim_galore.command
sh trim_galore.command
cd clean

##3 mapping
cd ~/3mapping
ln -s ~/2clean/*.fq.gz ~/3mapping
bowtie2_index= /home/reference/human_hg19/hg19

#单端测序文件
for i in `ls *.fq.gz`
do
i=${i/.fq.gz/}
bowtie2 -p 10 -x /home/reference/human_hg19/hg19 -U ${i}.fq.gz -S ${i}.sam
done 
#双端测序文件
for i in `ls *1_val_1.fq.gz`
do
i=${i/1_val_1.fq.gz/}
bowtie2 -p 20 -x $bowtie2_index -1 ${i}1_val_1.fq.gz -2 ${i}2_val_2.fq.gz -S ${i}.sam
done 

ls *.sam|cut -d"." -f 1|while read id;
do
##将sam文件转为bam文件--samtools view
samtools view -@ 8 -S $id.sam -1b -o $id.bam
##将得到的bam文件进行排序，建立索引
samtools sort -n -@ 20 -o $id.sort.bam $id.bam
##remap未配对的序列（若未配对数据部分为0B大小，也可省略这一步）
samtools fixmate -m -@ 20 $id.sort.bam $id.fixmate.bam
##重新写入
samtools sort -@ 30 -o $id.positionsort.bam $id.fixmate.bam
done

##这里看比对结果，可以判断是否有其他物种污染mapping应该达到99%，可以看PCR重复量

##4 rmdup，去除PCR重复
cd ~/4rmdup
ln -s ~/3mapping/*.positionsort.bam ~/datawpl/4rmdup
ls

ls *.positionsort.bam | while read id ;
do 
if [ ! -s ${id}.rmdup.bam ];
then
echo $id
(nohup samtools markdup -@ 10 -r $id $(basename $id ".positionsort.bam").rmdup.bam & )
fi
done

#将得到的去重后的bam文件转为bigwig/bw文件------bamCoverage
for i in `ls *.rmdup.bam`
do
i=${i/.rmdup.bam/}
samtools index $i.rmdup.bam
bamCoverage -b  $i.rmdup.bam -o ~/datawpl/4rmdup/$i.bw  --binSize 10 --normalizeUsing RPKM --scaleFactor 1 --numberOfProcessors 6
done
###生成的bigwig可导入IGV可视化peaks，有个程序可以合并生物学重复
##bigWigMerge（直接合并），将多个bigwig文件合并为单个bedGraph文件（支持求和、取均值等操作）
##WiggleTools，支持对bigwig文件进行统计计算（如均值、中位数、最大值等），再转为bigwig格式
#conda install -c bioconda wiggletools

#我个人觉得先完成peakcalling,把实验组与其对应的input做对比后，再进行merge会好点
####peakscalling，是对目标ChIP样本与input样品reads信号进行比较，input用于减少背景噪音
cd ~/5peakcalling
####H3K27ac的call peak的峰为narrow peak，如果是宽峰注意调整参数
ln -s ~/4rmdup/*.rmdup.bam ~/5peakcalling
ls
#SRR12515038_trimmed.rmdup.bam  SRR12515039_trimmed.rmdup.bam  SRR12515040_trimmed.rmdup.bam  SRR12515041_trimmed.rmdup.bam
#根据实验分组来input1 input2 rad21-1 rad21-2
macs2 callpeak -B -t SRR12515040_trimmed.rmdup.bam  -c SRR12515038_trimmed.rmdup.bam -f BAM -g hs -n R --outdir /home/5peakcalling -p 0.05 --cutoff-analysis --nomodel --extsize 147
macs2 callpeak -B -t SRR12515041_trimmed.rmdup.bam  -c SRR12515039_trimmed.rmdup.bam -f BAM -g hs -n P --outdir /home/5peakcalling -p 0.05 --cutoff-analysis --nomodel --extsize 147

###按照实际结果适当调整参数，默认-p 0.05 ，可以严格点-q 0.01 或者放宽条件-q 0.05 ，如果是宽峰就需要再调整参数broad-peak
#-f 选择输入格式，BAM是bam格式。本次分析为单端测序文件，故-f选项为BAM；若数据为双端测序文件，则需指定-f选项为BAMPE
#-t 可输入多个比对文件
#-c input文件
#-g hs是人，mm是小鼠。
#-n 输出的名字
#extsize指定片段大小
#-B/--bdg：是否输出bedgraph格式的文件

#正常情况下MACS2会输出以下几个结果文件：SAMPLE_peaks.narrowPeak，SAMPLE_peaks.xls，SAMPLE_summits.bed，SAMPLE_model.r，bdg文件
#narrowPeak：窄峰调用结果文件，包含峰值位置和统计信息;BED 6+4 格式
#broadPeak：宽峰调用结果文件（使用 --broad 选项时生成）
#.xls：峰值的详细信息，包含每个峰的统计数据（峰值位置、信噪比等）。
#.bedGraph：用于可视化的信号文件，表示峰值的富集程度。
#summits.bed：BED 格式，包含 peak 的 summits 位置，第 5 列是-log10pvalue，如果想找 motif，推荐使用此文件
#call peaks是每个处理样本与input或者igG最为control之间比较（可能有时还没有input），不能用call peaks来做不同处理的样本之间差异分析。

#macs2 callpeak -B -t rmdup.bam  -c input_rmdup.bam -f BAMPE -g hs -n R --outdir /peakcalling -p 0.05 --cutoff-analysis --nomodel --extsize 147
#有narrowPeak文件
#peak注释使用ChIPseeker，安装R语言，R包。narrowPeak文件
#得到的bed文件在R中用chipseeker包进行注释，注释还可以用linux上的annotate.pl函数进行注释
#获得gain.bed和loss.bed文件再计算bw文件中peaks的信号强度，并绘制热图和信号图谱

第四步 后续分析

cd /home/6peakmak
ln -s ~/4rmdup/*.bigwig ~/6peakmak
ln -s ~/5peakcalling/*.bed ~/6peakmak
ls
#P_summits.bed  R_summits.bed  SRR12515038_trimmed.bigwig  SRR12515039_trimmed.bigwig  SRR12515040_trimmed.bigwig  SRR12515041_trimmed.bigwig
computeMatrix reference-point --missingDataAsZero --referencePoint center -p 15 -b 2000 -a 2000 -R P_summits.bed  -S SRR12515040_trimmed.bigwig --skipZeros -o matrix_center.gz --outFileSortedRegions regions_center.bed 
plotHeatmap -m matrix_center.gz -out matrix_center.heatmap.png  --colorMap Blues --zMin 0 --zMax 5 --dpi 720 --heatmapWidth 10 
display matrix_center.heatmap.png

#########获得gain.bed和loss.bed文件再计算bw文件中peaks的信号强度，并绘制热图和信号图谱
###使用bedtools 的intersect函数对call peak后的bed文件交叉获得单独的unique.bed
bedtools intersect -a R_summits.bed -b P_summits.bed -v > want_to_show.bed    ###获得在P中有，而R中没有的就是R样本中获得的peak
#-wa -wb 可以对每一次做交集运算输出传给-a和-b的bed文件的相应信息
#-v 相当于集合运算的A-B，即不与B重叠的A部分
#-c 则可以对A包含多少个B进行统计

bedtools merge [OPTION] –i <bed/gff/vcf>
#并集运算。

#bedtools功能强大，包含众多子命令。它可试图解决多种问题：基因组区间运算，多文件操作（bam或者bed文件），BEDPE文件操作，各种文件格式转换， Fasta文件操作，BAM文件操作，以及其他一些小工具。
#intersect-查看两个文件不同区间是否有overlap；
#closest-#找到和目标区间最近的特征区间；
#coverage-计算特定区间的覆盖深度；

#deeptools 提供了computeMatrix命令以计算特定基因组区域的ChIP-seq信号，并通过plotHeatmap和plotProfile函数分别产生ChIP-seq信号热图与谱图。
#使用deeptools生成ChIP-seq信号热图与谱图 - 简书 https://www.jianshu.com/p/1f38674b2475
#官方说明书：computeMatrix — deepTools 3.5.6 documentation https://deeptools.readthedocs.io/en/develop/content/tools/computeMatrix.html#deepBlue%20arguments
###可视化peaks强度
#computeMatrix提供以下两种用法以不同的参考系计算ChIP-seq的信号：scale-regions和reference-point
#computeMatrix scale-regions ：以一段区域为参考系，将所有参考的基因组区域缩放至指定的区域长度，计算该区域内的ChIP-seq信号
#computeMatrix reference-point：以某一坐标为参考系，取上下游一定长度作为窗口计算区域内的ChIP-seq信号值
computeMatrix reference-point --missingDataAsZero --referencePoint center -p 15 -b 2000 -a 2000 -R gain.bed  -S P-H3K27ac.bw R-H3K27ac.bw --skipZeros -o matrix_center.gz --outFileSortedRegions regions_center.bed 

computeMatrix reference-point -S peak1.bw peak2.bw
                            -R want_to_show.bed
                            -b 2000 -a 2000 -m 5000
                            --skipZeros -o matrix_center.gz
                            --missingDataAsZero --referencePoint center -p 15
                            --outFileSortedRegions regions_center.bed

#--scoreFileName, -S：输入的bw文件，可以是校正后的bw文件,可输入多个文件，即生物学中不同处理组
#--regionsFileName, -R：BED/GTF 格式的文件，指定需要计算的基因组区域,即需要绘制的靶标区域
#--outFileName, -out, -o：输出的文件
#--regionBodyLength, -m：缩放后的基因组区域长度（bp），指定的基因组区域都会被缩放至该长度
#--beforeRegionStartLength, -b, --upstream：基因组区域起始位点上游需要延伸的长度（bp）
#--afterRegionStartLength, -a, --downstream：基因组区域终止位点下游需要延伸的长度（bp）
#--skipZeros：设置后跳过只有0的位点。
#reference-point模式大部分参数与scale-regions模式一致，不同的是：--referencePoint：指定参考点，可选：TSS, TES, center

plotHeatmap -m matrix_center.gz -out matrix_center.heatmap.png  --colorMap Blues --zMin 0 --zMax 5 --dpi 720 --heatmapWidth 10 

#--matrixFile, -m：computeMatrix输出的结果文件
#--outFileName, -out, -o：输出的图片名称，格式可选 “png”, “eps”, “pdf”, “svg”等
#--whatToShow：输出图片中的内容，可选： “plot and heatmap”, “heatmap only”, “heatmap and colorbar”

$ plotHeatmap -m matrix.mat.gz \
      -out ExampleHeatmap1.png \
      --whatToShow “plot, heatmap and colorbar”

#只想单独生成ChIP-seq信号谱图可以使用plotProfile命令
plotProfile -m matrix.mat.gz \
              -out ExampleProfile1.png \
              --numPlotsPerRow 2 \
              --plotTitle "Test data profile"


bamCoverage -b reads.bam -o coverage.bw
#–bam, -b：要处理的 BAM 文件
#–outFileName, -o：输出文件名
#–outFileFormat, -of：输出文件的格式，可以是 bigwig、bedgraph
#–normalizeUsing：选择：RPKM，CPM，BPM，RPGC，None

第五步 motif分析

#利用ChIP-Seq结果在基因组区域中寻找富集的Motifs
#homer进行motif分析
conda install -c bioconda homer 
#MACS2找到的peaks会存放在生成的*.bed文件里。homer软件找motif需要读取我们这里得到的bed格式的peaks文件。homer软件不仅可以找到motif，还可以注释peaks。
#Motif分析输入的数据一般是ChIP-seq流程过程中MACS2 进行callpeak之后的bed（narrawPeak或broadPeak）文件
#虽然conda安装好了Homer，但并没有包含参考序列或注释数据 。所以需要使用 configureHomer.pl下载各种数据包。

#使用configureHomer.pl完成Homer软件的配置
#Homer Software and Data Download http://homer.ucsd.edu/homer/motif/
#使用Homer进行motif分析需要先安装下载对应的基因组，安装基因组需要用到configureHomer.pl。
wget http://homer.salk.edu/homer/configureHomer.pl
# Installing the basic HOMER software
perl configureHomer.pl -install
# Download the hg19 version of the human genome
perl configureHomer.pl -install hg19
# 因为课题组人和小鼠的研究居多，这里只以这两者举例
perl configureHomer.pl -install mm10

perl /home/miniconda3/envs/chipseq/share/homer/.//configureHomer.pl -install hg19
#使用 Ctrl+Z 挂起进程 ，再通过 bg 后台运行，或 fg 将其调回前台。
#用这个下载，只有一个homer.package.zip，而且是在当前的文件夹
unzip homer.package.zip
#/home/miniconda3/envs/chipseq/share/homer/data
mv /home/6homer/data/genomes /home/miniconda3/envs/chipseq/share/homer/data
#安装完成后，HOMER 会在 share/homer/data/genomes/hg19/ 目录存储基因组数据，后续分析可直接调用 hg19 标识符。

mkdir ~/6homer
cd ~/6homer
ln -s  ~/5peakcalling/*.bed ~/6homer
ls
#P_summits.bed  R_summits.bed

## 提取对应的列给HOMER作为输入文件，格式如下 MACS_peak_1  chr1    1454086 1454256 +
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' P_summits.bed >homer_peaks.bed
head -n 5 homer_peaks.bed
findMotifsGenome.pl homer_peaks.bed /home/wangpeili/miniconda3/envs/chipseq/share/homer/data/genomes/hg19 ~/datawpl/6homer -size 150
#输入文件：处理好的 Peak/Positions file,bed文件
# 基因组版本：hg19，这里最好指定绝对路径，不然会报错
# 输出路径：out_motifDir
#-len：motif 大小设置，默认 8,10,12；越大需要的计算资源越多
#-size: 用于 motif 寻找得片段大小，默认 200bp；越大需要得计算资源越多
#-S：结果输出多少 motifs, 默认 25,motif的软件默认长度为8，10，12
#-mask: 使用repeated-mask序列

多组合绘图

#H3K4me1（组蛋白 H3 第 4 位赖氨酸单甲基化）：功能：通常与增强子区域相关联，是活性增强子标志之一。作用：标志潜在的增强子区域，参与调控基因的转录起始。
#H3K4me3（组蛋白 H3 第 4 位赖氨酸三甲基化）：功能：与启动子区域相关联，是活跃基因的启动子标志。作用：参与转录起始复合物的招募，促进基因转录。
#H3K27me3（组蛋白 H3 第 27 位赖氨酸三甲基化）：功能：与基因沉默相关联，是抑制性标志。作用：标志沉默基因区域，通过聚合酶抑制和染色质结构改变抑制基因表达
#H3K27ac（组蛋白 H3 第 27 位赖氨酸乙酰化）：功能：与活性增强子和启动子相关联。作用：标志活跃的增强子和启动子区域，促进开放染色质结构和基因表达。
#H3K36me3（组蛋白 H3 第 36 位赖氨酸三甲基化）：功能：与转录延伸相关联。作用：标志转录活跃区域，参与 RNA 聚合酶 II 的延伸过程。
#H3K9me3（组蛋白 H3 第 9 位赖氨酸三甲基化）：功能：与异染色质形成相关联，是沉默性标志。作用：参与异染色质形成，抑制基因表达，维持基因组稳定性。

bedtools intersect -a EtOH_S100A8.bed -b H3K4me3_peaks.narrowPeak > promoter.bed
bedtools intersect -a EtOH_S100A8.bed -b H3K4me1_peaks.narrowPeak > enhancer.bed
computeMatrix reference-point \ 
       --referencePoint center \
       -b 2000 -a 2000 \
       -R promoter.bed enhancer.bed \
       -S EtOH_S100A8.bw EtOH_H3K4me1.bw EtOH_H3K4me3.bw \
       --skipZeros \
       -o matrix1_TSS.gz \
       --outFileSortedRegions regions1_H3K4me3_l2r_genes.bed

plotHeatmap -m matrix1_TSS.gz \
     -out ExampleHeatmap.png \
     --whatToShow 'heatmap and colorbar'