文章目录

CHIP-Seq数据分析流程
- 相关软件安装
- 数据下载
- 将SRA转化为fastq文件
- 数据质控，过滤
- - 数据质控
  - 然后进行质控，过滤低质量reads，去接头
- 比对
- - 下载小鼠参考基因组的索引和注释文件
  - 比对
- sam文件转bam
- 为bam文件建立索引
- 载入IGV查看
- 用MACS call peak
- - 安装MACS
- 结果注释与可视化
- - 使用deeptools进行可视化
  - peaks注释

CHIP-Seq数据分析流程

数据下载

还是利用sratoolkit的prefetch下载，首先构建SRR_Acc_List.txt

for ((i=204;i<=209;i++))
do
echo SRR620$i >> SRR_Acc_List.txt
done
###文件内容如下
SRR620204
SRR620205
SRR620206
SRR620207
SRR620208
SRR620209
##接下来进行下载
wkd=/home/meiling/baiduyundisk/Chip-seq #设置工作目录
cat SRR_Acc_List.txt | while read id
do
nohup prefetch  ${id} &
done

将SRA转化为fastq文件

wkd=/home/meiling/baiduyundisk/Chip-seq #设置工作目录for i in $wkd/rawdata/*sra
doecho $ifastq-dump --split-3 --skip-technical --clip --gzip $i  ## 批量转换
done

数据质控，过滤

数据质控

fastqc生成质控报告，multiqc将各个样本的质控报告整合为一个

wkd=/home/meiling/baiduyundisk/Chip-seq #设置工作目录ls $wkd/rawdata/*gz | xargs fastqc -t 5
multiqc ./
##得到以下结果
ls -lh multiqc_data/
total 7.1M
-rw-r--r-- 1 meiling meiling 312K Nov 10 14:10 multiqc_data.json
-rw-r--r-- 1 meiling meiling 1.5K Nov 10 14:10 multiqc_fastqc.txt
-rw-r--r-- 1 meiling meiling  621 Nov 10 14:10 multiqc_general_stats.txt
-rw-r--r-- 1 meiling meiling  24K Nov 10 14:10 multiqc.log
-rw-r--r-- 1 meiling meiling 1.2M Nov 10 14:10 multiqc_report.html
-rw-r--r-- 1 meiling meiling  604 Nov 10 14:10 multiqc_sources.txt
-rw-r--r-- 1 meiling meiling 576K Nov 10 14:08 SRR620204_fastqc.html
-rw-r--r-- 1 meiling meiling 431K Nov 10 14:08 SRR620204_fastqc.zip

每个id_fastqc.html都是一个质量报告，multiqc_report.html是所有样本的整合报告

然后进行质控，过滤低质量reads，去接头

