ChIP-Seq数据分析上下游打通 / 开普饭

连续两次求贤令：曾经我给你带来了十万用户，但现在祝你倒闭，以及生信技能树知识整理实习生招募，让我走大运结识了几位优秀小伙伴！大家开始根据我的ngs组学视频进行一系列公共数据集分析实战，其中一个小伙伴让我非常惊喜，不需要怎么沟通和指导，就默默的完成了一个实战！

下面是温州医科大硕士“何物昂”的笔记

数据

数据使用GSE159370数据集中的SRF的ChIP-seq数据，分为serum starved (0.3% FBS) condition and 45 min after serum stimulation (15% FBS). control使用Input. （进入EBI数据库很容易拿到其地址）

数据下载使用迅雷（就4个fq文件）

## 分别为
### NIH-SRF_03FBS
### NIH-SRF_15FBS
### NIH-input_03FBS
### NIH-input_15FBS

ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR128/008/SRR12807908/SRR12807908.fastq.gz
ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR128/009/SRR12807909/SRR12807909.fastq.gz
ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR128/004/SRR12807904/SRR12807904.fastq.gz
ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR128/005/SRR12807905/SRR12807905.fastq.gz

## 过滤低质量reads和去除接头
mkdir clean
ls raw/*gz | while read id;do nohup trim_galore -q 25 --phred33 --length 36 --stringency 3 -o clean $id & done;

也可以使用aspera下载，具体参考 使用aspera从EBI下载fastq数据，抛弃NCBI的SRA数据库

比对

使用bowtie2进行比对，一个很简单的循环就可以批量处理啦！

值的注意的是小鼠参考基因组文件的bowtie2的index制作这里省略了哦：ref/bowtie2_idx/mm10

cd clean
mkdir align
ls *gz | while read id;do echo "bowtie2 -p 6 -x ref/bowtie2_idx/mm10 -U $id | samtools sort -O bam -o ../align/${id/_trimmed.fq.gz/.bam} &"; done;

## 去除pcr重复
cd ../align
ls *bam | while read id; do nohup samtools rmdup --output-fmt BAM -s $id $(basename $id ".bam").rmdup.bam & done

ll *.bam -h | cut -d " " -f 5-
2.5G Jul 7 20:51 align/NIH-input_03FBS.bam
2.3G Jul 6 19:51 align/NIH-input_03FBS.rmdup.bam
2.4G Jul 7 20:51 align/NIH-input_15FBS.bam
2.2G Jul 6 19:50 align/NIH-input_15FBS.rmdup.bam
1.5G Jul 7 20:44 align/NIH-SRF_03FBS.bam
642M Jul 6 19:46 align/NIH-SRF_03FBS.rmdup.bam
1.4G Jul 7 20:43 align/NIH-SRF_15FBS.bam
594M Jul 6 19:46 align/NIH-SRF_15FBS.rmdup.bam

call peak

peak calling 用MACS3。

软件安装：

pip install MACS3

运行：

cd ..
macs3 callpeak -t align/NIH-SRF_03FBS.rmdup.bam -c align/NIH-input_03FBS.rmdup.bam -f BAM -g mm -n callpeak/SRF_03FBS -B -q 0.05 &
macs3 callpeak -t align/NIH-SRF_15FBS.rmdup.bam -c align/NIH-input_15FBS.rmdup.bam -f BAM -g mm -n callpeak/SRF_15FBS -B -q 0.05 &

ll callpeak/* | cut -d " " -f 5-
2595368744 Jul 7 09:40 callpeak/SRF_03FBS_control_lambda.bdg
25167 Jul 6 12:58 callpeak/SRF_03FBS_model.pdf
99100 Jul 7 09:34 callpeak/SRF_03FBS_model.r
76187 Jul 7 09:40 callpeak/SRF_03FBS_peaks.narrowPeak
83419 Jul 7 09:40 callpeak/SRF_03FBS_peaks.xls
2595368744 Jul 7 09:40 callpeak/SRF_03FBS_rmdup_control_lambda.bdg
99106 Jul 7 09:33 callpeak/SRF_03FBS_rmdup_model.r
81629 Jul 7 09:40 callpeak/SRF_03FBS_rmdup_peaks.narrowPeak
88897 Jul 7 09:40 callpeak/SRF_03FBS_rmdup_peaks.xls
59556 Jul 7 09:40 callpeak/SRF_03FBS_rmdup_summits.bed
833082196 Jul 7 09:40 callpeak/SRF_03FBS_rmdup_treat_pileup.bdg
54114 Jul 7 09:40 callpeak/SRF_03FBS_summits.bed
833082196 Jul 7 09:40 callpeak/SRF_03FBS_treat_pileup.bdg
2489970108 Jul 7 09:41 callpeak/SRF_15FBS_control_lambda.bdg
24943 Jul 6 13:03 callpeak/SRF_15FBS_model.pdf
99106 Jul 7 09:34 callpeak/SRF_15FBS_model.r
155699 Jul 7 09:41 callpeak/SRF_15FBS_peaks.narrowPeak
169193 Jul 7 09:41 callpeak/SRF_15FBS_peaks.xls
2489970108 Jul 7 09:40 callpeak/SRF_15FBS_rmdup_control_lambda.bdg
99112 Jul 7 09:34 callpeak/SRF_15FBS_rmdup_model.r
166814 Jul 7 09:40 callpeak/SRF_15FBS_rmdup_peaks.narrowPeak
180344 Jul 7 09:40 callpeak/SRF_15FBS_rmdup_peaks.xls
122318 Jul 7 09:40 callpeak/SRF_15FBS_rmdup_summits.bed
757295072 Jul 7 09:40 callpeak/SRF_15FBS_rmdup_treat_pileup.bdg
111193 Jul 7 09:41 callpeak/SRF_15FBS_summits.bed
757295072 Jul 7 09:41 callpeak/SRF_15FBS_treat_pileup.bdg

wc -l callpeak/*.bed
907 callpeak/SRF_03FBS_rmdup_summits.bed
907 callpeak/SRF_03FBS_summits.bed
1852 callpeak/SRF_15FBS_rmdup_summits.bed
1852 callpeak/SRF_15FBS_summits.bed
5518 total

参数：

-t chip-seq treatment
-c chip-seq control

-f 输入BAM格式
-g 基因组size， mm表示mouse基因组size

输出文件格式参考https://github.com/macs3-project/MACS/blob/master/docs/callpeak.md

Peak 注释

peak注释使用ChIPseeker，安装：

BiocManager::install("ChIPseeker")library(TxDb.Mmusculus.UCSC.mm10.knownGene)

txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene

srf_03fbs.peak <- "./callpeak/SRF_03FBS_peaks.narrowPeak"
srf_15fbs.peak <- "./callpeak/SRF_15FBS_peaks.narrowPeak"

srf_03fbs.peak.anno <- annotatePeak(srf_03fbs.peak, tssRegion=c(-3000, 3000), TxDb=txdb)
srf_15fbs.peak.anno <- annotatePeak(srf_15fbs.peak, tssRegion=c(-3000, 3000), TxDb=txdb)

plotAnnoPie(srf_03fbs.peak.anno)
plotAnnoPie(srf_15fbs.peak.anno)

不知道是数据问题，还是SRF binding能力较弱。它在promoter区域的结合强度有点低吧

promoter <- getPromoters(TxDb=txdb,
upstream=3000, downstream=3000)

tagMatrix <- getTagMatrix(srf_03fbs.peak,
windows=promoter)

plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
xlab="Genomic Region (5'->3')",
ylab = "Read Count Frequency")

tagHeatmap(tagMatrix, xlim=c(-3000, 3000),
color="red")

同样的，也是批处理：

## peak注释详情
srf_03fbs.peak.anno.df <- as_tibble(srf_03fbs.peak.anno)

head(srf_03fbs.peak.anno.df)

# A tibble: 6 x 21
seqnames start end width strand V4 V5 V6 V7 V8 V9 V10 annotation geneChr geneStart geneEnd geneLength geneStrand geneId transcriptId distanceToTSS
<fct> <int> <int> <int> <fct> <chr> <int> <chr> <dbl> <dbl> <dbl> <int> <chr> <int> <int> <int> <int> <int> <chr> <chr> <dbl>
1 chr1 5.78e6 5.78e6 64 * callpe… 43 . 6.18 8.80 4.39 32 Distal Interg… 1 5913707 5.92e6 3692 2 226304 ENSMUST0000… 140234
2 chr1 1.04e7 1.04e7 64 * callpe… 75 . 7.76 12.1 7.51 39 Intron (ENSMU… 1 10472475 1.05e7 5349 2 329093 ENSMUST0000… 53902
3 chr1 1.26e7 1.26e7 105 * callpe… 724 . 33.6 79.6 72.5 52 Distal Interg… 1 12667563 1.27e7 5528 1 329096 ENSMUST0000… -22700
4 chr1 1.68e7 1.68e7 64 * callpe… 27 . 5.18 6.92 2.72 44 Distal Interg… 1 16688508 1.67e7 21030 1 17087 ENSMUST0000… 144271
5 chr1 1.89e7 1.89e7 68 * callpe… 69 . 7.55 11.5 6.99 24 Distal Interg… 1 19103022 1.92e7 63325 1 226896 ENSMUST0000… -170251
6 chr1 2.07e7 2.07e7 64 * callpe… 33 . 4.93 7.67 3.36 40 Distal Interg… 1 20730905 2.07e7 3592 1 16171 ENSMUST0000… -22733

### 挑选出promoter区域的peak，和保存peak坐标

srf_15fbs.peak.promoter.coor <- as_tibble(srf_15fbs.peak.anno) %>%
filter(str_detect(annotation, "Promoter")) %>%
select( seqnames, start, end)

srf_03fbs.peak.promoter.coor <- as_tibble(srf_03fbs.peak.anno) %>%
filter(str_detect(annotation, "Promoter")) %>%
select( seqnames, start, end)

write_tsv(srf_03fbs.peak.promoter.coor, col_names =F, file = "ref/srf_03_promoter.coor.tsv")
write_tsv(srf_15fbs.peak.promoter.coor, col_names =F, file = "ref/srf_15_promoter.coor.tsv")

Motif 富集

motif 使用meme-suite, meme-suite推出的不久的R包memes，并且也有一些新功能。这里就使用meme-suite的R包。注意，meme-suite的R包需要在R 4.1版本才能使用。

软件安装：

if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")

BiocManager::install("memes")

数据准备

使用meme-suite 做Motif富集，需要meme格式的motif数据,和peak fasta格式的序列。

meme格式的motif数据库可以在https://meme-suite.org/meme/doc/download.html 这个页面下载。

下载解压后

cd ref
wget https://meme-suite.org/meme/meme-software/Databases/motifs/motif_databases.12.21.tgz
tar xvf motif_databases.12.21.tgz
ll motif_databases/ | cut -d " " -f 5-
83 Jul 7 13:37 ARABD
16384 Jul 7 13:37 CIS-BP_1.02
32768 Jul 7 13:37 CIS-BP_2.00
28672 Jul 7 13:37 CISBP-RNA
73 Jul 7 13:37 ECOLI
4096 Jul 7 13:37 EUKARYOTE
166 Jul 7 13:37 FLY
200 Jul 7 13:37 HUMAN
4096 Jul 7 13:37 JASPAR
51 Jul 7 13:37 MALARIA
34 Jul 7 13:37 METHYLCYTOSINE
24576 Jul 7 13:37 MIRBASE
2277914 Dec 15 2020 motif_db.csv
229 Jul 7 13:37 MOUSE
89 Jul 7 13:37 PROKARYOTE
58 Jul 7 13:37 PROTEIN
4096 Jul 7 13:37 RNA
78 Jul 7 13:37 TFBSshape
40 Jul 7 13:37 WORM
178 Mar 28 2020 YEAST

ll motif_databases/MOUSE/ | cut -d " " -f 5-
13294 Dec 15 2020 chen2008.meme
389252 Dec 15 2020 HOCOMOCOv10_MOUSE_mono_meme_format.meme
312535 Dec 15 2020 HOCOMOCOv11_core_MOUSE_mono_meme_format.meme
452711 Dec 15 2020 HOCOMOCOv11_full_MOUSE_mono_meme_format.meme
543129 Dec 15 2020 uniprobe_mouse.meme

peak fasta序列提取。

## srf_03_promoter.coor.tsv 是之前在R里提取的peak坐标信息
## mm10.fa 是小鼠基因组序列
awk -v OFS='\t' '{printf "%s:%s-%s\n",$1,$2,$3}' srf_03_promoter.coor.tsv | while read id; do samtools faidx ref/mm10.fa $id | cat >> ref/srf_03.promoter.fa ;done
awk -v OFS='\t' '{printf "%s:%s-%s\n",$1,$2,$3}' srf_15_promoter.coor.tsv | while read id; do samtools faidx ref/mm10.fa $id | cat >> ref/srf_15.promoter.fa ;done

AME富集

我们采用meme-suite里的ame 富集方法来进行富集。以下是在R 4.1环境下运行的

library(memes)
library(Biostrings)

## 读取序列
srf_03_promoter.fa <- readBStringSet("ref/srf_03.promoter.fa")
srf_15_promoter.fa <- readBStringSet("ref/srf_15.promoter.fa")

## 设置motif db
meme_db_path <- "ref/motif_databases/MOUSE/HOCOMOCOv11_full_MOUSE_mono_meme_format.meme"
options(meme_db = meme_db_path)

srf.fa <- list("FBS_03" = srf_03_promoter.fa, "FBS_15" = srf_15_promoter.fa)
srf_promoter <- runAme(srf.fa)

## 返回的富集结果是个List

## 查看详情
View(srf_promoter$FBS_03)

head(srf_promoter$FBS_03)
rank motif_db motif_id motif_alt_id consensus pvalue adj.pvalue evalue tests fasta_max pos neg pwm_min tp
<int> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <int> <dbl> <int> <int> <dbl> <int>
1 1 /home/… SRF_MOUS… NA TTBCCWTA… 3.47e-33 4.65e-31 2.47e-28 134 138 138 1104 49 36
2 2 /home/ … ELF1_MOU… NA SRACCCGG… 2.91e-15 2.82e-13 1.5 e-10 97 138 138 1104 21.1 20
3 3 /home/ … PITX2_MO… NA TKGGATTA… 2.38e-14 3.71e-12 1.97e- 9 156 138 138 1104 33.2 14
4 4 /home/ … GABPA_MO… NA GGVRCCGG… 1.74e-13 2.08e-11 1.11e- 8 120 138 138 1104 41.1 16
5 5 /home/ … ELK1_MOU… NA RCCGGAAG… 1.74e-13 2.61e-11 1.38e- 8 150 138 138 1104 37.9 16
6 6 /home/ … FLI1_MOU… NA GGVRCCGG… 9.3 e-12 1.7 e- 9 9.03e- 7 183 138 138 1104 6.86 30

## ame 富集heatmap
srf_promoter %>%
dplyr::bind_rows(.id = "condition") %>%
dplyr::group_by(condition) %>%
dplyr:: slice_head(n = 25) %>% ## 选取了前25个motif
ungroup() %>%
plot_ame_heatmap(group = condition)

从图可以看出， SRF对应motif 在15%s FBS条件下相比3% FBS更为显著富集。

Peak 可视化

可以使用功能强大的IGV， 不过这里使用叫Gviz的R 包。使用IGV,使用倒是方便，但是保存图片，我只知道第三方截图软件，或者igv里自带的截图，只能保存成png。以下子在R 4.0 运行跑的。

BiocManager::install("Gviz")
BiocManager::install("rtracklayer")

library(Gviz)
library(rtracklayer)

### 这里用的数据，高峰值的peak很少，这里挑选了一个比较大的来可视化
chr_no <- "1"
chr_start <- 12636629
chr_end <- 12652990

fbs_03 <- import.bedGraph("callpeak/SRF_03FBS_treat_pileup.bdg") ## bdg 是macs3运行生成的
fbs_15 <- import.bedGraph("callpeak/SRF_15FBS_treat_pileup.bdg")

fbs_03_track <- DataTrack(range=fbs_03,
chromosome=chr_no,
ylim=c(0, 50),
name = "3% FBS",
background.title="lightgray"
)

fbs_15_track <- DataTrack(range=fbs_15,
chromosome=chr_no,
ylim=c(0, 50),
name = "15% FBS",
background.title="lightgray"
)

plotTracks(list(fbs_03_track, fbs_15_track),
from =chr_start,
to=chr_end,
type = "h")

参考

https://b23.tv/rBJhUh

https://mp.weixin.qq.com/s/_OPzvaEAbiMolCA2mqJXLw

https://master.bioconductor.org/packages/release/bioc/vignettes/memes/inst/doc/core_ame.html

https://jnicolaus.com/tutorials/gviztutorial.html

ChIP-Seq数据分析上下游打通

数据

比对