CHIP-Seq数据分析流程
文章目录
- CHIP-Seq数据处理流程
- 安装相关软件工具
-
获取原始SRA序列
-
将SRA序列转换为fastq格式文件
-
数据质量控制与初步筛选：随后对高通量测序数据进行质量控制和初步筛选
- 比对
-
- 下载小鼠参考基因组的索引和注释文件
- 比对
-
将SAM格式转换为BAM格式
对BAM文件建立索引
导入IGV软件查看
通过MACS工具识别高丰度峰
安装MACS软件包
* 结果注释与可视化
* * 使用deeptools进行可视化
* peaks注释CHIP-Seq数据分析流程
相关软件安装
所需的软件包括sratoolkit、fastqc、bowtie2、samtools、macs2、htseq-count以及bedtools等，并且这些功能与前面RNA-seq数据分析中常用的软件非常相似。因此已经完成安装工作。后续遇到未预先安装的则会进行补充安装。
###igv
axel -n 20 https://data.broadinstitute.org/igv/projects/downloads/2.8/IGV_Linux_2.8.12_WithJava.zip
unzip IGV_Linux_2.8.12_WithJava.zip
cd IGV_Linux_2.8.12/
./igv.sh
代码解读数据下载
还是利用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
do
echo $i
fastq-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
do
arr=(${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 times
7170824 (57.12%) aligned exactly 1 time
2233449 (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 ;do
samtools view -@ 8 -S $id.sam -1b -o $id.bam
samtools sort -@ 8 -l 5 -o $id.sort.bam $id.bam
done
代码解读为bam文件建立索引
ls *.sort.bam |xargs -i samtools index {}
代码解读载入IGV查看
通过官网上获取IGV软件包后，默认无需解压即可直接运行。
首先导入参考基因组；可以选择导入自己已下载的版本或使用IGV自带的标准参考基因组库中的Ref_Genome（需以FASTA格式存储）。
导入待比对的文件时，请确保该文件已事先经过排序（sort）和索引（index）处理。
比较结果显示，在对照组中未出现的特征点则被视为差异点。
用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-Sequencing
positional 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：Extendable Bed Data。它主要用于记录各峰值点的位置信息，并也即标记各峰值的位置坐标。此外，在进行结合位点分析时，MACS推荐使用此文件来识别结合位点对应的序列模式。值得注意的是，该BED文件可以直接导入UCSC BedTools浏览器进行直观查看，在其他软件处理时则需注意跳过首行数据以避免干扰。
该Excel文件以表格形式存储着peak数据。尽管其后缀为xls，在文本编辑器中即可打开，并与BED格式相似。需要注意的是，在此文件中坐标采用1为基线（即基于1编码），而BED文件则基于0编码（即BED files are 0-based）。关键点是将xls中的坐标值减去1即可得到对应的BED位置信息。
NAME_peaks中的narrowPeak和broadPeak具有相似性。以下列举了四个指标：显示分数（integer score for display）、变化倍数（fold-change）、−log₁₀ p值（−log10pvalue）、−log₁₀ q值（−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的功能分为三类: ● 标注：用于标注与基因最邻近的peaks及其所在位置 ● 分析：评估重叠区域在ChIP peak数据集中的统计显著性，并整合GEO数据库以促进当前发现与已有研究的有效对比 ● 可视化: 展示与peaks相关的染色体区域覆盖情况;通过整合TSS和peaks进行差异表达谱分析及热图展示;提供基因组注释功能;分析TSS距离特征以及peaks与基因之间的交集情况
使用deeptools进行可视化
deeptools采用bamCoverage和bamCompare功能模块实现数据格式转换功能。该系统为了对不同实验样本进行分析对比需求，在处理流程中首先将基因组序列均分为宽度相等的区间，并计算每个区间内的读取计数值。最终汇总生成描述性统计数据。对于多个样本而言其比值可表示为两个样本之间的相对比例关系而当处理单一实验样本时则可采用平均值标准差（SES）法进行标准化处理。
##deeptools安装
git clone https://github.com/fidelram/deepTools
cd deepTools
python setup.py install
代码解读bamCoverage的基本用法