wkd=/home/meiling/baiduyundisk/Chip-seq #设置工作目录
mkdir $wkd/cleandata
cd $wkd/cleandata
ls $wkd/rawdata/*.fastq.gz >config
##开始质控
wkd=/home/meiling/baiduyundisk/Chip-seq/cleandata #设置工作目录
cd $wkd
cat config |while read id
doarr=(${id})fq1=${arr[0]}
#        fq2=${arr[1]} trim_galore -q 25 --phred33 --length 36 --stringency 3 -o $wkd  $fq1
done

比对

下载小鼠参考基因组的索引和注释文件

还是选择在ensemble上下载小鼠参考基因组

选择primary进行下载

axel -n 20 ftp://ftp.ensembl.org/pub/release-101/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz
axel -n 20 ftp://ftp.ensembl.org/pub/release-101/gtf/mus_musculus/Mus_musculus.GRCm38.101.chr_patch_hapl_scaff.gtf.gz
##解压
gzip -d Mus_musculus.GRCm38.dna.primary_assembly.fa.gz
ls -lh
total 2.7G
-rw-r--r-- 1 meiling meiling  33M Nov  9 21:32 Mus_musculus.GRCm38.101.chr_patch_hapl_scaff.gtf.gz
-rw-r--r-- 1 meiling meiling 2.6G Nov  9 21:31 Mus_musculus.GRCm38.dna.primary_assembly.fa
#构建索引
nohup bowtie2-build ../Mus_musculus.GRCm38.dna.primary_assembly.fa mouse  &

比对

这里选择bowtie,或者bwa
将得到的fastq文件用bowtie2比对小鼠参考基因组上

wkd=/home/meiling/baiduyundisk/Chip-seq #设置工作目录
cd $wkd/cleandata
ls *gz|cut -d"_" -f 1 |sort -u |while read id;do
ls -lh ${id}_trimmed.fq.gz
bowtie2 -q -p 10 -x $wkd/reference/mouse/bowtie_index/mouse -U ${id}_trimmed.fq.gz -S $wkd/align/${id}.sam
done

比对结果输出

bash align.sh >align.log 2>&1
#打开align.log
-rw-r--r-- 1 meiling meiling 517M Nov 10 14:15 SRR620204_trimmed.fq.gz
12553187 reads; of these:12553187 (100.00%) were unpaired; of these:3148914 (25.08%) aligned 0 times7170824 (57.12%) aligned exactly 1 time2233449 (17.79%) aligned >1 times
74.92% overall alignment rate

sam文件转bam

wkd=/home/meiling/baiduyundisk/Chip-seq #设置工作目录
cd $wkd/align
ls *.sam|cut -d"." -f 1 |while read id ;dosamtools view -@ 8 -S $id.sam -1b -o $id.bamsamtools sort -@ 8 -l 5 -o $id.sort.bam $id.bam
done

为bam文件建立索引

ls *.sort.bam |xargs -i samtools index {}

载入IGV查看

从官网下载IGV, 解压即可使用，linux下 igv.sh打开IGV界面
首先载入参考基因组，可以载入自己下载好的参考基因组，也可选择IGV中含有的参考基因组，ref_Genome 必须是fasta格式。
载入比对的文件，比对的文件必须先经过sort 和 index, 才可加载。
比对可视化结果：其他组有的峰在对照组没有，即为peak.

用MACS call peak

安装MACS

#下载安装
axel -n 20 https://files.pythonhosted.org/packages/e2/61/85d30ecdd34525113e28cb0c5a9f393f93578165f8d848be5925c0208dfb/MACS2-2.2.7.1.tar.gz
tar -xvf MACS2-2.2.7.1.tar.gz
cd MACS2-2.2.7.1/
##安装
python setup.py install
##验证
macs2 -h
##出现以下结果
usage: macs2 [-h] [--version]{callpeak,bdgpeakcall,bdgbroadcall,bdgcmp,bdgopt,cmbreps,bdgdiff,filterdup,predictd,pileup,randsample,refinepeak}...macs2 -- Model-based Analysis for ChIP-Sequencingpositional arguments:{callpeak,bdgpeakcall,bdgbroadcall,bdgcmp,bdgopt,cmbreps,bdgdiff,filterdup,predictd,pileup,randsample,refinepeak}
##call peak
wkd=/home/meiling/baiduyundisk/Chip-seq #设置工作目录
cd $wkd/align
ls *.sort.bam|cut -d"." -f 1 |sort -u |while read id
do
macs2 callpeak -c SRR620208.bam -t ${id}.bam -q 0.05 -f BAM -g mm -n ${id} 2> ${id}.log &
done
# （在当前目录下）统计 *bed 的行数（peak数）
wc -l *bed
7765 SRR620204_summits.bed
1963 SRR620205_summits.bed
6775 SRR620206_summits.bed
0 SRR620207_summits.bed
0 SRR620208_summits.bed##对照组
0 SRR620209_summits.bed##对照组
16503 total
#RYBP SRR620207 无数据，估计是数据上传错误，可以下载作者上传的peak数据。 ● 下载RYBP的peak数据
wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE42nnn/GSE42466/suppl/GSE42466_RYBP_peaks_5.txt.gz
gzip -d GSE42466_RYBP_peaks_5.txt.gz
mv GSE42466_RYBP_peaks_5.txt SRR620207_summits.bed

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 文件。

结果注释与可视化

结果的注释用的是Y叔的 Chipseeker包。

ChIPseeker的功能分为三类: ● 注释：提取peak附近最近的基因，注释peak所在区域 ● 比较：估计ChIP peak数据集中重叠部分的显著性；整合GEO数据集，以便于将当前结果和已知结果比较 ● 可视化： peak的覆盖情况；TSS区域结合的peak的平均表达谱和热图；基因组注释；TSS距离；peak和基因的重叠。
后续在R语言中实现。

使用deeptools进行可视化

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

##deeptools安装
git clone https://github.com/fidelram/deepTools
cd deepTools
python setup.py install

bamCoverage的基本用法

peaks注释