从GEO下载数据

下载地址https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121424
以第一个样本GSM3436194为例,点进链接后最下面有SRA链接(第二张图 ARX4901015),再点进链接,看到SRR8073138编号。

Chip-seq 数据分析流程 - 图1
Chip-seq 数据分析流程 - 图2
Chip-seq 数据分析流程 - 图3
下载之前我们先写个文件记录要下载的文件与编号的对应关系,后面用的到。文件格式如下

  1. (base) [longfei@localhost GSE121424]$ cat info
  2. SRR8073138 WT_DMSO_Input
  3. SRR8073141 WT_DMSO_H3K4ME2
  4. SRR8073142 WT_DMSO_H3K27AC
  5. SRR8073143 WT_GSK_Input
  6. SRR8073146 WT_GSK_H3K4ME2
  7. SRR8073147 WT_GSK_H3K27AC

下载工具有多个,利用用R包GEOquery,使用aspera软件等。我们用最简单的方法,ncbi官方软件SRA Toolkit,请自行安装。

  1. $ cat info | while read -a line;do (nohup prefetch ${line[0]} -O sra &);done

根据SRR编号循环后台下载到新建的sra文件夹中。

  1. ls sra
  2. SRR8073138.sra
  3. SRR8073141.sra
  4. SRR8073142.sra
  5. SRR8073143.sra
  6. SRR8073146.sra
  7. SRR8073147.sra

提取fastq文件

同样使用SRA Toolkit 工具,提取到fq 文件夹中。

  1. cat info | while read -a line;do (nohup fastq-dump --gzip --split-3 -A srr/${line[0]}.sra -O fq &);done
  1. ls fq
  2. SRR8073138.sra_1.fastq.gz SRR8073141.sra_1.fastq.gz SRR8073142.sra_1.fastq.gz SRR8073143.sra_1.fastq.gz SRR8073146.sra_1.fastq.gz SRR8073147.sra_1.fastq.gz
  3. SRR8073138.sra_2.fastq.gz SRR8073141.sra_2.fastq.gz SRR8073142.sra_2.fastq.gz SRR8073143.sra_2.fastq.gz SRR8073146.sra_2.fastq.gz SRR8073147.sra_2.fastq.gz
  4. SRR8073138.sra.fastq.gz SRR8073141.sra.fastq.gz SRR8073142.sra.fastq.gz SRR8073143.sra.fastq.gz SRR8073146.sra.fastq.gz SRR8073147.sra.fastq.gz

fastqc 质控

我只挑了其中的一个fastq文件看了下,一般都是处理好的,包括去接头和去除低质量的碱基。

  1. fastqc fq/SRR8073138.sra_1.fastq.gz -o ./ --noextract

生成一个html文件和一个压缩包,传到自己的电脑上,打开html检查质量。

  1. SRR8073138.sra_1_fastqc.html SRR8073138.sra_1_fastqc.zip

Chip-seq 数据分析流程 - 图4

将fastq文件比对到基因组上

此处使用软件bowtie2,提前下载好参考基因组文件,这里用hg19。

  1. bowtie_ref=homo_ref/bowtie/hg19
  2. mkdir sam
  3. cat info | while read -a line;do ( nohup bowtie2 -x $bowtie_ref -p 4 -1 fq/${line[0]}.sra_1.fastq.gz -2 fq/${line[0]}.sra_2.fastq.gz -S sam/${line[1]}.sam & );done

将生成的sam文件通过排序变成bam文件,使用软件samtools。这两步可以通过管道合并成一步,省去生成sam文件占用空间浪费时间,但是我经常在这里遇到问题,所以分成两步。

  1. mkdir bam
  2. ls sam |while read sam;do (nohup samtools sort -O BAM -@ 4 -o bam/${sam%.*}.bam sam/$sam & );done

去除重复序列

使用软件Sambamba,可以查看我的另一篇文章Sambamba 去除重复工具

  1. mkdir rmdup
  2. ls bam | while read bam;do sambamba markdup -r -t 4 bam/${bam} rmdup/${bam%.*}.rmdup.bam ;done

使用 bamTools 工具质控

详细命令含义及结果含义移步https://www.jianshu.com/p/2fddb062c503
主要是看ChiP是否合格,想要的地方是否富集到了,区别于开头的测序数据质控。

plotCoverage

查看测序深度相关性

  1. mkdir qc
  2. plotCoverage -b rmdup/*bam \
  3. --plotFile qc \
  4. -n 10000000 \
  5. --plotTitle "example_coverage" \
  6. --outRawCounts coverage.tab \
  7. --ignoreDuplicates \
  8. --minMappingQuality 10

plotFingerprint

比较input和实验组,如果差距不大,说明的富集的不好。

  1. plotFingerprint \
  2. -b testFiles/*bam \
  3. --labels H3K27me3 H3K4me1 H3K4me3 H3K9me3 input \
  4. --minMappingQuality 30 \
  5. --skipZeros \
  6. --region 19 --numberOfSamples 50000 \
  7. -T "Fingerprints of different samples" \
  8. --plotFile fingerprints.png \
  9. --outRawCounts fingerprints.tab

—numberOfSamples 参数主要是随机从样本中抽取的箱数来计算相对的覆盖度

使用 deepTools 工具生成bw文件可视化

详细命令含义及结果含义移步 https://www.jianshu.com/p/e7e2c65183fd

bamCoverage

此命令用于单个样本标准化生成bw文件,使得样本之间可以互相比较。bw文件可用IGV软件可视化。

  1. mkdir bw
  2. ls rmdup/ | grep -v bai | while read bam ;do \
  3. (nohup bamCoverage --bam rmdup/$bam \
  4. -o bw/${bam%*rmdup}.bw \
  5. --binSize 10 \
  6. --normalizeUsing RPGC \
  7. --effectiveGenomeSize 2685511504 &);done

computeMatrix

将多个样本放入一个矩阵中,输入为上面生成的bw文件,为后续可视化作图做准备,如plotHeatmap,plotProfile作图。

  1. mkdir result
  2. nohup computeMatrix reference-point \
  3. --referencePoint TSS \
  4. -b 3000 -a 10000 \
  5. -R homo_ref/hg19.bed \
  6. -S bw/H3K4Me2_DMSO_a.bw bw/H3K4Me2_DMSO_b.bw bw/H3K4Me2_OG86_a.bw bw/H3K4Me2_OG86_b.bw \
  7. -o result/H3K4me2_TSS.gz \
  8. --outFileSortedRegions result/H3K4ME2_gene.bed \
  9. --skipZeros \
  10. -p 10 &

plotProfile

使用plotProfile生成gene在TSS区的富集图

  1. nohup plotProfile -m H3K4me2_TSS.gz -out H3K4me2_TSS.png --perGroup &

Chip-seq 数据分析流程 - 图5

使用 MACS2 call peak

参考
https://github.com/taoliu/MACS/
https://www.jianshu.com/p/0c272643f88b
https://www.jianshu.com/p/6a975f0ea65a

  1. nohup macs2 callpeak -t rmdup/H3K9me2K.rmdup.bam -c rmdup/Input-K.rmdup.bam -f BAMPE -g hs -n H3K9me2K -q 0.05 --nomodel &