ChIP-seq分析
创建小环境
conda create -n ChIP python=3
安装软件
conda install -y sra-tools
conda install -y trim-galore samtools
conda install -y deeptools homer meme
conda install -y macs2 bowtie bowtie2 bwa
数据下载
prefetch SRR8476707 -O ~
数据转换，sra格式转换为fastq格式
for i in *sra
do
echo $i
/path/sratoolkit.2.3.5-2-mac64/bin/fastq-dump --split-3 $i
done
|单个命令如下
fastq-dump --gzip --split-3 -O ./ ./SRR8476707.sra
过滤，用trim-galore软件，直接fastqc
因为文件多，所以写个脚本直接运行
vim trim_galore_batch.sh
输入以下内容
#!/bin/bash
# This is for trimming a batch data
for i in SRR8476709 SRR8476710 SRR8476717 SRR8476718 SRR8476719 SRR8476720 SRR8476721
do
trim_galore -q 25 --phred33 --length 25 --stringency 4 --gzip --paired ~/Data/H3K9me2-ChIP_seq_VC/SRR/fastq/${i}_1.fastq.gz ~/Data/H3K9me2-ChIP_seq_VC/SRR/fastq/${i}_2.fastq.gz -o ~/Data/H3K9me2-ChIP_seq_VC/clean_2 --fastqc
done
运行
bash trim_galore_batch.sh
建索引
bwa index -a bwtsw ~/genome/mouse/ensembl/Mus_musculus.GRCm39.dna.toplevel.fa
| 索引和gtf&gfa文件建立在一个文件夹下面，会生成5个文件，分别是.amb .ann .bwt .pac .sa
比对
比对一般用BWA或者Bowtie2
|也可以写成文件，运行文件，当时一直报错，就一个一个比对的
加上nohup ***** &，让代码在后台运行
nohup bwa mem -v 3 -t 4 ~/genome/mouse/ensembl/Mus_musculus.GRCm39.dna.toplevel.fa ~/Data/H3K9me2-ChIP_seq_VC/fastq_p/SRR8476710_1_val_1.fq ~/Data/H3K9me2-ChIP_seq_VC/fastq_p/SRR8476710_2_val_2.fq -o ~/Data/H3K9me2-ChIP_seq_VC/bwa/SRR8476710_bwa.sam &
转换为bam文件及排序，或者有时候需要合并bam文件（将多个重复合并到一个）
|因为bwa生成的是sam文件，所以要转换成bam文件以及进行排序（如果不排序后面用macs2 peak calling的时候可能会报错）
vim bam.sh
输入以下内容
for i in SRR8476707 SRR8476708 SRR8476709 SRR8476710 SRR8476711 SRR8476717 SRR8476718 SRR8476719 SRR8476720 SRR8476721
do
samtools view -bS -h ~/Data/H3K9me2-ChIP_seq_VC/bwa/${i}_bwa.sam -o ~/Data/H3K9me2-ChIP_seq_VC/bam/${i}_bwa.bam
查看比对成功率
samtools flagstat -@ 8 SRR8476707.sam
转为bam文件
vim bam.sh
输入以下内容
do
samtools view -bS -h ~/Data/H3K9me2-ChIP_seq_VC/bwa_NCBI/${i}_bwa.sam -o ~/Data/H3K9me2-ChIP_seq_VC/bam_NCBI/${i}_bwa.bam
done
bash bam.sh
对bam文件进行排序
vim bam.sort.sh
for i in SRR8476707 SRR8476708 SRR8476709 SRR8476710 SRR8476711 SRR8476717 SRR8476718 SRR8476719 SRR8476720 SRR8476721
do
samtools sort -@ 5 -o ~/Data/H3K9me2-ChIP_seq_VC/bam_NCBI/${i}_bwa.sort.bam ~/Data/H3K9me2-ChIP_seq_VC/bam_NCBI/${i}_bwa.bam
done
转换为bw模式，导入IGV查看
安装deeptools
conda install -c bioconda/label/cf201901 deeptools
在sorted bam文件下进行转换
ls *.bam |xargs -i samtools index {}
ls *.bam |while read id;do \
bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.bw & \
done
peakcalling
vim macs2_2.sh
输入以下内容
for i in SRR8476717 SRR8476718 SRR8476719 SRR8476720 SRR8476721
do
macs2 callpeak -f BAM -t ${i}_bwa.sort.bam -n ${i}_bwa.sort.bam -g mm --outdir ~/Data/H3K9me2-ChIP_seq_VC/mcs2_2 --bdg -q 0.05
done
运行代码
nohup bash macs2_2.sh>macs_2.log &
看下找到了多少个peak
wc -l *bed
有3个生物学重复的组蛋白修饰ChIP-Seq数据怎样合并比较好呢？文献中目前有2种方法：（1）直接合并bam文件，再call peak；（2）3个重复分别call peak，然后取至少在2个样品中同时出现的 peak 进行后续的分析。这两种方法各有哪些优缺点？哪种比较科学呢？期待老师的解答，谢谢！
callpeak会得到如下结果文件：
NAME_summits.bed：Browser Extensible Data，记录每个peak的peak summits，话句话说就是记录极值点的位置。MACS建议用该文件寻找结合位点的motif。能够直接载入UCSC browser，用其他软件分析时需要去掉第一行。
NAME_peaks.xls：以表格形式存放peak信息，虽然后缀是xls，但其实能用文本编辑器打开，和bed格式类似，但是以1为基，而bed文件是以0为基.也就是说xls的坐标都要减一才是bed文件的坐标。
NAME_peaks.narrowPeak，NAME_peaks.broadPeak 类似。后面4列表示为， integer score for display， fold-change，-log10pvalue，-log10qvalue，relative summit position to peak start。内容和NAME_peaks.xls基本一致，适合用于导入R进行分析。
NAME_model.r：能通过$ Rscript NAME_model.r作图，得到是基于你提供数据的peak模型。
.bdg：能够用 UCSC genome browser 转换成更小的 bigWig 文件。
————————————————
版权声明：本文为CSDN博主「fengzhibifang」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/weixin_41745858/article/details/82628775
